3d 灰度体积投影到 2D 平面

Question

我有一个对应于超声数据的 3-D 灰度体积。在 Matlab 中，这个 3-D 体积只是 MxNxP 的 3-D 矩阵。我感兴趣的结构不是沿z 轴，而是沿已知的局部坐标系（x'y'z'）。到目前为止，我所拥有的类似于下图，描绘了原始 (xyz) 和局部坐标系 (x'y'z')：

我想通过局部坐标系上的特定平面获得该体积（即图像）的二维投影，例如 z' = z0。我该怎么做？

如果体积沿 z 轴定向，则可以轻松实现此投影。即，如果在 Matlab 中的体积是 V，那么：

projection = sum(V,3);

因此，可以将投影计算为沿数组第三维的总和。然而，随着方向的改变，问题变得更加复杂。

我一直在研究 radon 变换（2D，仅适用于 2-D 图像而不是体积）并且还在考虑正交投影，但此时我不知道该怎么做！

感谢您的建议！

Answer 1

解决方案的新尝试：

按照教程http://blogs.mathworks.com/steve/2006/08/17/spatial-transformations-three-dimensional-rotation/ 并进行一些小的更改，我可能会有一些可以帮助您的东西。请记住，我对 MATLAB 中的体积数据几乎没有经验，因此实现起来非常麻烦。

在下面的代码中，我使用 tformarray() 在空间中旋转结构。首先，将数据居中，然后使用 rotationmat3D 进行旋转以产生空间变换，然后再将数据移回其原始位置。

由于我以前从未使用过 tformarray，因此我通过简单地用零填充数据矩阵 (NxMxP) 来处理旋转后落在定义区域之外的数据点。如果有人知道更好的方法，请告诉我们:)

代码：

    %Synthetic dataset, 25x50x25
blob = flow();

%Pad to allow for rotations in space. Bad solution, 
%something better might be possible to better understanding
%of tformarray()
blob = padarray(blob,size(blob));

f1 = figure(1);clf;
s1=subplot(1,2,1);
p = patch(isosurface(blob,1));
set(p, 'FaceColor', 'red', 'EdgeColor', 'none');
daspect([1 1 1]);
view([1 1 1])
camlight
lighting gouraud

%Calculate center
blob_center = (size(blob) + 1) / 2;

%Translate to origin transformation
T1 = [1 0 0 0
    0 1 0 0
    0 0 1 0
    -blob_center 1];

%Rotation around [0 0 1]
rot = -pi/3;
Rot = rotationmat3D(rot,[0 1 1]);
T2 = [ 1  0  0   0
       0  1  0   0
       0  0  1   0
       0  0  0   1];
T2(1:3,1:3) = Rot;   

%Translation back
T3 = [1 0 0 0
    0 1 0 0
    0 0 1 0
    blob_center 1];

%Total transform
T = T1 * T2 * T3;

%See http://blogs.mathworks.com/steve/2006/08/17/spatial-transformations-three-dimensional-rotation/
tform = maketform('affine', T);
R = makeresampler('linear', 'fill');
TDIMS_A = [1 2 3];
TDIMS_B = [1 2 3];
TSIZE_B = size(blob);
TMAP_B = [];
F = 0;
blob2 = ...
tformarray(blob, tform, R, TDIMS_A, TDIMS_B, TSIZE_B, TMAP_B, F);

s2=subplot(1,2,2);
p2 = patch(isosurface(blob2,1));
set(p2, 'FaceColor', 'red', 'EdgeColor', 'none');
daspect([1 1 1]);
view([1 1 1])
camlight
lighting gouraud

下面的任意可视化只是为了确认数据按预期旋转，当数据传递值“​​1”时绘制一个闭合曲面。使用 blob2，您应该知道能够使用简单的求和进行投影。

figure(2)
subplot(1,2,1);imagesc(sum(blob,3));
subplot(1,2,2);imagesc(sum(blob2,3));

Answer 2

假设您可以访问坐标基础 R=[x' y' z']，并且这些向量是正交的，您可以通过将数据与 3x3 矩阵 R 相乘来简单地提取此基础中的表示，其中x',y',z' 是列向量。

使用存储在 D (Nx3) 中的数据，您可以通过乘以 R 来获得表示： Dmarked = D*R；

现在 D = Dmarked*inv(R)，所以来回是顺势而为。

以下代码可能有助于查看转换。在这里，我创建了一个合成数据集，将其旋转，然后将其旋转回来。做 sum(DR(:,3)) 将是你沿 z' 的总和

%#Create synthetic dataset
N1 = 250;
r1 = 1;
dr1 = 0.1;
dz1 = 0;
mu1 = [0;0];
Sigma1 = eye(2);
theta1 = 0 + (2*pi).*rand(N1,1);
rRand1 = normrnd(r1,dr1,1,N1);
rZ1 = rand(N1,1)*dz1+1;
D = [([rZ1*0 rZ1*0] + repmat(rRand1',1,2)).*[sin(theta1) cos(theta1)] rZ1];

%Create roation matrix
rot = pi/8;
R = rotationmat3D(rot,[0 1 0]);
% R =       0.9239          0       0.3827
%           0               1.0000  0
%           -0.3827         0       0.9239

Rinv = inv(R);

%Rotate data
DR = D*R;

%#Visaulize data
f1 = figure(1);clf
subplot(1,3,1);
plot3(DR(:,1),DR(:,2),DR(:,3),'.');title('Your data')
subplot(1,3,2);
plot3(DR*Rinv(:,1),DR*Rinv(:,2),DR*Rinv(:,3),'.r');
view([0.5 0.5 0.2]);title('Representation using your [xmarked ymarked zmarked]');
subplot(1,3,3);
plot3(D(:,1),D(:,2),D(:,3),'.');
view([0.5 0.5 0.2]);title('Original data before rotation');

Answer 3

如果您有两个 标准化 3x1 向量 x2 和 y2 对应于您的本地坐标系（x' 和 y'）。

那么，对于位置P，其局部坐标将为xP=P'x2和yP=P'*y2。

所以您可以尝试使用 accumarray 投影您的音量：

[x y z]=ndgrid(1:M,1:N,1:P);
posP=[x(:) y(:) z(:)];
xP=round(posP*x2);
yP=round(posP*y2);
xP=xP+min(xP(:))+1;
yP=yP+min(yP(:))+1;
V2=accumarray([xP(:),yP(:)],V(:));

如果你提供你的数据，我会测试它。