3D灰度体投影到2D平面上

Question

3D灰度体投影到2D平面上

4

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

[image]

我想通过本地坐标系上的特定平面获得该体积的2D投影（即图像），比如在z'=z0处。我该如何做？

如果该体积沿z轴定向，这个投影可以很容易地实现。也就是说，如果该体积在Matlab中被表示为V，则：

projection = sum(V,3);

因此，投影可以通过计算数组沿第三维的总和来计算。但是，如果方向改变，问题就变得更加复杂。

我一直在研究Radon变换（2D，仅适用于2D图像而不是体积），并且也考虑了正交投影，但是目前我对该怎么做毫无头绪！

感谢任何帮助！

- CodificandoBits

到目前为止，我正在尝试使用三维仿射变换的方法，参考这篇博客（Matlab的图像处理专家Steve）。然而，我仍然缺少一个部分：如何根据我的正交基向量x' y' z'来定义旋转矩阵中的角度？ - CodificandoBits

3个回答

1

假设您可以访问坐标基础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');

enter image description here

- Vidar

我有几个问题：(i) 为什么数据向量_D_是Nx3？在我的情况下，数据指的是每个体素的灰度级别，因此对于MxNxP的体积，我有M.N.P个数据点，而不是Nx3；(ii) 数据在基础上的表示（即Dmarked = D*R）实际上并不是将D的点坐标（x，y，z）表示为基础[x' y' z']的表示吗？(iii) "rotationmat3D"函数代表什么？我在哪里可以找到它？...谢谢！ - CodificandoBits

嗯，我想我没有理解你的数据结构。你正在使用等间距空间样本进行工作，宽度为M点，深度为N点，高度为P点？Rotmat3D可以在这里找到：http://www.mathworks.com/matlabcentral/fileexchange/23417-compute-3d-rotation-matrix - Vidar

是的，我正在使用均匀间隔的样本进行工作：宽度为M个体素，深度为N，高度为P。 - CodificandoBits

请看下面的新回答。我应该删除之前的吗？ - Vidar

0

如果您有两个归一化的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(:));

如果您提供数据，我会进行测试。

- Oli

在这种情况下的问题是，在获取本地坐标xP和yP时，有时这些下标超出了边界（即为负数、零或不包含在体积中）。例如：V = randn(200,100,50); %合成体积
[M,N,P] = size(V);
[x,y,z] = meshgrid(1:N,1:M,1:P);
x2 = [0.7071 0.7071 -0.0035]'; %在x'方向上的基向量
y2 = [-0.6869 0.6857 -0.2410]'; %在y'方向上的基向量 使用这些数据时，yP值超出了限制（超出了我的200x100x50体积的范围）。 - CodificandoBits

我添加了以下代码行以处理小于零的值：xP=xP+min(xP(:))+1;yP=yP+min(yP(:))+1;。如果你真的需要边界，可以这样做：inds= xP>lowX & xP<highX & yP>lowY & yP<highY; xP=xP(inds);yP=yP(inds);V=V(inds); - Oli

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Vidar · Accepted Answer

新的解决方案尝试：

参考教程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));

enter image description here