vb_aniso_difuse

 anisotropic difusion
 v = vb_aniso_difuse(v,u,Nloop,c)

 2005/2/5 M.Sato

 v : 平滑化するマスクパターン(白質・頭蓋内)
 u : v の相補パターン(灰白質・頭蓋)
 Nloop : 拡散回数
 c     : 拡散係数

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function    v = vb_aniso_difuse(v,u,Nloop,c)
0002 % anisotropic difusion
0003 % v = vb_aniso_difuse(v,u,Nloop,c)
0004 %
0005 % 2005/2/5 M.Sato
0006 %
0007 % v : 平滑化するマスクパターン(白質・頭蓋内)
0008 % u : v の相補パターン(灰白質・頭蓋)
0009 % Nloop : 拡散回数
0010 % c     : 拡散係数
0011 %
0012 % Copyright (C) 2011, ATR All Rights Reserved.
0013 % License : New BSD License(see VBMEG_LICENSE.txt)
0014 
0015 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0016 %      異方性拡散緩和による境界面の平滑化
0017 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0018 
0019 [N1, N2, N3]=size(v);
0020 
0021 siz=[N1 N2 N3];
0022 
0023 % 隣接点インデックスリストの作成
0024 
0025 % West-East (X-axis) 隣接点インデックス
0026 jw = 1:(N1-1);
0027 je = 2:N1;
0028 % North-South (Y-axis) 隣接点インデックス
0029 jn = 1:(N2-1);
0030 js = 2:N2;
0031 % Up-Down (Z-axis) 隣接点インデックス
0032 jd = 1:(N3-1);
0033 ju = 2:N3;
0034 
0035 % 256x256x256　はメモリーエラーになる
0036 % メモリエラーを避けるため、2D-配列計算のループにする
0037 % 2D-配列計算のループの方が3D-配列計算より高速
0038 vw = zeros(N1-1, N2, 1);
0039 vn = zeros(N1, N2-1, 1);
0040 vd = zeros(N1, 1, N3-1);
0041 
0042 % 異方性拡散緩和による境界面の平滑化
0043 
0044 for n=1:Nloop,
0045     for k=1:N3,
0046         % West-east
0047         vw   = ( v(je,:,k) - v(jw,:,k) ).*( u(je,:,k) + u(jw,:,k) )*0.5;
0048         v(jw,:,k) = v(jw,:,k) + c*vw;
0049         v(je,:,k) = v(je,:,k) - c*vw;
0050         % North-south
0051         vn   = ( v(:,js,k) - v(:,jn,k) ).*( u(:,js,k) + u(:,jn,k) )*0.5;
0052         v(:,jn,k) = v(:,jn,k) + c*vn;
0053         v(:,js,k) = v(:,js,k) - c*vn;
0054     end;
0055         
0056         % Doun-up
0057     for k=1:N2,
0058         vd   = ( v(:,k,ju) - v(:,k,jd) ).*( u(:,k,ju) + u(:,k,jd) )*0.5;
0059         v(:,k,jd) = v(:,k,jd) + c*vd;
0060         v(:,k,ju) = v(:,k,ju) - c*vd;
0061     end;
0062     
0063 end;    
0064 
0065