3D 图像上的 PCA 以获得 3 个主轴

Question

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

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

该文件仅包含 0 和 1。

我可以用等值面正确地可视化它： isosurface(B);

我想在 j 中的某个坐标处切割体积（每个体积都不同）。

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

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

有人告诉我用 PCA 来做，但我不知道如何做，阅读帮助页面也没有帮助。欢迎任何方向！

编辑： 我一直在遵循建议，现在我得到了一个以零为中心的新卷，但我认为轴没有按照应有的解剖图像。这些是前后图像：

这是我用来创建图像的代码（我从答案和 cmets 的想法中获取了一些代码）：

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,J,K) 3D 体积。
• 根据@werner 的建议，我将其更改为 4 列矩阵 (x,y,z,v)，大小 (I*J*K, 4)。
• 当 v == 0 时消除坐标 (x,y,z)，并且 v 也是如此。所以现在，大小是（原始卷，3）。只有值为1的坐标，所以向量的长度就是图形的体积。
• 执行 PCA 以获得 COEFF 和 SCORE。
• COEFF 是一个 3x3 矩阵。它是单一且正交的，它是我的体数据的旋转矩阵。
• 我在 PCA1 轴上进行了编辑（即当 COEFF(1) 大于某个值时删除 al 值）。这是我的问题，我想垂直于主轴切割体积。
• 这对我来说已经足够了，扩孔坐标给了我想要的音量。但是：
• 我不需要返回，因为我只需要剩余的音量，但是
• 为了返回，我只需要重建原始坐标。所以首先我用 SCORE*COEFF' 转换了剩余的坐标。
• 然后我再次创建了一个 (I*J*K, 4) 矩阵，仅当转换后的坐标与新矩阵匹配时才使 v 列 = 1（使用 ismember，选项“行”）。
• 我使用 reshape(v, I, J, K) 创建了索引卷。
• 如果我将新卷向后可视化，它会按照我的需要垂直于图形的主要几何轴进行切割。

拜托，我真的很想听听有关该方法的任何评论或建议。