三次元データの任意の断面と各種立体画像

CTやMRIの三次元データの任意の断面を表示するプログラムを作成してみます。 当然のようにforループがたくさん出てきます。

Contents

三次元データの任意の断面を表示するhns_section.mの作成

三次元データを中心方向に平行移動し、voxelサイズの大きさを1mmの立方体に変換し、 方位角(Z軸方向)、天頂角(X軸方向)、回転角(Z軸方向)に回転させ、中心点から手前に○mmのところの画像を表示させるプログラムを作成しました。

functionからreturnまでの%を削除してreturnまでをコピーしてhns_section.mという名前で保存してください。 Dはint16の三次元配列、 voxelsizeは直方体の大きさ[mm]、 xdirとydirとzdirは出力画面のX軸・Y軸・Z軸方向の座標[mm]、 anglesは方位角(Z軸)、天頂角(X軸)、回転角(Z軸)の回転角度[degree]です。 回転する軸の順番はZ軸、X軸、Z軸でY軸は使いません。 返り値のRはint16の三次元配列です。

CTやMRIでは三次元配列の型がint16であることが多いので、引数はint16とすることにしました。 断面の座標が三次元配列をはみ出している場合は、0としました。

% function R=hns_section(D,voxelsize,xdir,ydir,zdir,angles)
% % Dはint16の三次元配列
% % voxelsizeは直方体の大きさ[mm]
% % xdirとydirとzdirは出力画面のX軸・Y軸・Z軸方向の座標[mm]
% % anglesは方位角(Z軸)、天頂角(X軸)、回転角(Z軸)の回転角度[degree]
% % 回転する軸の順番はZ軸、X軸、Z軸でY軸は使わない予定
%
% MAX=1024;
% assert(isa(D,'int16')); % int16にしました。
% assert(all(size(D)<=[MAX,MAX,MAX]));
% assert(isa(voxelsize,'double'));
% assert(all(size(voxelsize)==[1,3]));
% assert(isa(xdir,'double'));
% assert(all(size(xdir)<=[1,MAX]));
% assert(isa(ydir,'double'));
% assert(all(size(ydir)<=[1,MAX]));
% assert(isa(zdir,'double'));
% assert(all(size(zdir)<=[1,MAX]));
% assert(isa(angles,'double'));
% assert(all(size(angles)==[1,3]));
%
% center=size(D)/2;
% [xmax,ymax,zmax]=size(D);
% ang1=-angles(1)/180*pi;% Z軸方向の回転(方位角)
% ang2=-angles(2)/180*pi;% X軸方向の回転(天頂角)
% ang3=-angles(3)/180*pi;% Z軸方向の回転(回転角)
% Phi=[cos(ang1),-sin(ang1),0;sin(ang1),cos(ang1),0;0,0,1];
% Theta=[1,0,0;0,cos(ang2),-sin(ang2);0,sin(ang2),cos(ang2)];
% Lambda=[cos(ang3),-sin(ang3),0;sin(ang3),cos(ang3),0;0,0,1];
% M=Phi*Theta*Lambda;
% M=diag(1./voxelsize)*M;
% M=[M,center';0,0,0,1]';
%
% R=repmat(int16(0),[length(xdir),length(ydir),length(zdir)]);
% for m=1:length(zdir)
%     depth=zdir(m);
%     [xx,yy]=meshgrid(ydir,xdir);% meshgridはX軸とY軸が逆になる
%     P1=[xx(:),yy(:),ones(numel(xx),1)*[depth,1]];
%     P0=P1*M;
%     for n=1:size(P0,1)
%         x=P0(n,1);x0=floor(x);
%         if x0<1;P0(n,4)=0;continue;end;
%         x1=x0+1;
%         if x1>xmax;P0(n,4)=0;continue;end;
%
%         y=P0(n,2);y0=floor(y);y1=y0+1;dy=y-y0;
%         if y0<1;P0(n,4)=0;continue;end;
%         y1=y0+1;
%         if y1>ymax;P0(n,4)=0;continue;end;
%
%         z=P0(n,3);z0=floor(z);z1=z0+1;dz=z-z0;
%         if z0<1;P0(n,4)=0;continue;end;
%         z1=z0+1;
%         if z1>zmax;P0(n,4)=0;continue;end;
%
%         dx=x-x0;xd=1-dx;
%         dy=y-y0;yd=1-dy;
%         dz=z-z0;zd=1-dz;
%         v1=xd*yd*zd*double(D(x0,y0,z0));
%         v2=dx*yd*zd*double(D(x1,y0,z0));
%         v3=xd*dy*zd*double(D(x0,y1,z0));
%         v4=dx*dy*zd*double(D(x1,y1,z0));
%         v5=xd*yd*dz*double(D(x0,y0,z1));
%         v6=dx*yd*dz*double(D(x1,y0,z1));
%         v7=xd*dy*dz*double(D(x0,y1,z1));
%         v8=dx*dy*dz*double(D(x1,y1,z1));
%         P0(n,4)=v1+v2+v3+v4+v5+v6+v7+v8;
%     end;
%     v=reshape(P0(:,4),[length(xdir),length(ydir)]);
%     R(:,:,m)=int16(v);
% end;
% return;

hns_section.のMEX化とMRIデータの読み込み

予めuint8のDをint16に変換します。 MEX化することで処理時間が短くなりました。

codegen('hns_section.m');
load mri;% MATLABのサンプルデータ
D=int16(squeeze(D));

方位角45度、天頂角90度、回転角0度に回転した時の画像

hns\_section.mは本来X-Y平面用の前後方向で端っこのスライスはデータがないので真っ暗です。 MEX化したら処理時間は半分程度になりました。

voxelsize=[2,2,5];
xdir=linspace(-128,128,512);
ydir=linspace(-128,128,512);
zdir=linspace(-128,128,10);
angles=[45,90,0];
tic;
D1=hns_section(D,voxelsize,xdir,ydir,zdir,angles);
st1=toc;
tic;
D2=hns_section_mex(D,voxelsize,xdir,ydir,zdir,angles);
st2=toc;

M=length(zdir);
figure('color',[1,1,1]);
for n=1:2
    switch n
        case 1;D0=D1;st=sprintf('hns\\_section.mの処理時間は%0.4f秒です',st1);
        case 2;D0=D2;st=sprintf('hns\\_section\\_mexの処理時間は%0.4f秒です',st2);
    end;
    for m=1:M
        subplot(2*2,M/2,m+(n-1)*M);
        imagesc(xdir,ydir,squeeze(D0(:,:,m)));
        if m==2;title(st);end;
        axis tight;
        daspect([1,1,1]);
    end;
end;
colormap(gray);
set(findobj(gcf,'type','axes'),'clim',[0,100]);

MATLABのcontourslice関数を使った表現

なんか使えそうな気もするんですが・・・。

voxelsize=[2,2,5];
xdir=linspace(-128,128,512);
ydir=linspace(-128,128,256);
zdir=-110:20:110;
angles=[45,90,180];
D0=hns_section_mex(D,voxelsize,xdir,ydir,zdir,angles);
figure('color',[1,1,1]);
xx=1:size(D0,1);yy=1:size(D0,2);zz=1:size(D0,3);
[x,y,z]=meshgrid(yy,xx,zz);% yyとxxは逆にすること!
contourslice(x,y,z,D0,[],[],zz);
axis tight;daspect([0.5,1,0.02]);
colormap(gray);
view([40,20]);

最大値投射法(Maximum Intensity Projection)

MR血管撮影で使われる方法です。 生憎 Time of Flight法でなく、通常のSpin Echo法のようで血管の信号が小さいため 血管は描出できていませんが、原理は最大値!という単純なものです。

voxelsize=[2,2,5];
xdir=linspace(-128,128,256);
ydir=linspace(-128,128,256);
zdir=linspace(-128,128,256);
angles=[-45,70,0];
tic;
D0=hns_section_mex(D,voxelsize,xdir,ydir,zdir,angles);
toc
d1=max(D0,[],3);
figure('color',[1,1,1]);
imagesc(xdir,ydir,d1);
colormap(gray);
daspect([1,1,1]);
xlabel(sprintf('[%d %d %d]',angles));
colorbar;
title('最大値投射法');
経過時間は 1.162417 秒です。

平均値を用いた透視風画像の作成

平均値を使うと透視風の画像を作成することができます。CTだとたぶんうまくいく筈ですが、MRIなんで雰囲気だけです。

d2=mean(D0,3);
figure('color',[1,1,1]);
imagesc(xdir,ydir,d2);
colormap(gray);
daspect([1,1,1]);
xlabel(sprintf('[%d %d %d]',angles));
colorbar;
title('平均値を用いた透視風画像');

閾値を利用した透視風画像の作成

閾値以上のvoxelだけで透視風画像を作成しました。 forループが出てくるので処理時間は遅くなってきますが、画像はくっきりしてきます。 閾値の設定は画面中心[128,128]の前後方向のデータから、とりあえず10としました。

figure('color',[1,1,1]);
plot(squeeze(D0(128,128,:)));hold on;
threshold=10;
plot([1,length(zdir)],threshold*[1,1],'color',[1,0,0]);

D1=D0;
D1(D1<threshold)=0;
nx=length(xdir);
ny=length(ydir);
nz=length(zdir);
d3=zeros(nx,ny);

for y=1:ny
    for x=1:nx
        k=squeeze(double(D1(x,y,:)));
        d3(x,y)=max(k);
        n=(k>=threshold);
        nn=sum(n);
        if nn>2;% ごみ取り
            d3(x,y)=sum(k)/nn;
        end;
    end;
end;
figure('color',[1,1,1]);
imagesc(xdir,ydir,d3);
colormap(gray);
daspect([1,1,1]);
xlabel(sprintf('[%d %d %d]',angles));
colorbar;
title('閾値>10の透視風画像');

画像の前後方向の情報から画像を作成

閾値を設定して画像の前後方向の情報(Z buffer)だけを使って作図します。 最も近いところの画像と最も遠いところの画像の2つができます。

d4=zeros(nx,ny);
d5=d4;
for y=1:ny
    for x=1:nx
        k=squeeze(double(D1(x,y,:)));
        n=(k>threshold);
        nn=find(n==1);% ゴミ取り
        if ~isempty(nn);
            d4(x,y)=nn(end);% 最も手前
            d5(x,y)=nn(1);% 最も遠いところ
        end;
    end;
end;

 for n=1:2
    figure('color',[1,1,1]);
     switch n
         case 1;d=d4;st='前';
         case 2;d=d5;st='後';
     end;
     imagesc(xdir,ydir,d);
     daspect([1,1,1]);
     colorbar;
     xlabel(sprintf(' [%d %d %d]',angles));
     colormap(gray);
     title(['Z buffer画像 ',st]);
 end;

簡単なvolume rendering画像

手前方向のZbufferを使って、その点からN画素分深部のデータを平均することで volume rendering画像を作成しました。Nが20画素位になると、うっすら脳表がでてきます。

d6=zeros(nx,ny); % 5画素
d7=d6;           % 10画素
d8=d6;           % 20画素
for y=1:ny
    for x=1:nx
        k=squeeze(double(D1(x,y,:)));
        d2(x,y)=mean(k);
        n=(k>threshold);
        nn=find(n==1);
        if ~isempty(nn);
            nnn=nn(end)+(-4:0);% 5画素
            nnn(nnn<1)=[];
            d6(x,y)=mean(k(nnn));
            nnn=nn(end)+(-9:0);% 10画素
            nnn(nnn<1)=[];
            d7(x,y)=mean(k(nnn));
            nnn=nn(end)+(-19:0);% 20画素
            nnn(nnn<1)=[];
            d8(x,y)=mean(k(nnn));
        end;
    end;
end;

 for n=1:3
    figure('color',[1,1,1]);
     switch n
         case 1;d=d6;st='5画素';
         case 2;d=d7;st='10画素';
         case 3;d=d8;st='20画素';
     end;
     imagesc(xdir,ydir,d);
     daspect([1,1,1]);
     colorbar;
     xlabel(sprintf(' [%d %d %d]',angles));
     colormap(gray);
     title(['簡単なvolume rendering ',st]);
 end;

頭皮と脳の隙間の画素の閾値を利用して脳表画像作成の試み

頭皮と脳表の間の隙間の閾値が30くらいなので、この値を利用して脳表画像の作成を試みてみます。 頭皮除去は甘く、顔面頭蓋の軟部組織は全然除去できていません。 実際の脳表画像は脳実質のvoxelをマスクして、volume renderingして作成します。

figure('color',[1,1,1]);
k=squeeze(D0(128,128,:));
plot(k);hold on;
threshold=10;scalp=40;
plot([1,nz],threshold*[1,1],'color',[1,0,0]);
plot([1,nz],scalp*[1,1],'color',[0,1,0]);
title('[x,y]=[128,128]の前後方向の値');

figure('color',[1,1,1]);
kk=k;
kk(kk<scalp)=0;
kk(kk>0)=1;
thickness=10;
for m=(length(k)-1):-1:1
    if kk(m)>0
        kk(m)=kk(m+1)+1;
    end;
end;
plot(kk);hold on;
plot([1,nz],thickness*[1,1],'color',[1,0,0]);
title('脳表抽出のため頭皮>50 厚さscalp>10以上で検出');

d9=zeros(nx,ny);
for y=1:ny
    for x=1:nx
        k=squeeze(double(D1(x,y,:)));
        kk=k;
        kk(kk<scalp)=0;
        kk(kk>0)=1;
        for m=(length(k)-1):-1:1
            if kk(m)>0
                kk(m)=kk(m+1)+1;
            end;
        end;
        mm=find(kk==thickness);
        if ~isempty(mm);
            mmm=mm(end)+(-9:0)+thickness;
            mmm(mmm>nz)=[];
            mmm(mmm<1)=[];
            kk=mean(k(mmm));
            d9(x,y)=kk;
        end;
    end;
end;

figure('color',[1,1,1]);
imagesc(xdir,ydir,d9);
colormap(gray);
daspect([1,1,1]);
xlabel(sprintf('[%d %d %d]',angles));
colorbar;
title('頭皮>40 厚さ>10 で脳表画像作成');