簡易版脳ボクセル抽出

脳ボクセルの抽出は閾値と連続性を利用します。forループの塊です。 MEX化しないとかなり時間がかかります。

Contents

閾値別に表示する

DICOM値が40~60のボクセルを赤で表示します。 全スライスを表示すると煩雑なので6枚だけにします。

load mri;
D=squeeze(D);
D=int16(D);
[dx,dy,dz]=size(D);
threshold=[40,80];
zz=round(linspace(1,dz,8));% 適当なスライス表示用
zz([1,end])=[];
D0=zeros(dx,dy,3);
map=gray(256);
map2=map;
map2(:,2:3)=0;
clim=[min(D(:)),max(D(:))];
w=256/(double(clim(2)-clim(1)));
figure('color',[1,1,1]);
for n=1:length(zz)
    for x=1:dx
        for y=1:dy
            k=D(x,y,zz(n));
            kk=round((double(k)-clim(1))*w);
            if kk<1;kk=1;
            elseif kk>256;kk=256;
            end;
            if k>=threshold(1) && k<=threshold(2)
                D0(x,y,:)=map2(kk,:);
            else
                D0(x,y,:)=map(kk,:);
            end;
        end;
    end;
    subplot(2,3,n);
    image(D0);daspect([1,1,1]);
    title(sprintf('スライス%d枚目',zz(n)));
end;

閾値内のボクセルをマスク領域(青)に設定

閾値内の赤のボクセルを青にします。

% functionからreturnまでの%を削除してreturnまでコピーしてhns_MaskVoxel.mという名前で保存してください。
% 引数のDはint16の三次元配列、
% Dmaskは同じサイズのuint8の三次元配列、
% mは目標となるマスク用のビット数、
% thresholdはdouble 1x2 閾値、
% 返り値はuint8の三次元配列です。

% function Dmask=hns_MaskVoxel(D,Dmask,m,threshold)
% MAX=1024;
% assert(isa(D,'int16'));
% assert(all(size(D)<=[MAX,MAX,MAX]));
% assert(isa(Dmask,'uint8'));
% assert(all(size(Dmask)<=[MAX,MAX,MAX]));
% assert(isa(m,'double'));
% assert(isa(threshold,'double'));
% assert(all(size(threshold)==[1,2]));
% if ~all(size(D)==size(Dmask));return;end;
% [dx,dy,dz]=size(D);
% for x=1:dx
%     for y=1:dy
%         for z=1:dz
%             if (D(x,y,z)>=threshold(1)) && (D(x,y,z)<=threshold(2))
%                 Dmask(x,y,z)=uint8(bitset(Dmask(x,y,z),m,1));
%             else
%                 Dmask(x,y,z)=uint8(bitset(Dmask(x,y,z),m,0));
%             end;
%         end;
%     end;
% end;
%

Dmask=repmat(uint8(0),[dx,dy,dz]);
m=1;
tic;
Dmask0=hns_MaskVoxel(D,Dmask,m,threshold);
toc
codegen('hns_MaskVoxel.m');
tic;
Dmask0=hns_MaskVoxel_mex(D,Dmask,m,threshold);
toc
Dmask=Dmask0;
clear Dmask0;
経過時間は 0.874746 秒です。
経過時間は 0.133915 秒です。

青マスクを1ボクセル縮小する

forループと場合分けが出てきます。 以下の関数Mファイルを作成しました。

functionからreturnまでの%を削除してreturnまでコピーしてhns_ShrinkVoxel.mという名前で保存してください。 引数のDはuint8の三次元配列、 mは目標となるマスク用のビット数、 Nは1ボクセル縮小する処置を何回繰り返すかです。

% function R=hns_ShrinkVoxel(D,m,N)
% MAX=1024;
% assert(isa(D,'uint8'));
% assert(all(size(D)<=[MAX,MAX,MAX]));
% assert(isa(m,'double'));
% assert(isa(N,'double'));
% [dx,dy,dz]=size(D);
% R=repmat(uint8(0),[dx,dy,dz]);
% D3=max(D,[],3);
% D3(D3<2^(m-1))=0;
% D23=max(D3,[],2);
% D23(D23<2^(m-1))=0;
% x=find(D23>0);
% xmin=min(x);
% xmax=max(x);
% D13=max(D3,[],1);
% D13(D13<2^(m-1))=0;
% y=find(D13>0);
% ymin=min(y);
% ymax=max(y);
% D12=max(max(D,[],1),[],2);
% z=find(D12>0);
% zmin=min(z);
% zmax=max(z);
% for n=1:N
%     for x=xmin:xmax
%         for y=ymin:ymax
%             for z=zmin:zmax
%                 if x>1
%                     if y>1
%                         if z>1;if bitget(D(x-1,y-1,z-1),m)==0;R(x,y,z)=bitset(R(x,y,z),m,0);continue;end;end;
%                         if bitget(D(x-1,y-1,z),m)  ==0;R(x,y,z)=bitset(R(x,y,z),m,0);continue;end;
%                         if z<dz;if bitget(D(x-1,y-1,z+1),m)==0;R(x,y,z)=bitset(R(x,y,z),m,0);continue;end;end;
%                     end;
%                     if z>1;if bitget(D(x-1,y,  z-1),m)==0;R(x,y,z)=bitset(R(x,y,z),m,0);continue;end;end;
%                     if bitget(D(x-1,y,  z),m)  ==0;R(x,y,z)=bitset(R(x,y,z),m,0);continue;end;
%                     if z<dz;if bitget(D(x-1,y,  z+1),m)==0;R(x,y,z)=bitset(R(x,y,z),m,0);continue;end;end;
%                     if y<dy
%                         if z>1;if bitget(D(x-1,y+1,z-1),m)==0;R(x,y,z)=bitset(R(x,y,z),m,0);continue;end;end;
%                         if bitget(D(x-1,y+1,z),m)  ==0;R(x,y,z)=bitset(R(x,y,z),m,0);continue;end;
%                         if z<dz;if bitget(D(x-1,y+1,z+1),m)==0;R(x,y,z)=bitset(R(x,y,z),m,0);continue;end;end;
%                     end;
%                 end;
%                 if y>1
%                     if z>1;if bitget(D(x,  y-1,z-1),m)==0;R(x,y,z)=bitset(R(x,y,z),m,0);continue;end;end;
%                     if bitget(D(x,  y-1,z),m)  ==0;R(x,y,z)=bitset(R(x,y,z),m,0);continue;end;
%                     if z<dz;if bitget(D(x,  y-1,z+1),m)==0;R(x,y,z)=bitset(R(x,y,z),m,0);continue;end;end;
%                 end;
%                 if z>1;if bitget(D(x,  y,  z-1),m)==0;R(x,y,z)=bitset(R(x,y,z),m,0);continue;end;end;
%                 if bitget(D(x,  y,  z),m)  ==0;R(x,y,z)=bitset(R(x,y,z),m,0);continue;end;
%                 if z<dz;if bitget(D(x,  y,  z+1),m)==0;R(x,y,z)=bitset(R(x,y,z),m,0);continue;end;end;
%                 if y<dy
%                     if z>1;if bitget(D(x,  y+1,z-1),m)==0;R(x,y,z)=bitset(R(x,y,z),m,0);continue;end;end;
%                     if bitget(D(x,  y+1,z),m)  ==0;R(x,y,z)=bitset(R(x,y,z),m,0);continue;end;
%                     if z<dz;if bitget(D(x,  y+1,z+1),m)==0;R(x,y,z)=bitset(R(x,y,z),m,0);continue;end;end;
%                 end;
%                 if x<dx
%                     if y>1
%                         if z>1;if bitget(D(x+1,y-1,z-1),m)==0;R(x,y,z)=bitset(R(x,y,z),m,0);continue;end;end;
%                         if bitget(D(x+1,y-1,z),m)  ==0;R(x,y,z)=bitset(R(x,y,z),m,0);continue;end;
%                         if z<dz;if bitget(D(x+1,y-1,z+1),m)==0;R(x,y,z)=bitset(R(x,y,z),m,0);continue;end;end;
%                     end;
%                     if z>1;if bitget(D(x+1,y,  z-1),m)==0;R(x,y,z)=bitset(R(x,y,z),m,0);continue;end;end;
%                     if bitget(D(x+1,y,  z),m)  ==0;R(x,y,z)=bitset(R(x,y,z),m,0);continue;end;
%                     if z<dz;if bitget(D(x+1,y,  z+1),m)==0;R(x,y,z)=bitset(R(x,y,z),m,0);continue;end;end;
%                     if y<dy
%                         if z>1;if bitget(D(x+1,y+1,z-1),m)==0;R(x,y,z)=bitset(R(x,y,z),m,0);continue;end;end;
%                         if bitget(D(x+1,y+1,z),m)  ==0;R(x,y,z)=bitset(R(x,y,z),m,0);continue;end;
%                         if z<dz;if bitget(D(x+1,y+1,z+1),m)==0;R(x,y,z)=bitset(R(x,y,z),m,0);continue;end;end;
%                     end;
%                 end;
%                 R(x,y,z)=bitset(R(x,y,z),m,1);
%             end;
%         end;
%     end;
%     k=find(bitget(R,m)==0);
%     D(k)=bitset(D(k),m,0);
%     if numel(k)==dx*dy*dz;break;end;
% end;
% R=D;
% return;

tic;
Dmask0=hns_ShrinkVoxel(Dmask,1,1);
toc
codegen('hns_ShrinkVoxel.m');
tic;
Dmask0=hns_ShrinkVoxel_mex(Dmask,1,1);%
toc
Dmask=Dmask0;
clear Dmask0;% メモリ節約のため
map3=map;
map3(:,1:2)=0;
figure('color',[1,1,1]);
for n=1:length(zz)
    for x=1:dx
        for y=1:dy
            k=D(x,y,zz(n));
            kk=round((double(k)-clim(1))*w);
            if kk<1;kk=1;
            elseif kk>256;kk=256;
            end;
            if bitget(Dmask(x,y,zz(n)),1)==1
                D0(x,y,:)=map3(kk,:);
            else
                if k>=threshold(1) && k<=threshold(2)
                    D0(x,y,:)=map2(kk,:);
                else
                    D0(x,y,:)=map(kk,:);
                end;
            end;
        end;
    end;
    subplot(2,3,n);
    image(D0);daspect([1,1,1]);hold on;
    title(sprintf('スライス%d枚目',zz(n)));
end;
経過時間は 3.099521 秒です。
経過時間は 0.314167 秒です。

領域選択

青色の中でimfreehand関数で領域を選択します。マウスボタンを離し、領域をダブルクリックすれば、 選択した領域が水色になります。 とりあえず5枚目のスライスで始めます。

haxes=findobj(gcf,'type','axes');
haxes=sort(haxes);
figure('color',[1,1,1]);
copyobj(haxes(5),gcf);hold on;
set(gca,'outerposition',[0,0,1,1]);
h=imfreehand;
position=wait(h);

hCloseLineTop=findobj(h,'type','line','tag','close line top');% imfreehandでマウスを放した後、閉じられた線
hTopLine=findobj(h,'type','line','tag','top line');% imfreehandで引いた線
P=[get(hTopLine,'xdata'),get(hCloseLineTop,'xdata');
   get(hTopLine,'ydata'),get(hCloseLineTop,'ydata')];
k=fill(P(1,:),P(2,:),[0,1,1]);

選択した領域をマスク

領域内外判定には色々な方法が有ると思いますが、ここではベクトルの外積を使います。 隣接する領域の点と任意の点を頂点とする三角形の任意の点側の角度の総和が±2π、±6π、±10π・・・になることを利用します。 個々の三角形の角度は±90度以下とします。 そこで以下のような関数Mファイルを作成しました。

functionからreturnまでの%を削除してreturnまでコピーしてhns_RegionCheck.mという名前で保存してください。 引数のPは領域の協会のNX座標とY座標個のからなるdoubleの2×Nの行列です。 dxとdyは元の画像のX軸方向とY軸方向のボクセル数です。 返り値で奇数が領域内、偶数が領域外となります。 尚境界線上は領域内としています。

imfreehand関数で設定した領域が緑で表示されています。

% function R=hns_RegionCheck(P,dx,dy)
% MAX=1024;
% assert(isa(P,'double'));
% assert(size(P,1)==2);
% assert(size(P,2)<=MAX);
% assert(isa(dx,'double'));
% assert(isa(dy,'double'));
%
% R=zeros(dx,dy);
% nP=size(P,2);
% NP=ones(1,nP);
% ang=zeros(1,nP);
% xmin=floor(min(P(1,:)));
% if xmin<1;xmin=1;end;
% xmax=ceil(max(P(1,:)));
% if xmax>dx;xmax=dx;end;
% ymin=floor(min(P(2,:)));
% if ymin<1;ymin=1;end;
% ymax=ceil(max(P(2,:)));
% if ymax>dy;ymax=dy;end;
% Norm=zeros(1,nP);
% for x=xmin:xmax
%     for y=ymin:ymax
%         p=P-[x;y]*NP;
%         for n=1:nP
%             Norm(n)=sqrt(sum(p(:,n).^2));
%         end;
%         for n=1:(nP-1)
%             NN=Norm(n)*Norm(n+1);
%             if NN<1e-14;R(x,y)=2*pi;continue;end;
%             ang(n)=asin((p(1,n)*p(2,n+1)-p(1,n+1)*p(2,n))/NN);
%         end;
%         R(x,y)=sum(ang);
%     end;
% end;
% R=R'/(2*pi);% 何故か逆
% return;

hf=gcf;
tic;
K=hns_RegionCheck_mex(P,dx,dy);
toc
codegen('hns_RegionCheck.m');
tic;
K=hns_RegionCheck_mex(P,dx,dy);
toc
figure('color',[1,1,1]);
subplot(121);imagesc(K);colorbar;title('領域の内外判定');
K=mod(abs(round(K)),2);
subplot(122);imagesc(K);colorbar;title('領域の内外判定 二分法');
m=2;
for x=1:dx
    for y=1:dy
        if K(x,y)==1
            Dmask(x,y,zz(5))=uint8(bitset(Dmask(x,y,zz(5)),m,1));
        else
            Dmask(x,y,zz(5))=uint8(bitset(Dmask(x,y,zz(5)),m,0));
        end;
    end;
end;

map4=map;
map4(:,[1,3])=0;
figure('color',[1,1,1]);
for n=1:length(zz)
    for x=1:dx
        for y=1:dy
            k=D(x,y,zz(n));
            kk=round((double(k)-clim(1))*w);
            if kk<1;kk=1;
            elseif kk>256;kk=256;
            end;
            if bitget(Dmask(x,y,zz(n)),2)==1
                D0(x,y,:)=map4(kk,:);
            elseif bitget(Dmask(x,y,zz(n)),1)==1
                D0(x,y,:)=map3(kk,:);
            else
                if k>=threshold(1) && k<=threshold(2)
                    D0(x,y,:)=map2(kk,:);
                else
                    D0(x,y,:)=map(kk,:);
                end;
            end;
        end;
    end;
    subplot(2,3,n);
    image(D0);daspect([1,1,1]);hold on;
    title(sprintf('スライス%d枚目',zz(n)));
end;
経過時間は 0.228467 秒です。
経過時間は 0.328338 秒です。

青マスク内で緑マスクの拡張

青内で緑のマスクを拡張します。1ボクセルずつ拡張して、 拡張するボクセルがなくなったら終了です。 そこで以下のような関数Mファイルを作成しました。

functionからreturnまでの%を削除してreturnまでコピーしてhns_ExpandVoxel.mという名前で保存してください。 引数のDはuint8の三次元配列です。8ビットを想定し、n bit内でm bitをN回拡張することを想定しています。 n,m,Nはdoubleです。拡張回数を1000位にすると、途中で拡張の余地がないと判断され収束します。

% function R=hns_ExpandVoxel(D,n,m,N)
% MAX=1024;
% assert(isa(D,'uint8'));
% assert(all(size(D)<=[MAX,MAX,MAX]));
% assert(isa(n,'double'));% n bit内で
% assert(isa(m,'double'));% m bitを拡張
% assert(isa(N,'double'));% N回拡張を繰り返す
% [dx,dy,dz]=size(D);
% R=D;
% n=uint8(n);
% m=uint8(m);
% for nn=1:N
%     for x=1:dx
%         for y=1:dy
%             for z=1:dz
%                 k=D(x,y,z);
%                 if bitget(k,m)==1;continue;end;
%                 if bitget(k,n)==0;continue;end;
%                 if x>1
%                     if y>1
%                         if z>1;if bitget(D(x-1,y-1,z-1),m); R(x,y,z)=uint8(bitset(k,m,1));continue;end;end;
%                         if bitget(D(x-1,y-1,z),m);          R(x,y,z)=uint8(bitset(k,m,1));continue;end;
%                         if z<dz;if bitget(D(x-1,y-1,z+1),m);R(x,y,z)=uint8(bitset(k,m,1));continue;end;end;
%                     end;
%                     if z>1;if bitget(D(x-1,y,z-1),m); R(x,y,z)=uint8(bitset(k,m,1));continue;end;end;
%                     if bitget(D(x-1,y,z),m);          R(x,y,z)=uint8(bitset(k,m,1));continue;end;
%                     if z<dz;if bitget(D(x-1,y,z+1),m);R(x,y,z)=uint8(bitset(k,m,1));continue;end;end;
%                     if y<dy
%                         if z>1;if bitget(D(x-1,y+1,z-1),m); R(x,y,z)=uint8(bitset(k,m,1));continue;end;end;
%                         if bitget(D(x-1,y+1,z),m);          R(x,y,z)=uint8(bitset(k,m,1));continue;end;
%                         if z<dz;if bitget(D(x-1,y+1,z+1),m);R(x,y,z)=uint8(bitset(k,m,1));continue;end;end;
%                     end;
%                 end;
%                 if y>1
%                     if z>1;if bitget(D(x,y-1,z-1),m); R(x,y,z)=uint8(bitset(k,m,1));continue;end;end;
%                     if bitget(D(x,y-1,z),m);          R(x,y,z)=uint8(bitset(k,m,1));continue;end;
%                     if z<dz;if bitget(D(x,y-1,z+1),m);R(x,y,z)=uint8(bitset(k,m,1));continue;end;end;
%                 end;
%                 if z>1;if bitget(D(x,y,z-1),m); R(x,y,z)=uint8(bitset(k,m,1));continue;end;end;
%                 if bitget(D(x,y,z),m);          R(x,y,z)=uint8(bitset(k,m,1));continue;end;
%                 if z<dz;if bitget(D(x,y,z+1),m);R(x,y,z)=uint8(bitset(k,m,1));continue;end;end;
%                 if y<dy
%                     if z>1;if bitget(D(x,y+1,z-1),m); R(x,y,z)=uint8(bitset(k,m,1));continue;end;end;
%                     if bitget(D(x,y+1,z),m);          R(x,y,z)=uint8(bitset(k,m,1));continue;end;
%                     if z<dz;if bitget(D(x,y+1,z+1),m);R(x,y,z)=uint8(bitset(k,m,1));continue;end;end;
%                 end;
%                 if x<dx
%                     if y>1
%                         if z>1;if bitget(D(x+1,y-1,z-1),m); R(x,y,z)=uint8(bitset(k,m,1));continue;end;end;
%                         if bitget(D(x+1,y-1,z),m);          R(x,y,z)=uint8(bitset(k,m,1));continue;end;
%                         if z<dz;if bitget(D(x+1,y-1,z+1),m);R(x,y,z)=uint8(bitset(k,m,1));continue;end;end;
%                     end;
%                     if z>1;if bitget(D(x+1,y,z-1),m); R(x,y,z)=uint8(bitset(k,m,1));continue;end;end;
%                     if bitget(D(x+1,y,z),m);          R(x,y,z)=uint8(bitset(k,m,1));continue;end;
%                     if z<dz;if bitget(D(x+1,y,z+1),m);R(x,y,z)=uint8(bitset(k,m,1));continue;end;end;
%                     if y<dy
%                         if z>1;if bitget(D(x+1,y+1,z-1),m); R(x,y,z)=uint8(bitset(k,m,1));continue;end;end;
%                         if bitget(D(x+1,y+1,z),m);          R(x,y,z)=uint8(bitset(k,m,1));continue;end;
%                         if z<dz;if bitget(D(x+1,y+1,z+1),m);R(x,y,z)=uint8(bitset(k,m,1));continue;end;end;
%                     end;
%                 end;
%             end;
%         end;
%     end;
%     check=0;
%     for x=1:dx
%         for y=1:dy
%             for z=1:dz
%                 if bitget(D(x,y,z),m)==0
%                     if bitget(R(x,y,z),m)==1
%                         D(x,y,z)=R(x,y,z);
%                         check=check+1;
%                     end;
%                 end;
%             end;
%         end;
%     end;
%     if check==0;break;end;
% end;
%
% return;

tic;
Dmask1=hns_ExpandVoxel(Dmask,1,2,1000);
toc
codegen('hns_ExpandVoxel.m');
tic;
Dmask1=hns_ExpandVoxel_mex(Dmask,1,2,1000);
toc
figure('color',[1,1,1]);
for n=1:length(zz)
    for x=1:dx
        for y=1:dy
            k=D(x,y,zz(n));
            kk=round((double(k)-clim(1))*w);
            if kk<1;kk=1;
            elseif kk>256;kk=256;
            end;
            if bitget(Dmask1(x,y,zz(n)),2)==1
                D0(x,y,:)=map4(kk,:);
            elseif bitget(Dmask1(x,y,zz(n)),1)==1
                D0(x,y,:)=map3(kk,:);
            else
                if k>=threshold(1) && k<=threshold(2)
                    D0(x,y,:)=map2(kk,:);
                else
                    D0(x,y,:)=map(kk,:);
                end;
            end;
        end;
    end;
    subplot(2,3,n);
    image(D0);daspect([1,1,1]);hold on;
    title(sprintf('スライス%d枚目',zz(n)));
end;
経過時間は 51.455075 秒です。
経過時間は 0.817151 秒です。

閾値内で緑マスクの拡張

再度閾値を青でマスクし、青マスク内で緑マスクを1ボクセル拡張します。 これで緑マスクが脳ボクセルです。とりあえず、これで脳ボクセルの切り出しが終了します。 実際には緑マスクが脳ボクセル以外を含んでいたり、脳ボクセルが含まれていなかったりします。 閾値の設定や、青マスクのボクセル縮小・拡張をの回数を調整して、試行錯誤しながら脳の切り出しを行うことになります。

Dmask1=hns_MaskVoxel(D,Dmask1,1,threshold);
Dmask1=hns_ExpandVoxel_mex(Dmask1,1,2,1);

figure('color',[1,1,1]);
for n=1:length(zz)
    for x=1:dx
        for y=1:dy
            k=D(x,y,zz(n));
            kk=round((double(k)-clim(1))*w);
            if kk<1;kk=1;
            elseif kk>256;kk=256;
            end;
            if bitget(Dmask1(x,y,zz(n)),2)==1
                D0(x,y,:)=map4(kk,:);
            elseif bitget(Dmask1(x,y,zz(n)),1)==1
                D0(x,y,:)=map3(kk,:);
            else
                if k>=threshold(1) && k<=threshold(2)
                    D0(x,y,:)=map2(kk,:);
                else
                    D0(x,y,:)=map(kk,:);
                end;
            end;
        end;
    end;
    subplot(2,3,n);
    image(D0);daspect([1,1,1]);hold on;
    title(sprintf('スライス%d枚目',zz(n)));
end;

脳表画像の作成

緑マスク以外を0にして脳ボクセルを抽出し、方位角-135度、天頂角70度、回転角0度で回転させ、 Volume renderingを作成してみます。 背景を白にするため背景を-1にしてcolrobarの最上位を白にし、背景を最上位の値に変更しています。

[dx,dy,dz]=size(D);
DD=D;
tic;
for x=1:dx
    for y=1:dy
        for z=1:dz
            if bitget(Dmask1(x,y,z),2)==0
                DD(x,y,z)=0;
            end;
        end;
    end;
end;
toc
voxelsize=[2,2,5];
xdir=linspace(-128,128,256);
ydir=linspace(-128,128,256);
zdir=linspace(-128,128,256);
angles=[-135,70,0];
D1=hns_section_mex(D, voxelsize,xdir,ydir,zdir,angles);
D2=hns_section_mex(DD,voxelsize,xdir,ydir,zdir,angles);

ny=length(ydir);
nx=length(xdir);

d1=zeros(nx,ny);
d2=zeros(nx,ny);
map1=gray;
map2=map1;map2(:,3)=0;
map2(end,:)=1;

for m=1:2
    for y=1:ny
        for x=1:nx
            if m==1;
                k=squeeze(double(D1(x,y,:)));
                n=(k>threshold(1));
            else
                k=squeeze(double(D2(x,y,:)));
                n=(k>0);
            end;
            nn=find(n==1);
            if ~isempty(nn);
                nnn=nn(end)+(-9:0);% 10画素
                nnn(nnn<1)=[];
                d(x,y)=mean(k(nnn));
            else
                d(x,y)=-1;
            end;
        end;
    end;
    figure('color',[1,1,1]);
    d0=max(d(:));
    d(d==-1)=d0+1;
    imagesc(xdir,ydir,d);
    if m==1;
        colormap(map1);
    else
        colormap(map2);
    end;
    daspect([1,1,1]);
    colorbar;
end;
経過時間は 6.363111 秒です。