create cylinders in 3D volumetric data

Question

I'm trying to create a dataset of raw volumetric data consisting of geometrical shapes. The point is to use volume ray casting to project them in 2D but first I want to create the volume manually.

The geometry is consisting of one cylinder that is in the middle of the volume, along the Z axis and 2 smaller cylinders that are around the first one, deriving from rotations around the axes.

Here is my function so far:

function cyl= createCylinders(a, b, c, rad1, h1, rad2, h2)

% a : data width
% b : data height
% c : data depth
% rad1: radius of the big center cylinder
% rad2: radius of the smaller cylinders
% h1: height of the big center cylinder
% h2: height of the smaller cylinders

[Y X Z] =meshgrid(1:a,1:b,1:c);  %matlab saves in a different order so X must be Y
centerX = a/2;
centerY = b/2;
centerZ = c/2;

theta = 0; %around y
fi = pi/4; %around x

% First cylinder
cyl = zeros(a,b,c);

% create for infinite height
R = sqrt((X-centerX).^2 + (Y-centerY).^2);

startZ = ceil(c/2) - floor(h1/2);
endZ = startZ + h1 - 1;

% then trim it to height = h1
temp = zeros(a,b,h1);
temp( R(:,:,startZ:endZ)<rad1 ) = 255;
cyl(:,:,startZ:endZ) = temp;

% Second cylinder

cyl2 = zeros(a,b,c);

A = (X-centerX)*cos(theta) + (Y-centerY)*sin(theta)*sin(fi) + (Z-centerZ)*cos(fi)*sin(theta);
B = (Y-centerY)*cos(fi) - (Z-centerZ)*sin(fi);

% create again for infinite height
R2 = sqrt(A.^2+B.^2);
cyl2(R2<rad2) = 255;

%then use 2 planes to trim outside of the limits
N = [ cos(fi)*sin(theta) -sin(fi) cos(fi)*cos(theta) ];

P = (rad2).*N + [ centerX centerY centerZ];
T = (X-P(1))*N(1) + (Y-P(2))*N(2) + (Z-P(3))*N(3);
cyl2(T<0) = 0;

P = (rad2+h2).*N + [ centerX centerY centerZ];
T = (X-P(1))*N(1) + (Y-P(2))*N(2) + (Z-P(3))*N(3);
cyl2(T>0) = 0;

% Third cylinder
% ...

cyl = cyl + cyl2;

cyl = uint8(round(cyl));

% ...

The concept is that the first cylinder is created and then "cut" according to the z-axis value, to define its height. The other cylinder is created using the relation A² + B ² = R² where A and B are rotated accordingly using the rotation matrices only around x and y axes, using R_y(θ)R_x(φ) as described here.

Until now everything seems to be working, because I have implemented code (tested that it works well) to display the projection and the cylinders seem to have correct rotation when they are not "trimmed" from infinite height.

I calculate N which is the vector [0 0 1] aka z-axis rotated in the same way as the cylinder. Then I find two points P of the same distances that I want the cylinder's edges to be and calculate the plane equations T according to that points and normal vector. Lastly, I trim according to that equality. Or at least that's what I think I'm doing, because after the trimming I usually don't get anything (every value is zero). Or, the best thing I could get when I was experimenting was cylinders trimmed, but the planes of the top and bottom where not oriented well.

I would appreciate any help or corrections at my code, because I've been looking at the geometry equations and I can't find where the mistake is.

Edit: This is a quick screenshot of the object I'm trying to create. NOTE that the cylinders are opaque in the volume data, all the inside is considered as homogeneous material.

Answer 1

I think instead of:

T = (X-P(1))*N(1) + (Y-P(2))*N(2) + (Z-P(3))*N(3);

you should try the following at both places:

T = (X-P(1)) + (Y-P(2)) + (Z-P(3));

Multiplying by N is to account for the direction of the axis of the 2nd cylinder which you have already done just above that step.