在3D容积数据中创建柱面

Question

我正在尝试创建由几何形状组成的原始体积数据的数据集。重点是使用体积射线投射以2D方式投射它们，但首先我要手动创建体积。

几何体由位于体积中间的一个圆柱体沿着Z轴和围绕第一个圆柱体的两个较小圆柱体组成，这些圆柱体围绕轴线旋转。

这里是到目前为止我的功能：

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)); 

% ... 

的概念是，第一汽缸被创建，然后“切割”根据z轴值，以限定其高度。使用关系创建其他缸内A + B = R 其中A和B相应地旋转利用旋转矩阵只围绕x和y轴，，使用R _Ý（θ） R _x（φ）如here所述。

到现在为止，所有东西似乎都在工作，因为我已经实现了代码（测试它很好地工作）来显示投影，并且当它们没有从无限高度“修剪”时，柱面似乎有正确的旋转。

我计算N这是矢量[0 0 1]又名Z轴旋转方式与圆柱相同。然后我找到两个相同距离的点P，我想要圆柱体的边缘，并根据这些点和法向量计算平面方程T。最后，我根据这种平等来调整。或者至少这就是我所认为的我在做，因为修剪后我通常不会得到任何东西（每个值都为零）。或者，当我试验时，我能得到的最好的东西是气瓶修剪，但顶部和底部的飞机没有定向好。

我将不胜感激我的代码任何帮助或更正，因为我一直在寻找的几何公式​​，我找不到在哪里的错误是。

编辑： 这是我试图创建的对象的快速截图。请注意，体积数据中的圆柱体是不透明的，所有内部都被视为均质材料。

Answer 1

我认为不是的：

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

你应该尝试在这两个地方如下：

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

乘以N是考虑到轴的方向第二个气缸，你已经完成了该步骤。