mesh2vol

 [mask weight]=mesh2vol(node,face,Nxyz)
 [mask weight]=mesh2vol(node,face,[Nx,Ny,Nz])
 [mask weight]=mesh2vol(node,face,xi,yi,zi,hf)
   or
 newval=mesh2vol(node_val,face,...)

 fast rasterization of a 3D mesh to a volume with tetrahedron index labels
 
 author: Qianqian Fang <fangq at nmr.mgh.harvard.edu>
 date for initial version: Feb 10,2014

 input:
      node: node coordinates, dimension N by 2 or N by 3 array
      nodeval: a 4-column array, the first 3 columns are the node coordinates, 
            the last column denotes the values associated with each node
      face: a triangle surface, N by 3 or N by 4 array
      Nx,Ny,Nxy: output image in x/y dimensions, or both
      xi,yi: linear vectors for the output pixel center positions in x/y
      hf: the handle of a pre-created figure window for faster rendering

 output:
      mask: a 3D image, the value of each pixel is the index of the
            enclosing triangle, if the pixel is outside of the mesh, NaN
      weight: (optional) a 3 by Nx by Ny array, where Nx/Ny are the dimensions
            for the mask
      newval: when the node has 4 columns, the last column represents the
            values (color) at each node, the output newval is the rasterized
            mesh value map over the specified grid.

 note: This function only works for matlab

 example:

   [no,el]=meshgrid6(0:5,0:5,0:3);
   mask=mesh2vol(no,el(:,1:4),0:0.1:5,0:0.1:5,0:0.1:4);
   imagesc(mask(:,:,3))

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

0001 function [mask weight]=mesh2vol(node,elem,xi,yi,zi)
0002 %
0003 % [mask weight]=mesh2vol(node,face,Nxyz)
0004 % [mask weight]=mesh2vol(node,face,[Nx,Ny,Nz])
0005 % [mask weight]=mesh2vol(node,face,xi,yi,zi,hf)
0006 %   or
0007 % newval=mesh2vol(node_val,face,...)
0008 %
0009 % fast rasterization of a 3D mesh to a volume with tetrahedron index labels
0010 %
0011 % author: Qianqian Fang <fangq at nmr.mgh.harvard.edu>
0012 % date for initial version: Feb 10,2014
0013 %
0014 % input:
0015 %      node: node coordinates, dimension N by 2 or N by 3 array
0016 %      nodeval: a 4-column array, the first 3 columns are the node coordinates,
0017 %            the last column denotes the values associated with each node
0018 %      face: a triangle surface, N by 3 or N by 4 array
0019 %      Nx,Ny,Nxy: output image in x/y dimensions, or both
0020 %      xi,yi: linear vectors for the output pixel center positions in x/y
0021 %      hf: the handle of a pre-created figure window for faster rendering
0022 %
0023 % output:
0024 %      mask: a 3D image, the value of each pixel is the index of the
0025 %            enclosing triangle, if the pixel is outside of the mesh, NaN
0026 %      weight: (optional) a 3 by Nx by Ny array, where Nx/Ny are the dimensions
0027 %            for the mask
0028 %      newval: when the node has 4 columns, the last column represents the
0029 %            values (color) at each node, the output newval is the rasterized
0030 %            mesh value map over the specified grid.
0031 %
0032 % note: This function only works for matlab
0033 %
0034 % example:
0035 %
0036 %   [no,el]=meshgrid6(0:5,0:5,0:3);
0037 %   mask=mesh2vol(no,el(:,1:4),0:0.1:5,0:0.1:5,0:0.1:4);
0038 %   imagesc(mask(:,:,3))
0039 %
0040 % -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
0041 %
0042 
0043 nodeval=[];
0044 if(size(node,2)==4)
0045    nodeval=node(:,4);
0046    node(:,4)=[];
0047 end
0048 
0049 if(nargin==3 && length(xi)==1 && xi>0)
0050     mn=min(node);
0051     mx=max(node);
0052     df=(mx(1:min(3,size(node,2)))-mn(1:min(3,size(node,2))))/xi;
0053 elseif(nargin==3 && length(xi)==3 && all(xi>0))
0054     mn=min(node);
0055     mx=max(node);
0056     df=(mx(1:min(3,size(node,2)))-mn(1:min(3,size(node,2))))./(xi(:)');
0057 elseif(nargin==5)
0058     mx=[max(xi) max(yi) max(zi)];
0059     mn=[min(xi) min(yi) min(zi)];
0060     df=[min(diff(xi(:))) min(diff(yi(:))) min(diff(zi(:)))];
0061 else
0062     error('you must give at least xi input');
0063 end
0064 
0065 xi=mn(1):df(1):mx(1);
0066 yi=mn(2):df(2):mx(2);
0067 zi=mn(3):df(3):mx(3);
0068 
0069 if(size(node,2)~=3 || size(elem,2)<=3)
0070     error('node must have 3 columns; face can not have less than 3 columns');
0071 end
0072 
0073 nz=length(zi);
0074 mask=zeros(length(xi)-1,length(yi)-1,length(zi)-1);
0075 if(nargout>1 || ~isempty(nodeval))
0076    weight=zeros(4,length(xi)-1,length(yi)-1,length(zi)-1);
0077 end
0078 
0079 hf=figure('visible','on');
0080 for i=1:nz
0081     if(~isempty(nodeval))
0082         [cutpos,cutvalue,facedata,elemid]=qmeshcut(elem,node,nodeval,['z=' num2str(zi(i))]);
0083     else
0084         [cutpos,cutvalue,facedata,elemid]=qmeshcut(elem,node,node(:,1),['z=' num2str(zi(i))]);
0085     end
0086     if(isempty(cutpos))
0087         continue;
0088     end
0089     if(nargout>1 || ~isempty(nodeval))
0090         [maskz, weightz]=mesh2mask(cutpos,facedata,xi,yi,hf);
0091         weight(:,:,:,i)=weightz;
0092     else
0093         maskz=mesh2mask(cutpos,facedata,xi,yi,hf);
0094     end
0095     idx=find(~isnan(maskz));
0096     if(~isempty(nodeval))
0097         eid=facedata(maskz(idx),:);
0098         maskz(idx)=cutvalue(eid(:,1)).*weightz(1,idx)'+cutvalue(eid(:,2)).*weightz(2,idx)'+...
0099             cutvalue(eid(:,3)).*weightz(3,idx)'+cutvalue(eid(:,4)).*weightz(4,idx)';
0100     else
0101         maskz(idx)=elemid(maskz(idx));
0102     end
0103     mask(:,:,i)=maskz;
0104 end
0105 close(hf);