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

Question

我有一个解剖体积图像（B），它是一个索引图像I，J，K：

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

的文件只包含0和1。

我可以与等值面正确地进行可视化： 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 

任何人都可以提供什么可能发生的提示吗？看起来问题可能是重塑命令，是否有可能取消以前完成的工作？

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 

下面是表示切片平面渲染为线框等值面，除了结果两个轴对齐和一个旋转。我们可以看到，在等值面的内部体积空间在外面

Answer 2

我不认为PCA解决您的问题等于1，和值0。如果将PCA应用于您的数据，它将为您提供三个新轴。这些轴被称为主要组件（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）矩阵，只有当变换后的坐标匹配新矩阵（带有ismember，选项'row'）时，v列= 1。
• 我使用reshape（v，I，J，K）创建了索引音量。
• 如果我想象新的体积，它会垂直于图形的主要几何轴线切割，就像我需要的一样。

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