vb_erosion_2d

 Erosion : 境界から半径 R 以内の近傍点を削除
 外部境界点を抽出してから近傍点操作を行う

 B     : 3D-マスクパターン
 B = 1 : 内部点
   = 0 : 外部点
 R     : 近傍点の半径

 2005-1-18  M. Sato 

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

0001 function    B = vb_erosion_2d(B, R)
0002 % Erosion : 境界から半径 R 以内の近傍点を削除
0003 % 外部境界点を抽出してから近傍点操作を行う
0004 %
0005 % B     : 3D-マスクパターン
0006 % B = 1 : 内部点
0007 %   = 0 : 外部点
0008 % R     : 近傍点の半径
0009 %
0010 % 2005-1-18  M. Sato
0011 %
0012 % Copyright (C) 2011, ATR All Rights Reserved.
0013 % License : New BSD License(see VBMEG_LICENSE.txt)
0014 
0015 [N1,N2] = size(B);
0016 
0017 Bo = zeros(N1,N2);    % (外部点) outer point mask
0018 Bx = zeros(N1,N2);    % (内部境界点) inside boundary point mask
0019 
0020 % neighbor index list
0021 % (X-axis) 隣接点インデックス
0022 j1d = 1:(N1-1);
0023 j1u = 2:N1;
0024 % (Y-axis) 隣接点インデックス
0025 j2d = 1:(N2-1);
0026 j2u = 2:N2;
0027 
0028 % Make outer point mask (外部点)
0029 iz       = find( B == 0 );
0030 Bo(iz) = 1;
0031 
0032 % 内部点マスクを前後左右上下(26近傍点方向)に動かす
0033 Bx(j1d, : ) = Bx(j1d, : ) + B(j1u, : );
0034 Bx(j1u, : ) = Bx(j1u, : ) + B(j1d, : );
0035 Bx( : ,j2d) = Bx( : ,j2d) + B( : ,j2u);
0036 Bx( : ,j2u) = Bx( : ,j2u) + B( : ,j2d);
0037 
0038 % 外部点マスクと内部点マスクを前後左右上下に動かしたものとの積
0039 Bx = Bx.*Bo;
0040 
0041 % 外部境界点の抽出
0042 iz = find( Bx > 0 );
0043 Nz = length(iz);
0044 
0045 % 3D-subscript of external point
0046 [j1,j2] = ind2sub([N1,N2], iz );
0047 
0048 % Make neighbor index within R
0049 Rmax = ceil(R);
0050 x     = -Rmax:Rmax;
0051 
0052 [x1, x2 ] = ndgrid(x);
0053 
0054 jx = find( (x1.^2 + x2.^2 ) <= R^2 );
0055 
0056 x1 = x1(jx);
0057 x2 = x2(jx);
0058 
0059 NN = length(jx);
0060 
0061 % 内部境界点 index
0062 k1 = zeros(Nz,1);
0063 k2 = zeros(Nz,1);
0064 
0065 for n=1:NN,
0066     % Make mask point subscript with fixed neighbor index
0067     k1 = j1 + x1(n);
0068     k2 = j2 + x2(n);
0069 
0070     % Check index limit
0071     k1 = max(k1,1);
0072     k2 = max(k2,1);
0073     
0074     k1 = min(k1,N1);
0075     k2 = min(k2,N2);
0076     
0077     % Add neighbor points to mask
0078     inx = k1+ N1*(k2 - 1);
0079     
0080     B(inx) = 0;
0081 end;