vb_surf_to_mask

 Make surface to mask image
   B = vb_surf_to_mask(V,F,vstep,step,Dim)

 V(NV,3) : vertex point coordinate on the surface
 F(NF,3) : patch index
 vstep   : voxcel size
 step    : minimum length for edge
 Dim = [NX, NY, NZ] : Dimension of mask image

 B(NX, NY, NZ) : mask image

 Made by M. Sato 2004-3-28

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

0001 function    B = vb_surf_to_mask(V,F,vstep,step,Dim)
0002 % Make surface to mask image
0003 %   B = vb_surf_to_mask(V,F,vstep,step,Dim)
0004 %
0005 % V(NV,3) : vertex point coordinate on the surface
0006 % F(NF,3) : patch index
0007 % vstep   : voxcel size
0008 % step    : minimum length for edge
0009 % Dim = [NX, NY, NZ] : Dimension of mask image
0010 %
0011 % B(NX, NY, NZ) : mask image
0012 %
0013 % Made by M. Sato 2004-3-28
0014 %
0015 % Copyright (C) 2011, ATR All Rights Reserved.
0016 % License : New BSD License(see VBMEG_LICENSE.txt)
0017 
0018 
0019 NX    = Dim(1);
0020 NY    = Dim(2);
0021 NZ    = Dim(3);
0022 NXY = NX*NY;
0023 
0024 NV    = size(V,1);
0025 NF    = size(F,1);
0026 
0027 % Voxcel image
0028 B    = zeros(NX, NY, NZ);
0029 V   = V/vstep;
0030 step = step/vstep;
0031 %
0032 % Set voxcel at vertex point 'ON'
0033 % 三角形頂点のボクセルを ON にする
0034 % change coordinate to index
0035 %     V = [-1/2,-1/2,-1/2] - [1/2,1/2,1/2]
0036 % <=> J = [1,1,1]
0037 
0038 %J  = floor(V) + 1;
0039 J  = floor(V + 0.5) + 1;
0040 ix = J(:,1)+ NX*(J(:,2) - 1)+ NXY*(J(:,3) - 1);
0041 
0042 B(ix) = 1;
0043 
0044 %
0045 % Find triangle larger than step size
0046 % 格子サイズより長い辺を持つ三角形を探す
0047 % triangle vertex　三角形頂点
0048 x1 = V(F(:,1),:);
0049 x2 = V(F(:,2),:);
0050 x3 = V(F(:,3),:);
0051 
0052 % triangle edge　三角形辺
0053 z1 = x1 - x2;
0054 z2 = x2 - x3;
0055 z3 = x3 - x1;
0056 
0057 % edge length　辺の長さ
0058 d1 = sqrt(sum(z1.^2,2));
0059 d2 = sqrt(sum(z2.^2,2));
0060 d3 = sqrt(sum(z3.^2,2));
0061 
0062 % Largest edge length in each triangle　最長辺
0063 [dd ,imax] = max([d1 d2 d3],[],2);
0064 
0065 % find edge larger than step　格子サイズより長い辺を探す
0066 idx = find( dd > step );
0067 ND  = length(idx);
0068 
0069 if ND == 0, return; end;
0070 
0071 % Max length
0072 dmax = max(dd);
0073 
0074 % 2D index to specify vertex point inside the triangle
0075 % 三角形内部の分割点を表すインデックス生成
0076 Nmax = ceil(dmax/step);
0077 NN   = (Nmax+1)*(Nmax+2)/2;
0078 inx1 = zeros(NN,1);
0079 inx2 = zeros(NN,1);
0080 id   = 0;
0081 
0082 % 0 <= k2 <= k1 <= Nmax
0083 % 0 <= Nmax - k1 + k2 <= Nmax
0084 for k1=0:Nmax
0085     for k2=0:k1
0086         id = id + 1;
0087         inx1(id) = k1;
0088         inx2(id) = k2;
0089     end;
0090 end
0091 
0092 % large triangle loop
0093 for m = 1:ND
0094     %
0095     n  = idx(m);        % triangle index
0096     d  = dd(n);            % largest edge length
0097     Nd = ceil(d/step);    % # of edge division
0098     
0099     % triangle vertex
0100     y1 = x1(n,:);
0101     y2 = x2(n,:);
0102     y3 = x3(n,:);
0103     
0104     % 0 <= Nd - k1 + k2 <= Nd
0105     % 0 <=  kk1 + kk2   <= Nd
0106     NT  = (Nd+1)*(Nd+2)/2;
0107     kk1 = Nd - inx1(1:NT);    % kk1 = Nd - k1
0108     kk2 = inx2(1:NT);        % kk2 = k2
0109 
0110     % area coordinate index
0111     % kk1 + kk2 + kk3 = Nd
0112     kk3 = Nd - kk1 - kk2;    % kk1 + kk2 + kk3 = Nd
0113     
0114     % triangle inner coordinate
0115     j1 = floor(( kk1*y1(1) + kk2*y2(1) + kk3*y3(1) )/Nd ) + 1;
0116     j2 = floor(( kk1*y1(2) + kk2*y2(2) + kk3*y3(2) )/Nd ) + 1;
0117     j3 = floor(( kk1*y1(3) + kk2*y2(3) + kk3*y3(3) )/Nd ) + 1;
0118 
0119     % 1D-index corresponding to (j1, j2, j3)
0120     ii = j1+ NX*(j2 - 1)+ NXY*(j3 - 1);
0121     
0122     %
0123     B(ii) = 1;
0124     
0125 end
0126