在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文档中此演示的启发。

% 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