surf2vol

 [img, v2smap]=surf2vol(node,face,xi,yi,zi,'options',values,...)

 convert a triangular surface to a shell of voxels in a 3D image

 author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)

 input:
     node: node list of the triangular surface, 3 columns for x/y/z
     face: triangle node indices, each row is a triangle
              if face contains the 4th column, it indicates the label of
              the face triangles (each face componment must be closed); if
              face contains 5 columns, it stores a tetrahedral mesh with
              labels, where the first 4 columns are the element list and 
              the last column is the element label;
     xi,yi,zi: x/y/z grid for the resulting volume
        options: 'fill', if set to 1, the enclosed voxels are labeled by 1
                 'label', if set to 1, the enclosed voxels are labeled by
                          the corresponding label of the face or element;
                          setting 'label' to 1 also implies 'fill'.

 output:
     img: a volumetric binary image at position of ndgrid(xi,yi,zi)
        v2smap (optional): a 4x4 matrix denoting the Affine transformation to map
             the voxel coordinates back to the mesh space.

 -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)

0001 function [img, v2smap]=surf2vol(node,face,xi,yi,zi,varargin)
0002 %
0003 % [img, v2smap]=surf2vol(node,face,xi,yi,zi,'options',values,...)
0004 %
0005 % convert a triangular surface to a shell of voxels in a 3D image
0006 %
0007 % author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
0008 %
0009 % input:
0010 %     node: node list of the triangular surface, 3 columns for x/y/z
0011 %     face: triangle node indices, each row is a triangle
0012 %              if face contains the 4th column, it indicates the label of
0013 %              the face triangles (each face componment must be closed); if
0014 %              face contains 5 columns, it stores a tetrahedral mesh with
0015 %              labels, where the first 4 columns are the element list and
0016 %              the last column is the element label;
0017 %     xi,yi,zi: x/y/z grid for the resulting volume
0018 %        options: 'fill', if set to 1, the enclosed voxels are labeled by 1
0019 %                 'label', if set to 1, the enclosed voxels are labeled by
0020 %                          the corresponding label of the face or element;
0021 %                          setting 'label' to 1 also implies 'fill'.
0022 %
0023 % output:
0024 %     img: a volumetric binary image at position of ndgrid(xi,yi,zi)
0025 %        v2smap (optional): a 4x4 matrix denoting the Affine transformation to map
0026 %             the voxel coordinates back to the mesh space.
0027 %
0028 % -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
0029 %
0030 
0031 fprintf(1,'converting a closed surface to a volumetric binary image ...\n');
0032 
0033 opt=varargin2struct(varargin{:});
0034 label=jsonopt('label',0,opt);
0035 
0036 elabel=[1];
0037 if(size(face,2)>=4)
0038         elabel=unique(face(:,end));
0039         if(size(face,2)==5)
0040             label=1;
0041             el=face;
0042             face=[];
0043             for i=1:length(elabel)
0044                     fc=volface(el(el(:,5)==elabel(i),1:4));
0045                     fc(:,4)=elabel(i);
0046                     face=[face ; fc];
0047             end
0048         end
0049 else
0050         fc=face;
0051 end
0052 
0053 img=zeros(length(xi),length(yi),length(zi));
0054 
0055 for i=1:length(elabel)
0056         if(size(face,2)==4)
0057             fc=face(face(:,4)==elabel(i),1:3);
0058         end
0059         im=surf2volz(node(:,1:3),fc(:,1:3),xi,yi,zi);
0060         im=im | shiftdim(surf2volz(node(:,[3 1 2]),fc(:,1:3),zi,xi,yi),1);
0061         im=im | shiftdim(surf2volz(node(:,[2 3 1]),fc(:,1:3),yi,zi,xi),2);
0062 
0063         v2smap=[];
0064 
0065         % here we assume the grid is uniform; surf2vol can handle non-uniform grid,
0066         % but the affine output is not correct in this case
0067 
0068         if(jsonopt('fill',0,opt) || label)
0069                 im=imfill(im,'holes');
0070                 if(label)
0071                     im=im*elabel(i);
0072                 end
0073         end
0074         img=max(im,img);
0075 end
0076  
0077 if(nargout>1) 
0078         dlen=abs([xi(2)-xi(1) yi(2)-yi(1) zi(2)-zi(1)]);
0079         p0=min(node);
0080         offset=p0;
0081         v2smap=diag(abs(dlen));
0082         v2smap(4,4)=1;
0083         v2smap(1:3,4)=offset';
0084 end