surf2volz

 img=surf2volz(node,face,xi,yi,zi)

 convert a triangular surface to a shell of voxels in a 3D image
 along the z-axis

 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
     xi,yi,zi: x/y/z grid for the resulting volume

 output:
     img: a volumetric binary image at position of ndgrid(xi,yi,zi)

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

0001 function img=surf2volz(node,face,xi,yi,zi)
0002 %
0003 % img=surf2volz(node,face,xi,yi,zi)
0004 %
0005 % convert a triangular surface to a shell of voxels in a 3D image
0006 % along the z-axis
0007 %
0008 % author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
0009 %
0010 % input:
0011 %     node: node list of the triangular surface, 3 columns for x/y/z
0012 %     face: triangle node indices, each row is a triangle
0013 %     xi,yi,zi: x/y/z grid for the resulting volume
0014 %
0015 % output:
0016 %     img: a volumetric binary image at position of ndgrid(xi,yi,zi)
0017 %
0018 % -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
0019 %
0020 
0021 ne=size(face,1);
0022 img=zeros(length(xi),length(yi),length(zi),'uint8');
0023 
0024 dx0=min(abs(diff(xi)));
0025 dx=dx0/2;
0026 dy0=min(abs(diff(yi)));
0027 dy=dy0/2;
0028 dz0=min(abs(diff(zi)));
0029 
0030 dl=sqrt(dx*dx+dy*dy);
0031 
0032 minz=min(node(:,3));
0033 maxz=max(node(:,3));
0034 iz=hist([minz,maxz],zi);
0035 hz=find(iz);
0036 iz=hz(1):min(length(zi),hz(end)+1);
0037 
0038 for i=1:length(iz);
0039     plane=[0 100 zi(iz(i)); 100 0 zi(iz(i)); 0 0 zi(iz(i))];
0040     [bcutpos,bcutvalue,bcutedges]=qmeshcut(face(:,1:3),node,node(:,1),plane);
0041     if(isempty(bcutpos))
0042         continue;
0043     end
0044     enum=length(bcutedges);
0045     for j=1:enum
0046         e0=bcutpos(bcutedges(j,1),1:2);
0047         e1=bcutpos(bcutedges(j,2),1:2);
0048         len=ceil(sum(abs(e1-e0))/(abs(dx)+abs(dy)))+1;
0049         dd=(e1-e0)/len;
0050 
0051         posx= floor(((e0(1)+(0:len)*dd(1)-xi(1)))/dx0)';
0052         posy= floor(((e0(2)+(0:len)*dd(2)-yi(1)))/dy0)';
0053         pos=[posx, posy];
0054         pos(find(posx>length(xi) | posy>length(yi) | posx<=0|posy<=0), :)=[];
0055 
0056         if(length(pos)>0)
0057             zz=floor(((zi(iz(i))-zi(1)))/dz0);
0058             for k=1:size(pos,1)
0059                 img(pos(k,1),pos(k,2),zz)=1;
0060             end
0061           %img(sub2ind(size(img),pos(:,1),pos(:,2),i*ones(size(pos,1),1)))=1;
0062         end
0063     end
0064 end
0065