vb_dilation_3d

 Dilation : Add boundary neighbor points with radius R
  Bout = vb_dilation_3d(B, R)

 境界から半径 R 以内の近傍点を加える
 内部境界点を抽出してから近傍点操作を行う

 Bout  : Output mask pattern
 B     : 3D-マスクパターン (3D-mask pattern)
 B = 1 : 内部点 (interior point)
   = 0 : 外部点 (outer point)
 R     : 近傍点の半径 (Radius of eraser)

 2005-1-18  M. Sato 
 内部点マスクを前後左右上下に動かし、内部境界点を抽出
 
 2005-7-2   M. Sato 
 Z方向にループを回し、2D イメージ処理をするように変更

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

0001 function    Bout = vb_dilation_3d(B, R)
0002 % Dilation : Add boundary neighbor points with radius R
0003 %  Bout = vb_dilation_3d(B, R)
0004 %
0005 % 境界から半径 R 以内の近傍点を加える
0006 % 内部境界点を抽出してから近傍点操作を行う
0007 %
0008 % Bout  : Output mask pattern
0009 % B     : 3D-マスクパターン (3D-mask pattern)
0010 % B = 1 : 内部点 (interior point)
0011 %   = 0 : 外部点 (outer point)
0012 % R     : 近傍点の半径 (Radius of eraser)
0013 %
0014 % 2005-1-18  M. Sato
0015 % 内部点マスクを前後左右上下に動かし、内部境界点を抽出
0016 %
0017 % 2005-7-2   M. Sato
0018 % Z方向にループを回し、2D イメージ処理をするように変更
0019 %
0020 % Copyright (C) 2011, ATR All Rights Reserved.
0021 % License : New BSD License(see VBMEG_LICENSE.txt)
0022 
0023 if vb_matlab_version >= 7
0024     Bout = vb_dilation_3d_int(B, R);
0025     return;
0026 end
0027 
0028 [N1,N2,N3] = size(B);
0029 
0030 Bout = B;    % Output mask
0031 
0032 % Make outer point mask (外部点)
0033 Bo     = zeros(N1,N2,N3);    
0034 iz       = find( B == 0 );
0035 Bo(iz) = 1;
0036 
0037 % neighbor index list for boundary detection
0038 % (X-axis) 隣接点インデックス
0039 j1d = 1:(N1-1);
0040 j1u = 2:N1;
0041 % (Y-axis) 隣接点インデックス
0042 j2d = 1:(N2-1);
0043 j2u = 2:N2;
0044 
0045 % Make neighbor index within R
0046 Rmax  = ceil(R);
0047 x      = -Rmax:Rmax;
0048 
0049 [x1, x2, x3 ] = ndgrid(x);
0050 
0051 jx = find( (x1.^2 + x2.^2 +x3.^2 ) <= R^2 );
0052 
0053 x1 = x1(jx);
0054 x2 = x2(jx);
0055 x3 = x3(jx);
0056 
0057 NN   = length(jx);
0058 N1N2 = N1*N2;
0059 
0060 Bx   = zeros(N1,N2);    % (内部境界点) inside boundary point mask
0061 
0062 for zz = 1:N3
0063     
0064     Bx = zeros(N1,N2);    % (内部境界点) inside boundary point mask
0065     
0066     % 外部点マスクを前後左右上下(6近傍点方向)に動かす
0067     
0068     Bx(j1d, : ) = Bx(j1d, : ) + Bo(j1u, : , zz );
0069     Bx(j1u, : ) = Bx(j1u, : ) + Bo(j1d, : , zz );
0070     Bx( : ,j2d) = Bx( : ,j2d) + Bo( : ,j2u, zz );
0071     Bx( : ,j2u) = Bx( : ,j2u) + Bo( : ,j2d, zz );
0072     
0073     if zz > 1,
0074         Bx = Bx + Bo(:,:,zz-1);
0075     end
0076     if zz < N3,
0077         Bx = Bx + Bo(:,:,zz+1);
0078     end
0079     
0080     % 内部点マスクと外部点マスクを前後左右に動かしたものとの積
0081     Bx = Bx.*B(:,:,zz);
0082     
0083     % 内部境界点の抽出
0084     iz = find( Bx > 0 );
0085     Nz = length(iz);
0086     
0087     % 2D-subscript of inner boundary point
0088     % ----- [j1,j2] = ind2sub([N1,N2], iz );
0089     j3  = zz;
0090     j2  = floor((iz-1)/N1)+1;
0091     j1  = rem((iz-1),N1)+1;
0092     
0093     % 内部境界点 index
0094     k1 = zeros(Nz,1);
0095     k2 = zeros(Nz,1);
0096     
0097     for n=1:NN,
0098         % Make mask point subscript with fixed neighbor index
0099         k1 = j1 + x1(n);
0100         k2 = j2 + x2(n);
0101         k3 = j3 + x3(n);
0102     
0103         % Check index limit
0104         k1 = max(k1,1);
0105         k2 = max(k2,1);
0106         k3 = max(k3,1);
0107         
0108         k1 = min(k1,N1);
0109         k2 = min(k2,N2);
0110         k3 = min(k3,N3);
0111         
0112         % Add neighbor points to mask
0113         %---- ind2sub ----
0114         inx = k1+ N1*(k2 - 1)+ N1N2*(k3 - 1);
0115         
0116         Bout(inx) = 1;
0117     end;
0118 end