matlab - PCA 在 3D 图像上获得 3 个主轴

Question

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

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

该文件仅包含 0 和 1。

我可以用等值面正确地可视化它:等值面(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 标量体积数据，并在倾斜平面(非轴对齐)切割体积。代码的灵感来自 this demo在 MATLAB 文档中。

% 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