在3D图像上进行PCA以获得3个主轴

Question

我有一个解剖体积图像(B)，它是一个索引图像i，j，k：

B(1,1,1)=0 %for example.

该文件仅包含0和1。

我可以用isosurface正确地将其可视化：isosurface(B);

我想在j中的某个坐标上切割体积(每个体积是不同的)。

问题是体积是垂直倾斜的，它可能有45%的角度，所以切割不会跟随解剖体积。

我想为数据获得一个新的正交坐标系，这样我在坐标j中的平面将以更精确的方式切割解剖体积。

有人告诉我使用PCA来做这件事，但我不知道如何做，而且阅读帮助页面也没有帮助。任何方向都是受欢迎的！

编辑:我一直在遵循建议，现在我有了一个新的卷，以零为中心，但我认为轴没有遵循解剖学图像的应有之处。以下是前图和后图：

​

​

​

这是我用来创建图像的代码(我从答案中提取了一些代码，从评论中提取了一些想法)：

clear all; close all; clc;
hippo3d = MRIread('lh_hippo.mgz');
vol = hippo3d.vol;

[I J K] = size(vol);


figure;
isosurface(vol);

% customize view and color-mapping of original volume
daspect([1,1,1])
axis tight vis3d; 
camlight; lighting gouraud
camproj perspective
colormap(flipud(jet(16))); caxis([0 1]); colorbar
xlabel x; ylabel y; zlabel z
box on

% create the 2D data matrix
a = 0;
for i=1:K
    for j=1:J
        for k=1:I
            a = a + 1;
            x(a) = i;
            y(a) = j;
            z(a) = k;
            v(a) = vol(k, j, i);
        end
    end
end

[COEFF SCORE] = princomp([x; y; z; v]');

% check that we get exactly the same image when going back
figure;
atzera = reshape(v, I, J, K);
isosurface(atzera);
% customize view and color-mapping for the check image
daspect([1,1,1])
axis tight vis3d; 
camlight; lighting gouraud
camproj perspective
colormap(flipud(jet(16))); caxis([0 1]); colorbar
xlabel x; ylabel y; zlabel z
box on

% Convert all columns from SCORE
xx = reshape(SCORE(:,1), I, J, K);
yy = reshape(SCORE(:,2), I, J, K);
zz = reshape(SCORE(:,3), I, J, K);
vv = reshape(SCORE(:,4), I, J, K);

% prepare figure
%clf
figure;
set(gcf, 'Renderer','zbuffer')

% render isosurface at level=0.5 as a wire-frame
isoval = 0.5;
[pf,pv] = isosurface(xx, yy, zz, vv, isoval);
p = patch('Faces',pf, 'Vertices',pv, 'FaceColor','none', 'EdgeColor',[0.5 1 0.5]);

% customize view and color-mapping
daspect([1,1,1])
axis tight vis3d;view(-45,35);
camlight; lighting gouraud
camproj perspective

colormap(flipud(jet(16))); caxis([0 1]); colorbar
xlabel x; ylabel y; zlabel z
box on

有没有人能提供一个提示，可能会发生什么？问题可能出在reshape命令上，有没有可能是我取消了之前完成的工作？

Answer 1

不确定PCA，但这里有一个例子展示了如何可视化3D标量体积数据，并在倾斜的平面上切割体积(非轴对齐)。代码的灵感来自MATLAB文档中的this demo。

% volume data
[x,y,z,v] = flow();
vv = double(v < -3.2);  % threshold to get volume with 0/1 values
vv = smooth3(vv);       % smooth data to get nicer visualization :)

xmn = min(x(:)); xmx = max(x(:));
ymn = min(y(:)); ymx = max(y(:));
zmn = min(z(:)); zmx = max(z(:));

% let create a slicing plane at an angle=45 about x-axis,
% get its coordinates, then immediately delete it
n = 50;
h = surface(linspace(xmn,xmx,n), linspace(ymn,ymx,n), zeros(n));
rotate(h, [-1 0 0], -45)
xd = get(h, 'XData'); yd = get(h, 'YData'); zd = get(h, 'ZData');
delete(h)

% prepare figure
clf
set(gcf, 'Renderer','zbuffer')

% render isosurface at level=0.5 as a wire-frame
isoval = 0.5;
[pf,pv] = isosurface(x, y, z, vv, isoval);
p = patch('Faces',pf, 'Vertices',pv, ...
    'FaceColor','none', 'EdgeColor',[0.5 1 0.5]);
isonormals(x, y, z, vv, p)

% draw a slice through the volume at the rotated plane we created
hold on
h = slice(x, y, z, vv, xd, yd, zd);
set(h, 'FaceColor','interp', 'EdgeColor','none')

% draw slices at axis planes
h = slice(x, y, z, vv, xmx, [], []);
set(h, 'FaceColor','interp', 'EdgeColor','none')
h = slice(x, y, z, vv, [], ymx, []);
set(h, 'FaceColor','interp', 'EdgeColor','none')
h = slice(x, y, z, vv, [], [], zmn);
set(h, 'FaceColor','interp', 'EdgeColor','none')

% customize view and color-mapping
daspect([1,1,1])
axis tight vis3d; view(-45,35);
camlight; lighting gouraud
camproj perspective
colormap(flipud(jet(16))); caxis([0 1]); colorbar
xlabel x; ylabel y; zlabel z
box on

下面的结果显示了渲染为线框的等值面，以及轴向对齐和旋转的切片平面。我们可以看到，等值面内部的体积空间的值等于1，而外部的体积空间值为0

​

​

Answer 2

我不认为PCA能解决你的问题。如果你将PCA应用于你的数据，它会给你三个新的轴。这些轴称为主成分(PC)。它们具有这样的特性，即当数据被投影到第一台PC上时，它具有最大的方差。当数据投影到第二个PC上时，第二个PC也必须具有最大的方差，约束条件是它应该与第一个PC正交，第三个PC是类似的。

现在的问题是，当您将数据投影到新的坐标系(由PC定义)时，它是否与解剖体积相匹配？不一定，也很可能不会。数据的右轴没有PCA的优化目标。

Answer 3

抱歉，我试着回答@Tevfik-Aytekin，但答案太长了。

希望这个答案对某些人有用：

嗨@Tevfik，谢谢你的回答。我已经为同样的问题挣扎了好几天，我想我现在已经解决了。

我认为与你所说的不同之处在于，我使用的是坐标。当我在我的坐标上执行PCA时，我得到了我的数据的3x3转换矩阵(COEFF矩阵，它是么正的，它只是一个旋转矩阵)，所以我知道我在转换时得到的体积是完全相同的。

以下是我遵循的步骤：

• I有一个(I，J，K)，3D体积。
• 根据@werner的建议，我把它改成了一个4列矩阵(x，y，z，v)，大小(I*J*K，4)。当v == 0和v时，
• 消除了坐标(x，y，z)。所以现在，大小是(原始卷，3)。只有值为1的坐标，所以向量的长度就是图形的体积。
• 执行主成分分析以获得系数和得分。
• 系数是一个3x3矩阵。它是么正的，它是我的体积数据的旋转矩阵。
• I在PCA1轴上进行了编辑(即，当系数(1)大于某个值时，删除al值)。这是我的问题，我想要切割垂直于主轴的体积，

，

• ，这对我来说已经足够了，旋转坐标给了我想要的体积。但是，我不需要返回，因为我只需要剩余的体积，但是，

，

• ，为了返回，我只需要重建原始坐标。所以首先我用reshape(v，I，J，K)转换剩余的坐标，我再次创建了一个(I*J*K，4)矩阵，只有当转换后的坐标与新的矩阵匹配时，v列才=1(使用ismember，选项‘SCORE*COEFF'.
• Then’)。
• 我使用reshape(v，i，J，K)创建了索引体积。
• 如果我可视化新的体积回来，它被垂直于图形的主要几何轴切割，就像我需要的那样。

请，我真的很想听到任何意见或建议的方法。