How to match values of different meshgrids?

Question

I am trying to create a grayscale stack of images of a cylinder starting from the information about the radius, size and orientation of a cylinder.

This cylinder should be contained in a cubic 3D mesh, where each gridpoint will represent a pixel.

The code would be something like this.

  #Define meshgrid
  x_ = np.linspace(0,Lx,int(Lx/vox))
  y_ = np.linspace(0,Ly,int(Ly/vox))
  z_ = np.linspace(0,Lz,int(Lz/vox))

  X,Y,Z = np.meshgrid(x_,y_,z_,indexing='ij')

Where vox is the voxel size I input. This would define the 3D mesh, the Matrix from which we will derive the grayscale stack will be defined as

  M = np.zeros((x_.size,y_.size,z_.size),dtype=int)

Now we are given a Cylinder whose center axis goes from point p1 to p2, and has a radius r. So based on the following link Numpy mask from cylinder coordinates, I create an additional meshgrid that contains all points belonging to such cylinder.

        #Vector
        v = p2 - p1

        #Normalize vector
        lenght = scipy.linalg.norm(v)
        v = v / lenght

        # make some vector not in the same direction as v
        not_v = np.array([1.0, 0, 0])
        if (v == not_v).all():
              not_v = np.array([0, 1.0, 0])
        # make vector perpendicular to v
        n1 = np.cross(v, not_v)
        # normalize n1
        n1 = n1 / scipy.linalg.norm(n1)
        # make unit vector perpendicular to v and n1
        n2 = np.cross(v, n1)

        #Define gridpoints for the cilinder
        l_ =      np.linspace(0,lenght,100)
        r_ =      np.linspace(0,r,10)
        theeta_ = np.linspace(0,2*np.pi,10)

        #define meshgrid for cilinder
        L,R,Theeta = np.meshgrid(l_,r_,theeta_,indexing='ij')

Then I would get the xyz coordinates from these cylindrical coordinates.

#Transform to x, y, z coordinates
Xc, Yc, Zc = [p1[i] + v[i] * L + R * np.sin(Theeta) * n1[i] + r * np.cos(Theeta) * n2[i] for i in [0, 1, 2]]

So now the question is the following. I have defined the gridpoints that constitute the given cylinder in cartesian coordinates. However, the problem I have now is with trying to match this meshgrid for the cylinder with the meshgrid defined to obtain the grayscale stack of images.

So what I would like to do is to take this Xc, Yc, Zc, coordinates and see to which X,Y,Z coordinates would they correspond, and light up the corresponding pixels in the matrix M. But I cannot see an obvious way of doing it.

Regards,

Answer 1

You are approaching the problem the wrong way, what you should do is directly find which points in your coordinate grid fall within your cylinder. This is how you can do that:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Problem parameters
Lx, Ly, Lz = 50, 20, 30
vox = 2.0
p1 = np.array([12., 8., 4.])
p2 = np.array([37., 14., 22.])
r = 5.
#Define meshgrid
x_ = np.linspace(0, Lx, int(Lx / vox))
y_ = np.linspace(0, Ly, int(Ly / vox))
z_ = np.linspace(0, Lz, int(Lz / vox))
# Grid coordinates
X, Y, Z = np.meshgrid(x_, y_, z_, indexing='ij')
# Stack into an array
coords = np.stack([X, Y, Z], axis=-1)
# Compute distance from each point to cylinder axis
v = p2 - p1
t = np.dot(coords - p1, v) / np.dot(v, v)
p = p1 + np.expand_dims(t, axis=-1) * v
dist = np.linalg.norm(coords - p, axis=-1)
# Select points within cylinder distance and bounds
mask = (dist <= r) & (0 <= t) & (t <= 1)
print(mask.shape)
# (10, 4, 6)
# Select coordinates in cylinder
cyl_coords = coords[mask]
# Plot
ax = plt.figure().add_subplot(111, projection='3d')
ax.scatter3D(cyl_coords[:, 0], cyl_coords[:, 1], cyl_coords[:, 2], s=3)
# Set plot limits for uniform aspect ratio
ax.set_xlim(0, 50)
ax.set_ylim(-15, 35)
ax.set_zlim(-10, 40)
ax.figure.tight_layout()
ax.figure.show()

Output: