vb_get_boundary_mesh

 get boundary surface of mask image
  [xyz,F] = vb_get_boundary_mesh(v,vlevel)
 マスク境界面の抽出
    v       : 3D-ボクセル上のマスクイメージ
   vlevel   : 閾値

   xyz      : 境界点座標
   F        : 三角面インデックス

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

0001 function    [xyz,F] = vb_get_boundary_mesh(v,vlevel)
0002 % get boundary surface of mask image
0003 %  [xyz,F] = vb_get_boundary_mesh(v,vlevel)
0004 % マスク境界面の抽出
0005 %    v       : 3D-ボクセル上のマスクイメージ
0006 %   vlevel   : 閾値
0007 %
0008 %   xyz      : 境界点座標
0009 %   F        : 三角面インデックス
0010 %
0011 % Copyright (C) 2011, ATR All Rights Reserved.
0012 % License : New BSD License(see VBMEG_LICENSE.txt)
0013 
0014 if nargin < 2, vlevel = 0.5; end;
0015 
0016 [N1, N2, N3]=size(v);
0017 
0018 siz=[N1 N2 N3];
0019 
0020 % 隣接点インデックスリストの作成
0021 
0022 % West-East (X-axis) 隣接点インデックス
0023 jw = 1:(N1-1);
0024 je = 2:N1;
0025 % North-South (Y-axis) 隣接点インデックス
0026 jn = 1:(N2-1);
0027 js = 2:N2;
0028 % Up-Down (Z-axis) 隣接点インデックス
0029 jd = 1:(N3-1);
0030 ju = 2:N3;
0031 
0032 % Boundary lattice flag
0033 vb = zeros(N1+1, N2+1, N3+1);
0034 vx = zeros(N1+1, N2+1, N3+1);
0035 vy = zeros(N1+1, N2+1, N3+1);
0036 vz = zeros(N1+1, N2+1, N3+1);
0037 
0038 F =[];
0039 % 256x256x256　はメモリーエラーになる
0040 % メモリエラーを避けるため、2D-配列計算のループにする
0041 % 2D-配列計算のループの方が3D-配列計算より高速
0042 vw = zeros(N1+1, N2+1, 1);
0043 vn = zeros(N1+1, N2+1, 1);
0044 vd = zeros(N1+1,    1, N3+1);
0045 
0046 % Extract boundary points
0047 
0048 for n=1:N3,
0049     % West-East (X-axis) 隣接点比較
0050     iw            = find( (v(jw,:,n)-vlevel) .* (v(je,:,n)-vlevel) <= 0 );
0051     
0052     if ~isempty(iw),
0053         % 2D-subscript of boundary point
0054         [j1,j2,j3]    = ind2sub([N1-1,N2,1], iw );
0055         jw1            = sub2ind([N1+1,N2+1], j1+1,j2  );
0056         jw2            = sub2ind([N1+1,N2+1], j1+1,j2+1);
0057         vw            = zeros(N1+1, N2+1, 1);
0058         % boundary surface base point for Y-Z plane
0059         vw(jw1)        = 1;
0060         vx(:,:,n)    = vx(:,:,n) + vw;
0061         % boundary lattice point
0062         vw(jw2)        = vw(jw2) + 1;
0063         vb(:,:,n)    = vb(:,:,n)   + vw;
0064         vb(:,:,n+1) = vb(:,:,n+1) + vw;
0065     end;
0066     
0067     % North-South (Y-axis) 隣接点比較
0068     iw           = find( (v(:,jn,n)-vlevel) .* (v(:,js,n)-vlevel) <= 0 );
0069     
0070     if ~isempty(iw),
0071         % 2D-subscript of boundary point
0072         [j1,j2,j3]    = ind2sub([N1,N2-1,1], iw );
0073         jw1            = sub2ind([N1+1,N2+1], j1,  j2+1);
0074         jw2            = sub2ind([N1+1,N2+1], j1+1,j2+1);
0075         vw            = zeros(N1+1, N2+1, 1);
0076         % boundary surface base point for Z-X plane
0077         vw(jw1)        = 1;
0078         vy(:,:,n)    = vy(:,:,n) + vw;
0079         % boundary lattice point
0080         vw(jw2)        = vw(jw2) + 1;
0081         vb(:,:,n)    = vb(:,:,n)   + vw;
0082         vb(:,:,n+1) = vb(:,:,n+1) + vw;
0083     end;
0084 end;
0085 
0086 for n=1:N2,
0087     % Up-Down (Z-axis) 隣接点比較
0088     id           = find( (v(:,n,jd)-vlevel) .* (v(:,n,ju)-vlevel) <= 0 );
0089     
0090     if ~isempty(id),
0091         % 2D-subscript of boundary point
0092         [j1,j2,j3]    = ind2sub([N1, 1, N3-1], id );
0093         j2(:)        = 1;
0094         jw1            = sub2ind([N1+1,1,N3+1], j1,  j2,j3+1);
0095         jw2            = sub2ind([N1+1,1,N3+1], j1+1,j2,j3+1);
0096         % boundary surface base point for X-Y plane
0097         vd            = zeros(N1+1,1,N3+1);
0098         vd(jw1)        = 1;
0099         vz(:,n,:)    = vz(:,n,:) + vd;
0100         % boundary lattice point
0101         vd(jw2)        = vd(jw2) + 1;
0102         vb(:,n  ,:)    = vb(:,n  ,:) + vd;
0103         vb(:,n+1,:) = vb(:,n+1,:) + vd;
0104     end;
0105 end;
0106 
0107 % Make triangle index
0108 
0109 % Boundary lattice points
0110 ix       = find( vb > 0 );
0111 NB       = length(ix);
0112 
0113 % index transformation table
0114 vb(ix) = 1:NB;
0115 
0116 % Boundary lattice point coordinate
0117 [kx,ky,kz] = ind2sub([N1+1, N2+1, N3+1], ix );
0118 
0119 xyz=[kx-1 ,ky-1 ,kz-1 ];
0120 
0121 % West-East (X-axis) 隣接境界
0122 iw    = find( vx > 0 );
0123 % 3D-subscript of boundary surface base point
0124 [j1,j2,j3]    = ind2sub([N1+1, N2+1, N3+1], iw );
0125 
0126 % other corner points in the square
0127 jw1            = sub2ind([N1+1,N2+1,N3+1], j1,j2+1,j3  );
0128 jw2            = sub2ind([N1+1,N2+1,N3+1], j1,j2  ,j3+1);
0129 jw3            = sub2ind([N1+1,N2+1,N3+1], j1,j2+1,j3+1);
0130 
0131 % add triangle index
0132 F=[F; vb(iw),vb(jw1),vb(jw2) ; vb(jw1),vb(jw2),vb(jw3)];
0133 
0134 % North-South (Y-axis) 隣接境界
0135 iw    = find( vy > 0 );
0136 % 3D-subscript of boundary surface base point
0137 [j1,j2,j3]    = ind2sub([N1+1, N2+1, N3+1], iw );
0138 
0139 % other corner points in the square
0140 jw1            = sub2ind([N1+1,N2+1,N3+1], j1+1,j2,j3  );
0141 jw2            = sub2ind([N1+1,N2+1,N3+1], j1  ,j2,j3+1);
0142 jw3            = sub2ind([N1+1,N2+1,N3+1], j1+1,j2,j3+1);
0143 
0144 % add triangle index
0145 F=[F; vb(iw),vb(jw1),vb(jw2) ; vb(jw1),vb(jw2),vb(jw3)];
0146 
0147 % Up-Down (Z-axis) 隣接境界
0148 iw    = find( vz > 0 );
0149 % 3D-subscript of boundary surface base point
0150 [j1,j2,j3]    = ind2sub([N1+1, N2+1, N3+1], iw );
0151 
0152 % other corner points in the square
0153 jw1            = sub2ind([N1+1,N2+1,N3+1], j1+1,j2  ,j3);
0154 jw2            = sub2ind([N1+1,N2+1,N3+1], j1  ,j2+1,j3);
0155 jw3            = sub2ind([N1+1,N2+1,N3+1], j1+1,j2+1,j3);
0156 
0157 % add triangle index
0158 F=[F; vb(iw),vb(jw1),vb(jw2) ; vb(jw1),vb(jw2),vb(jw3)];