vb_erosion_3d_int

 Erosion : Erase boundary points with radius R
 Use int8-type for MATLAB ver.7

 境界から半径 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_erosion_3d_int(B, R)
0002 % Erosion : Erase boundary points with radius R
0003 % Use int8-type for MATLAB ver.7
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 % 整数型データタイプ宣言
0024 B = int8(B);
0025 
0026 [N1,N2,N3] = size(B);
0027 
0028 Bout = B;    % Output mask
0029 
0030 % Make outer point mask (外部点)
0031 % 整数型データタイプ宣言
0032 Bo     = zeros(N1,N2,N3,'int8');
0033 iz       = find( B == 0 );
0034 Bo(iz) = 1;
0035 
0036 % neighbor index list
0037 % (X-axis) 隣接点インデックス
0038 j1d = 1:(N1-1);
0039 j1u = 2:N1;
0040 % (Y-axis) 隣接点インデックス
0041 j2d = 1:(N2-1);
0042 j2u = 2:N2;
0043 
0044 % Make neighbor index within R
0045 Rmax = ceil(R);
0046 x     = -Rmax:Rmax;
0047 
0048 [x1, x2, x3 ] = ndgrid(x);
0049 
0050 jx = find( (x1.^2 + x2.^2 +x3.^2 ) <= R^2 );
0051 x1 = x1(jx);
0052 x2 = x2(jx);
0053 x3 = x3(jx);
0054 
0055 NN   = length(jx);
0056 N1N2 = N1*N2;
0057 
0058 Bx = zeros(N1,N2, 'int8');    % (内部境界点) inside boundary point mask
0059 
0060 for zz = 1:N3
0061     
0062     Bx = zeros(N1,N2, 'int8');
0063     
0064     % 内部点マスクを前後左右上下(26近傍点方向)に動かす。
0065     
0066     Bx(j1d, : ) = Bx(j1d, : ) + B(j1u, : , zz );
0067     Bx(j1u, : ) = Bx(j1u, : ) + B(j1d, : , zz );
0068     Bx( : ,j2d) = Bx( : ,j2d) + B( : ,j2u, zz );
0069     Bx( : ,j2u) = Bx( : ,j2u) + B( : ,j2d, zz );
0070     
0071     if zz > 1,
0072         Bx = Bx + B(:,:,zz-1);
0073     end
0074     if zz < N3,
0075         Bx = Bx + B(:,:,zz+1);
0076     end
0077     
0078     % 内部点マスクと外部点マスクを前後左右に動かしたものとの積
0079     Bx = Bx.*Bo(:,:,zz);
0080     
0081     % 内部境界点の抽出
0082     iz = find( Bx > 0 );
0083     Nz = length(iz);
0084     
0085     % 2D-subscript of inner boundary point
0086     % ----- [j1,j2] = ind2sub([N1,N2], iz );
0087     j3  = zz;
0088     j2  = floor((iz-1)/N1)+1;
0089     j1  = rem((iz-1),N1)+1;
0090     
0091     % 内部境界点 index
0092     k1 = zeros(Nz,1);
0093     k2 = zeros(Nz,1);
0094     
0095     for n=1:NN,
0096         % Make mask point subscript with fixed neighbor index
0097         k1 = j1 + x1(n);
0098         k2 = j2 + x2(n);
0099         k3 = j3 + x3(n);
0100     
0101         % Check index limit
0102         k1 = max(k1,1);
0103         k2 = max(k2,1);
0104         k3 = max(k3,1);
0105         
0106         k1 = min(k1,N1);
0107         k2 = min(k2,N2);
0108         k3 = min(k3,N3);
0109         
0110         % Delete neighbor points from mask
0111         %---- ind2sub ----
0112         inx = k1+ N1*(k2 - 1)+ N1N2*(k3 - 1);
0113         
0114         Bout(inx) = 0;
0115     end;
0116 end