4

インデックス付き画像 i、j、k である解剖学的ボリューム画像 (B) があります。

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

ファイルには 0 と 1 のみが含まれます。

等値面で正しく視覚化できます。 isosurface(B);

jのある座標でボリュームをカットしたいと思います（ボリュームごとに異なります）。

問題は、ボリュームが垂直に傾いていることです。おそらく 45% の角度を持っているため、カットは解剖学的ボリュームに沿っていません。

データの新しい直交座標系を取得したいので、座標 j の平面はより正確な方法で解剖学的ボリュームを切断します。

PCAでやるように言われましたが、やり方が分からず、ヘルプページを読んでも役に立ちませんでした。どんな方向でも大歓迎です！

編集：私は推奨事項に従ってきましたが、今では中心がゼロの新しいボリュームを取得しましたが、軸が解剖学的イメージに従っていないと思います。これらは前後の画像です。 元の画像

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 コマンドにあるようですが、以前に行った作業をキャンセルしている可能性はありますか?

4

3 に答える 3

5

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 であることがわかります。

ボリュームの視覚化

于 2013-08-27T20:19:07.543 に答える
3

PCAがあなたの問題を解決するとは思いません。データに PCA を適用すると、3 つの新しい軸が得られます。これらの軸は主成分 (PC) と呼ばれます。それらには、データが投影されたときに最初の PC が最大の分散を持つという特性があります。2 番目の PC は、最初の PC と直交する必要があるという制約の下で、データが投影されるときに最大の分散を持たなければならず、3 番目の PC も同様です。

問題は、データを新しい座標系 (PC によって定義される) に投影したときに、解剖学的ボリュームと一致するかどうかです。必ずしもそうとは限りません。データの適切な軸には、PCA の最適化目的がありません。

于 2013-08-31T14:20:07.517 に答える
0

申し訳ありませんが、@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) がある値よりも大きい場合はすべての値を削除します)。これが私の問題でした。主軸に対して垂直にボリュームをカットしたかったのです。
  • 私にはこれで十分でした。リマイニング座標により、必要なボリュームが得られました。しかし：
  • 残りのボリュームが必要だったので、戻る必要はありませんでしたが、
  • 戻るには、元の座標を再構築する必要がありました。最初に残りの座標を SCORE*COEFF' で変換しました。
  • 次に、(I*J*K, 4) 行列を再度作成し、変換された座標が新しい行列と一致した場合にのみ v 列 = 1 にしました (ismember、オプション 'row' を使用)。
  • reshape(v, I, J, K) を使用してインデックス付きボリュームを作成しました。
  • 新しいボリュームを視覚化すると、必要に応じて、図の主要な幾何学的軸に対して垂直にカットされます。

どうぞ、この方法に関するご意見やご提案をお待ちしております。

于 2013-09-15T21:43:43.050 に答える