c++ - 如何沿一定体积绘制斜半球？

Question

我想做什么:

制作一个空的 3D 图像(在本例中为 .dcm)，图像方向为[1,0,0;
0,1,0；
0,0,1]。在这张图片中，我插入了一个倾斜的轨迹，它本质上代表了一个长方体。现在我想在这个长方体中插入一个空心半球(所有白色像素的长方体 - 常数值，半球可以是任何东西，但不可微)，以便它沿着轨迹的轴对齐。

我得到了什么

所以我使用了球体的通用公式:

x = x0 + r*cos(theta)*sin(alpha)
y = y0 + r*sin(theta)*sin(alpha)
z = z0 + r*cos(alpha)

其中，对于半球，0 <= theta <= 2 * pi，0 <= alpha <= pi/2。

我为实现这一目标所做的努力

1. 所以我首先想到的是获取图像坐标系和轨迹坐标系之间的旋转矩阵，然后将球体上的所有点与其相乘。由于缩放和平移了旋转的球体，这并没有给我想要的结果。当我自己检查这些点时，我不明白为什么会这样。

2. 然后我想为什么不用一个球体做一个半球体，这个球体被一个平行于轨迹坐标系的 y、z 平面的平面切割。为此，我计算了图像的 x、y 和 z 轴与轨迹的夹角。然后，我开始获取 theta_rotated 和 alpha_rotated 的半球坐标。这也不起作用，因为我得到的不是半球，而是一个相当奇怪的球体。

这是没有任何转换

这是角度变换(第二次尝试)

供引用，

轨迹坐标系:

[-0.4744，-0.0358506，-0.8553；
-0.7049、0.613244、0.3892；
-0.5273, -0.787537, 0.342;];

给出角度:

x_axis 角度 2.06508 piy_axis 角度 2.2319 piz_axis 角度 1.22175 pi

生成长方体的代码

Vector3d getTrajectoryPoints(std::vector<Vector3d> &trajectoryPoints, Vector3d &target1, Vector3d &tangent1){
    double distanceFromTarget = 10;
    int targetShift = 4;
    target -= z_vector;
    target -= (tangent * targetShift);
    Vector3d vector_x = -tangent;
    y_vector = z_vector.cross(vector_x);
    target -= y_vector;
    Vector3d start = target - vector_x * distanceFromTarget;

    std::cout << "target = " << target << "start = " << start << std::endl;
    std::cout << "x " << vector_x << " y " << y_vector << " z " << z_vector << std::endl;
    double height = 0.4;
    while (height <= 1.6)
    {
        double width = 0.4;
        while (width <= 1.6){
            distanceFromTarget = 10;
            while (distanceFromTarget >= 0){
                Vector3d point = target + tangent * distanceFromTarget;
                //std::cout << (point + (z_vector*height) - (y_vector * width)) << std::endl;
                trajectoryPoints.push_back(point + (z_vector * height) + (y_vector * width));
                distanceFromTarget -= 0.09;
            }
            width += 0.09;
        }
        height += 0.09;
    }
}

相对于体素间距递增的高度和宽度。

你们知道如何实现这一目标吗？我做错了什么？如果您需要任何其他信息，请告诉我。

编辑 1在@Dzenan 回答后，我尝试了以下操作:

目标 = { -14.0783, -109.8260, -136.2490 }, 切线 = { 0.4744, 0.7049, 0.5273 };

typedef itk::Euler3DTransform<double> TransformType;
    TransformType::Pointer transform = TransformType::New();

    double centerForTransformation[3];

const double pi = std::acos(-1);
try{
        transform->SetRotation(2.0658*pi, 1.22175*pi, 2.2319*pi);
        // transform->SetMatrix(transformMatrix);
    }
    catch (itk::ExceptionObject &excp){
        std::cout << "Exception caught ! " << excp << std::endl;
        transform->SetIdentity();
    }
transform->SetCenter(centerForTransformation);

然后我循环遍历半球中的所有点并使用以下方法转换它们，

point = transform->TransformPoint(point);

虽然，我更愿意给出等于轨迹坐标系(如上所述)的矩阵，但该矩阵不是正交的，它不会接受它。必须要说的是，我使用相同的矩阵对该图像进行重采样并提取长方体，这很好。因此，我找到了 x_image - x_trajectory、y_image - y_trajectory 和 z_image - z_trajectory 之间的角度，并使用 SetRotation 代替，这给了我以下结果(仍然不正确):

编辑 2

我试图在不实际使用极坐标的情况下获取球体坐标。在与@jodag 讨论之后，这就是我想出的:

Vector3d center = { -14.0783, -109.8260, -136.2490 };

    height = 0.4;
    while (height <= 1.6)
    {
        double width = 0.4;
        while (width <= 1.6){
            distanceFromTarget = 5;
            while (distanceFromTarget >= 0){
// Make sure the point lies along the cuboid direction vectors
                Vector3d point = center + tangent * distanceFromTarget + (z_vector * height) + (y_vector * width);
                double x = std::sqrt((point[0] - center[0]) * (point[0] - center[0]) + (point[1] - center[1]) * (point[1] - center[1]) + (point[2] - center[2]) * (point[2] - center[2]));
                if ((x <= 0.5) && (point[2] >= -136.2490 ))
                    orientation.push_back(point);
                distanceFromTarget -= 0.09;
            }
            width += 0.09;
        }
        height += 0.09;
    }

但这似乎也行不通。

这是输出

Answer 1

我对你的第一个图有点困惑，因为显示的点似乎没有在图像坐标中定义。我在下面发布的示例假设体素必须是图像坐标系的一部分。

下面的代码使用逆变换将图像空间中的体素坐标转换为轨迹空间。然后栅格化一个以 0,0,0 为中心的 2x2x2 立方体和一个沿 xy 轴切片的半径为 0.9 的半球。

与其在评论中继续进行长时间的讨论，我决定发布这个。如果您正在寻找不同的东西，请发表评论。

% define trajectory coordinate matrix
R = [-0.4744, -0.0358506, -0.8553;
     -0.7049, 0.613244, 0.3892;
     -0.5273, -0.787537, 0.342]

% initialize 50x50x50 3d image
[x,y,z] = meshgrid(linspace(-2,2,50));
sz = size(x);
x = reshape(x,1,[]);
y = reshape(y,1,[]);
z = reshape(z,1,[]);
r = ones(size(x));
g = ones(size(x));
b = ones(size(x));

blue = [0,1,0];
green = [0,0,1];

% transform image coordinates to trajectory coordinates
vtraj = R\[x;y;z];
xtraj = vtraj(1,:);
ytraj = vtraj(2,:);
ztraj = vtraj(3,:);

% rasterize 2x2x2 cube in trajectory coordinates
idx = (xtraj <= 1 & xtraj >= -1 & ytraj <= 1 & ytraj >= -1 & ztraj <= 1 & ztraj >= -1);
r(idx) = blue(1);
g(idx) = blue(2);
b(idx) = blue(3);

% rasterize radius 0.9 hemisphere in trajectory coordinates
idx = (sqrt(xtraj.^2 + ytraj.^2 + ztraj.^2) <= 0.9) & (ztraj >= 0);
r(idx) = green(1);
g(idx) = green(2);
b(idx) = green(3);

% plot just the blue and green voxels
green_idx = (r == green(1) & g == green(2) & b == green(3));
blue_idx = (r == blue(1) & g == blue(2) & b == blue(3));

figure(1); clf(1);
plot3(x(green_idx),y(green_idx),z(green_idx),' *g')
hold('on');
plot3(x(blue_idx),y(blue_idx),z(blue_idx),' *b')

axis([1,100,1,100,1,100]);
axis('equal');
axis('vis3d');