create a 3d cylinder inside 3d volume

Question

I have 3d volume. Which has shape of (399 x 512 x 512). And It has voxel spacing of 0.484704 x 0.484704 x 0.4847

Now, I want to define a cylinder inside this volume with length 5mm, diameter 1mm, intensity 1 inside, intensity 0 outside.

I saw an example to define a cylinder in internet like this code:

from mpl_toolkits.mplot3d import Axes3D  
def data_for_cylinder_along_z(center_x,center_y,radius,height_z):
    z = np.linspace(0, height_z, 50)
    theta = np.linspace(0, 2*np.pi, 50)
    theta_grid, z_grid=np.meshgrid(theta, z)
    x_grid = radius*np.cos(theta_grid) + center_x
    y_grid = radius*np.sin(theta_grid) + center_y
    return x_grid,y_grid,z_grid

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

Xc,Yc,Zc = data_for_cylinder_along_z(0.2,0.2,0.05,0.1)
ax.plot_surface(Xc, Yc, Zc, alpha=0.5)

plt.show()

However, I don't know how to define the cylinder inside the 3d volume keeping all the conditions (length 5mm, diameter 1mm, intensity 1 inside, intensity 0 outside) true. I also want to define the center of cylinder automatically. So that I can define the cylinder at any place of inside the 3d volume keeping the other condition true. Can anyone show or provide any example? Thanks a lot in advance.

score 0 · Accepted Answer · answered Nov 11 '20 at 13:31

One simple way of solving this would be to perform each of the checks individually and then just keep the voxels that satisfy all of your constraints.

If you build a grid with all of the centers of the voxels: P (399 x 512 x 512 x 3), each voxel at (i,j,k) will be associated with its real-world position (x,y,z).

That's a little tricky, but it should look something like this:

np.stack(np.meshgrid(np.arange(0, shape[0]),
                     np.arange(0, shape[1]),
                     np.arange(0, shape[2]), indexing='ij'), axis=3)

If you subtract the cylinder's center (center_x,center_y, center_z), you're left with the relative positions of each (i,j,k) voxel P_rel (399 x 512 x 512 x 3)

When you have that, you can apply each of your tests one after the other. For a Z-oriented cylinder with a radius and height_z it would look something like:

# constrain the Z-axis
not_too_high = P_rel[:,:,:,2]<= (0.5*height_z)
not_too_low = P_rel[:,:,:,2]>= (-0.5*height_z)

# constrain the radial direction
not_too_far = np.linalg.norm(P_rel[:,:,:,:2],axis=3)<=radius

voxels_in_cyl = not_too_high & not_too_low & not_too_far

I haven't tested the code, but you get the idea.

If you wanted to have an cylinder with an arbitrary orientation you would have to project P_rel into axial and radial components and then do an analogous check without "hard-coding" the indices as I did in this example

I actually worked with the idea not the code snippet. but it works :) — Tokai, Nov 13 '20 at 09:48
Hi, tested actual code snippet. When print `voxels_in_cyl` I get boolean values! Is it possible to give a complete example to define a cylinder, given a center, height and radius ? Would be very helpful. — Tokai, Dec 10 '20 at 21:50
The voxels inside the cylinder are 1 and the ones outside are 0 — that's sounds like a "complete" definition to me! You can use that as a mask to get the coordinates of the voxels inside (`P_rel[voxels_in_cyl]`) or convert it to whatever your use-case requires — imochoa, Dec 11 '20 at 10:05

create a 3d cylinder inside 3d volume

1 Answers1