vb_get_boundary

 get boundary surface point index of mask image
  [xyz, BZ] = vb_get_boundary(B,vlevel)
   B      : 3D-mask image
 vlevel   : Threshold
 xyz : voxel coordinate of boundary
   BZ     : boundary mask image

 2007/06/14 Masa-aki Sato

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

0001 function    [xyz, BZ] = vb_get_boundary(B,vlevel)
0002 % get boundary surface point index of mask image
0003 %  [xyz, BZ] = vb_get_boundary(B,vlevel)
0004 %   B      : 3D-mask image
0005 % vlevel   : Threshold
0006 % xyz : voxel coordinate of boundary
0007 %   BZ     : boundary mask image
0008 %
0009 % 2007/06/14 Masa-aki Sato
0010 %
0011 % Copyright (C) 2011, ATR All Rights Reserved.
0012 % License : New BSD License(see VBMEG_LICENSE.txt)
0013 
0014 if ~exist('vlevel','var'), vlevel = 0.5; end;
0015 
0016 [N1,N2,N3] = size(B);
0017 
0018 BZ = zeros(N1,N2,N3);
0019 
0020 % West-East (X-axis) 隣接点インデックス
0021 NN1 = N1-1;
0022 jw  = 1:NN1;
0023 je  = 2:N1;
0024 % North-South (Y-axis) 隣接点インデックス
0025 NN2 = N2-1;
0026 jn  = 1:NN2;
0027 js  = 2:N2;
0028 
0029 NX2 = N1*N2;
0030 
0031 for j3 = 1:N3
0032     % West-East (X-axis) 隣接点比較
0033     ix = find( (B(jw,:,j3)-vlevel) .* (B(je,:,j3)-vlevel) <= 0 );
0034     
0035     % 2D-subscript of inner boundary point
0036     % ----- [j1,j2] = ind2sub([NN1,N2], ix );
0037     j2 = floor((ix-1)/NN1)+1;
0038     j1 = rem((ix-1),NN1)+1;
0039     jx = j1+ N1*(j2 - 1)+ NX2*(j3 - 1);
0040     
0041     % set flag on the boundary points
0042     BZ(jx) = 1;
0043 
0044     % North-South (Y-axis) 隣接点比較
0045     ix = find( (B(:,jn,j3)-vlevel) .* (B(:,js,j3)-vlevel) <= 0 );
0046 
0047     % ----- [j1,j2] = ind2sub([N1,N2], ix );
0048     j2 = floor((ix-1)/N1)+1;
0049     j1 = rem((ix-1),N1)+1;
0050     jx = j1+ N1*(j2 - 1)+ NX2*(j3 - 1);
0051     
0052     % set flag on the boundary points
0053     BZ(jx) = 1;
0054 
0055     % Up-Down (Z-axis) 隣接点比較
0056     if j3 < N3,
0057         ix = find( (B(:,:,j3)-vlevel) .* (B(:,:,j3+1)-vlevel) <= 0 );
0058 
0059         % ----- [j1,j2] = ind2sub([N1,N2], ix );
0060         j2 = floor((ix-1)/N1)+1;
0061         j1 = rem((ix-1),N1)+1;
0062         jx = j1+ N1*(j2 - 1)+ NX2*(j3 - 1);
0063         
0064         % set flag on the boundary points
0065         BZ(jx) = 1;
0066     end
0067 end
0068 
0069 ix = find( BZ(:) > 0 );
0070 
0071 % ----- [j1,j2,j3] = ind2sub([N1,N2,N3], ix );
0072 j3  = floor((ix-1)/NX2)+1;
0073 ix  = rem((ix-1),NX2) +1;
0074 
0075 j2 = floor((ix-1)/N1)+1;
0076 j1 = rem((ix-1),N1)+1;
0077 
0078 xyz = [j1, j2, j3];