vb_dilation_3d_int

 Dilation : Add boundary neighbor 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_dilation_3d_int(B, R)
0002 % Dilation : Add boundary neighbor points with radius R
0003 % Use int8-type for MATLAB ver.7
0004 % 境界から半径 R 以内の近傍点を加える
0005 % 内部境界点を抽出してから近傍点操作を行う
0006 %
0007 % Bout  : Output mask pattern
0008 % B     : 3D-マスクパターン (3D-mask pattern)
0009 % B = 1 : 内部点 (interior point)
0010 %   = 0 : 外部点 (outer point)
0011 % R     : 近傍点の半径 (Radius of eraser)
0012 %
0013 % 2005-1-18  M. Sato
0014 % 内部点マスクを前後左右上下に動かし、内部境界点を抽出
0015 %
0016 % 2005-7-2   M. Sato
0017 % Z方向にループを回し、2D イメージ処理をするように変更
0018 %
0019 % Copyright (C) 2011, ATR All Rights Reserved.
0020 % License : New BSD License(see VBMEG_LICENSE.txt)
0021 
0022 % 整数型データタイプ宣言
0023 B = int8(B);
0024 
0025 [N1,N2,N3] = size(B);
0026 
0027 Bout = B;    % Output mask
0028 
0029 % Make outer point mask (外部点)
0030 % 整数型データタイプ宣言
0031 Bo     = zeros(N1,N2,N3,'int8');
0032 iz       = find( B == 0 );
0033 Bo(iz) = 1;
0034 
0035 % neighbor index list for boundary detection
0036 % (X-axis) 隣接点インデックス
0037 j1d = 1:(N1-1);
0038 j1u = 2:N1;
0039 % (Y-axis) 隣接点インデックス
0040 j2d = 1:(N2-1);
0041 j2u = 2:N2;
0042 
0043 % Make neighbor index within R
0044 Rmax  = ceil(R);
0045 x      = -Rmax:Rmax;
0046 
0047 [x1, x2, x3 ] = ndgrid(x);
0048 
0049 jx = find( (x1.^2 + x2.^2 +x3.^2 ) <= R^2 );
0050 
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     % 外部点マスクを前後左右上下(6近傍点方向)に動かす
0065     
0066     Bx(j1d, : ) = Bx(j1d, : ) + Bo(j1u, : , zz );
0067     Bx(j1u, : ) = Bx(j1u, : ) + Bo(j1d, : , zz );
0068     Bx( : ,j2d) = Bx( : ,j2d) + Bo( : ,j2u, zz );
0069     Bx( : ,j2u) = Bx( : ,j2u) + Bo( : ,j2d, zz );
0070     
0071     if zz > 1,
0072         Bx = Bx + Bo(:,:,zz-1);
0073     end
0074     if zz < N3,
0075         Bx = Bx + Bo(:,:,zz+1);
0076     end
0077     
0078     % 内部点マスクと外部点マスクを前後左右に動かしたものとの積
0079     Bx = Bx.*B(:,:,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         % Add neighbor points to mask
0111         %---- ind2sub ----
0112         inx = k1+ N1*(k2 - 1)+ N1N2*(k3 - 1);
0113         
0114         Bout(inx) = 1;
0115     end;
0116 end