surf2mesh

 [node,elem,face]=surf2mesh(v,f,p0,p1,keepratio,maxvol,regions,holes,forcebox)

 create quality volumetric mesh from isosurface patches

 author: Qianqian Fang, <q.fang at neu.edu>
 date: 2007/11/24

 input parameters:
      v: input, isosurface node list, dimension (nn,3)
         if v has 4 columns, the last column specifies mesh density near each node
      f: input, isosurface face element list, dimension (be,3)
      p0: input, coordinates of one corner of the bounding box, p0=[x0 y0 z0]
      p1: input, coordinates of the other corner of the bounding box, p1=[x1 y1 z1]
      keepratio: input, percentage of elements being kept after the simplification
      maxvol: input, maximum tetrahedra element volume
      regions: list of regions, specifying by an internal point for each region
      holes: list of holes, similar to regions
      forcebox: 1: add bounding box, 0: automatic

 outputs:
      node: output, node coordinates of the tetrahedral mesh
      elem: output, element list of the tetrahedral mesh
      face: output, mesh surface element list of the tetrahedral mesh 
             the last column denotes the boundary ID

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

0001 function [node,elem,face]=surf2mesh(v,f,p0,p1,keepratio,maxvol,regions,holes,forcebox,method)
0002 %
0003 % [node,elem,face]=surf2mesh(v,f,p0,p1,keepratio,maxvol,regions,holes,forcebox)
0004 %
0005 % create quality volumetric mesh from isosurface patches
0006 %
0007 % author: Qianqian Fang, <q.fang at neu.edu>
0008 % date: 2007/11/24
0009 %
0010 % input parameters:
0011 %      v: input, isosurface node list, dimension (nn,3)
0012 %         if v has 4 columns, the last column specifies mesh density near each node
0013 %      f: input, isosurface face element list, dimension (be,3)
0014 %      p0: input, coordinates of one corner of the bounding box, p0=[x0 y0 z0]
0015 %      p1: input, coordinates of the other corner of the bounding box, p1=[x1 y1 z1]
0016 %      keepratio: input, percentage of elements being kept after the simplification
0017 %      maxvol: input, maximum tetrahedra element volume
0018 %      regions: list of regions, specifying by an internal point for each region
0019 %      holes: list of holes, similar to regions
0020 %      forcebox: 1: add bounding box, 0: automatic
0021 %
0022 % outputs:
0023 %      node: output, node coordinates of the tetrahedral mesh
0024 %      elem: output, element list of the tetrahedral mesh
0025 %      face: output, mesh surface element list of the tetrahedral mesh
0026 %             the last column denotes the boundary ID
0027 %
0028 % -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
0029 %
0030 
0031 fprintf(1,'generating tetrahedral mesh from closed surfaces ...\n');
0032 
0033 if(nargin<10)
0034    method = 'tetgen';
0035 end
0036 
0037 exesuff=getexeext;
0038 exesuff=fallbackexeext(exesuff,method);
0039 
0040 if(keepratio>1 | keepratio<0)
0041    warn(['The "keepratio" parameter is required to be between 0 and 1. '...
0042          'Your input is out of this range. surf2mesh will not perform '...
0043      'simplification. Please double check to correct this.']);
0044 end
0045 
0046 % first, resample the surface mesh with cgal
0047 if(keepratio<1-1e-9 & ~iscell(f))
0048     fprintf(1,'resampling surface mesh ...\n');
0049     [no,el]=meshresample(v(:,1:3),f(:,1:3),keepratio);
0050     el=unique(sort(el,2),'rows');
0051 
0052     % then smooth the resampled surface mesh (Laplace smoothing)
0053 
0054     %% edges=surfedge(el);  % disable on 12/05/08, very slow on octave
0055     %% mask=zeros(size(no,1),1);
0056     %% mask(unique(edges(:)))=1;  % =1 for edge nodes, =0 otherwise
0057     %[conn,connnum,count]=meshconn(el,length(no));
0058     %no=smoothsurf(no,mask,conn,2);
0059 
0060     % remove end elements (all nodes are edge nodes)
0061     %el=delendelem(el,mask);
0062 else
0063     no=v;
0064     el=f;
0065 end
0066 if(nargin==6)
0067     regions=[];
0068     holes=[];
0069 elseif(nargin==7)
0070     holes=[];
0071 end
0072 
0073 if(size(regions,2)>=4 && ~isempty(maxvol))
0074     warning('you specified both maxvol and the region based volume constraint,the maxvol setting will be ignored');
0075     maxvol=[];
0076 end
0077 
0078 dobbx=0;
0079 if(nargin>=9)
0080     dobbx=forcebox;
0081 end
0082 
0083 % dump surface mesh to .poly file format
0084 if(~iscell(el) & ~isempty(no) & ~isempty(el))
0085     saveoff(no(:,1:3),el(:,1:3),mwpath('post_vmesh.off'));
0086 end
0087 deletemeshfile(mwpath('post_vmesh.mtr'));
0088 savesurfpoly(no,el,holes,regions,p0,p1,mwpath('post_vmesh.poly'),dobbx);
0089 
0090 moreopt='';
0091 if(size(no,2)==4)
0092    moreopt=[moreopt ' -m '];
0093 end
0094 % call tetgen to create volumetric mesh
0095 deletemeshfile(mwpath('post_vmesh.1.*'));
0096 fprintf(1,'creating volumetric mesh from a surface mesh ...\n');
0097 
0098 try
0099     cmdopt=evalin('caller','ISO2MESH_TETGENOPT');
0100 catch
0101     try
0102         cmdopt=evalin('base','ISO2MESH_TETGENOPT');
0103     catch
0104         cmdopt='';
0105     end
0106 end
0107 if(isempty(cmdopt))
0108   system([' "' mcpath(method,exesuff) '" -A -q1.414a' num2str(maxvol) ' ' moreopt ' "' mwpath('post_vmesh.poly') '"']);
0109 else
0110   system([' "' mcpath(method,exesuff) '" ' cmdopt ' "' mwpath('post_vmesh.poly') '"']);
0111 end
0112 
0113 % read in the generated mesh
0114 [node,elem,face]=readtetgen(mwpath('post_vmesh.1'));
0115 
0116 fprintf(1,'volume mesh generation is complete\n');
0117