mesh2mask

 [mask weight]=mesh2mask(node,face,Nxy)
   or
 [mask weight]=mesh2mask(node,face,[Nx,Ny])
   or
 [mask weight]=mesh2mask(node,face,xi,yi,hf)

 fast rasterization of a 2D mesh to an image with triangle index labels
 
 author: Qianqian Fang <fangq at nmr.mgh.harvard.edu>
 date for initial version: July 18,2013

 input:
      node: node coordinates, dimension N by 2 or N by 3 array
      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: (optional) the handle of a pre-created figure window, for faster 
          rendering

 output:
      mask: a 2D 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

 note: This function only works in MATLAB when the DISPLAY is not 
       disabled. The maximum size of the mask output is limited by the 
       screen size.

 example:

   [no,fc]=meshgrid6(0:5,0:5);
   [mask weight]=mesh2mask(no,fc,-1:0.1:5,0:0.1:5);
   imagesc(mask);

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

0001 function [mask, weight]=mesh2mask(node,face,xi,yi,hf)
0002 %
0003 % [mask weight]=mesh2mask(node,face,Nxy)
0004 %   or
0005 % [mask weight]=mesh2mask(node,face,[Nx,Ny])
0006 %   or
0007 % [mask weight]=mesh2mask(node,face,xi,yi,hf)
0008 %
0009 % fast rasterization of a 2D mesh to an image with triangle index labels
0010 %
0011 % author: Qianqian Fang <fangq at nmr.mgh.harvard.edu>
0012 % date for initial version: July 18,2013
0013 %
0014 % input:
0015 %      node: node coordinates, dimension N by 2 or N by 3 array
0016 %      face: a triangle surface, N by 3 or N by 4 array
0017 %      Nx,Ny,Nxy: output image in x/y dimensions, or both
0018 %      xi,yi: linear vectors for the output pixel center positions in x/y
0019 %      hf: (optional) the handle of a pre-created figure window, for faster
0020 %          rendering
0021 %
0022 % output:
0023 %      mask: a 2D image, the value of each pixel is the index of the
0024 %            enclosing triangle, if the pixel is outside of the mesh, NaN
0025 %      weight: (optional) a 3 by Nx by Ny array, where Nx/Ny are the dimensions for
0026 %            the mask
0027 %
0028 % note: This function only works in MATLAB when the DISPLAY is not
0029 %       disabled. The maximum size of the mask output is limited by the
0030 %       screen size.
0031 %
0032 % example:
0033 %
0034 %   [no,fc]=meshgrid6(0:5,0:5);
0035 %   [mask weight]=mesh2mask(no,fc,-1:0.1:5,0:0.1:5);
0036 %   imagesc(mask);
0037 %
0038 % -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
0039 %
0040 
0041 if(nargin==3 && length(xi)==1 && xi>0)
0042     mn=min(node);
0043     mx=max(node);
0044     df=(mx(1:2)-mn(1:2))/xi;
0045 elseif(nargin==3 && length(xi)==2 && all(xi>0))
0046     mn=min(node);
0047     mx=max(node);
0048     df=(mx(1:2)-mn(1:2))./xi;
0049 elseif(nargin==4 || nargin==5)
0050     mx=[max(xi) max(yi)];
0051     mn=[min(xi) min(yi)];
0052     df=[min(diff(xi(:))) min(diff(yi(:)))];
0053 else
0054     error('you must give at least xi input');
0055 end
0056 if(size(node,2)<=1 || size(face,2)<=2)
0057     error('node must have 2 or 3 columns; face can not have less than 2 columns');
0058 end
0059 
0060 if(nargin<5)
0061     hf=figure('visible','on');
0062 else
0063     clf(hf);
0064 end
0065 patch('Vertices',node,'Faces',face,'linestyle','none','FaceColor','flat',...
0066  'FaceVertexCData',(1:size(face,1))','CDataMapping', 'scaled');
0067 set(gca, 'Position', [0 0 1 1]);
0068 cm=jet(size(face,1));
0069 colormap(cm);
0070 axis off
0071 set(gca,'xlim',[mn(1) mx(1)]);
0072 set(gca,'ylim',[mn(2) mx(2)]);
0073 set(gca,'clim',[1 size(face,1)]);
0074 
0075 output_size = round((mx(1:2)-mn(1:2))./df);%Size in pixels
0076 
0077 if(isoctavemesh || isempty(getenv('DISPLAY')))
0078     resolution = 300; %Resolution in DPI
0079     set(gcf,'PaperPositionMode','manual')
0080     set(gcf,'paperunits','inches','paperposition',[0 0 output_size/resolution]);
0081     deletemeshfile(mwpath('post_mesh2mask.png'));
0082     print(mwpath('post_mesh2mask.png'),'-dpng',['-r' num2str(resolution)]);
0083     mask=imread(mwpath('post_mesh2mask.png'));
0084 else
0085     pos=get(hf,'position');
0086     pos(3:4)=max(pos(3:4),output_size+20);
0087     set(hf,'position',pos);
0088     set(gca, 'Units','pixels','position',[1, 1, output_size(1), output_size(2)]);
0089     mask=getframe(gca);
0090     if(any(size(mask.cdata)<[output_size([2 1]) 3]))
0091         error('the requested rasterization grid is larger than the screen resolution');
0092     end
0093     mask=mask.cdata(1:output_size(2),1:output_size(1),:);
0094 end
0095 if(nargin<5)
0096     close(hf);
0097 end
0098 mask=int32(reshape(mask,[size(mask,1)*size(mask,2) size(mask,3)]));
0099 [isfound,locb]=ismember(mask,floor(cm*255),'rows');
0100 locb(isfound==0)=nan;
0101 
0102 mask=rot90(reshape(locb,output_size([2 1]))');
0103 
0104 if(nargout>=2)
0105     xi=mn(1)+df(1)/2:df(1):mx(1);
0106     yi=mn(2)+df(2)/2:df(2):mx(2);
0107     weight=barycentricgrid(node,face,xi,yi,mask);
0108     if(size(face,2)>=4)
0109         badidx=find(weight(1,:,:)<0 | weight(2,:,:)<0 | weight(3,:,:)<0);
0110         badidx=badidx(face(mask(badidx),3)~=face(mask(badidx),4));
0111         weight2=barycentricgrid(node,face(:,[1 3 4]),xi,yi,mask);
0112         weight(:,badidx)=0;
0113         weight([1 3 4],badidx)=weight2(:,badidx);
0114     end
0115 end
0116 
0117 function weight=barycentricgrid(node,face,xi,yi,mask)
0118 [xx,yy]=meshgrid(xi,yi);
0119 idx=find(~isnan(mask));
0120 eid=mask(idx);
0121 t1=node(face(:,1),:);
0122 t2=node(face(:,2),:);
0123 t3=node(face(:,3),:);
0124 tt=(t2(:,2)-t3(:,2)).*(t1(:,1)-t3(:,1))+(t3(:,1)-t2(:,1)).*(t1(:,2)-t3(:,2));
0125 w(:,1)=(t2(eid,2)-t3(eid,2)).*(xx(idx)-t3(eid,1))+(t3(eid,1)-t2(eid,1)).*(yy(idx)-t3(eid,2));
0126 w(:,2)=(t3(eid,2)-t1(eid,2)).*(xx(idx)-t3(eid,1))+(t1(eid,1)-t3(eid,1)).*(yy(idx)-t3(eid,2));
0127 w(:,1)=w(:,1)./tt(eid);
0128 w(:,2)=w(:,2)./tt(eid);
0129 w(:,3)=1-w(:,1)-w(:,2);
0130 weight=zeros(3,size(mask,1),size(mask,2));
0131 ww=zeros(size(mask));
0132 ww(idx)=w(:,1);
0133 weight(1,:,:)=ww;
0134 ww(idx)=w(:,2);
0135 weight(2,:,:)=ww;
0136 ww(idx)=w(:,3);
0137 weight(3,:,:)=ww;