cgalv2m

 [node,elem,face]=cgalv2m(vol,opt,maxvol)

 wrapper for CGAL 3D mesher (CGAL 3.5 or up)
 convert a binary (or multi-valued) volume to tetrahedral mesh

 http://www.cgal.org/Manual/3.5/doc_html/cgal_manual/Mesh_3/Chapter_main.html

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

 input:
     vol: a volumetric binary image
     ix,iy,iz: subvolume selection indices in x,y,z directions
     opt: parameters for CGAL mesher, if opt is a structure, then
         opt.radbound: defines the maximum surface element size
         opt.angbound: defines the miminum angle of a surface triangle
         opt.distbound: defines the maximum distance between the 
         center of the surface bounding circle and center of the 
         element bounding sphere
         opt.reratio:  maximum radius-edge ratio
         if opt is a scalar, it only specifies radbound.
     maxvol: target maximum tetrahedral elem volume

 output:
     node: output, node coordinates of the tetrahedral mesh
     elem: output, element list of the tetrahedral mesh, the last 
          column is the region id
     face: output, mesh surface element list of the tetrahedral mesh
          the last column denotes the boundary ID
          note: each triangle will appear twice in the face list with each
            one attaches to each side of the interface. one can remove
            the redundant triangles by unique(face(:,1:3),'rows')

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

0001 function [node,elem,face]=cgalv2m(vol,opt,maxvol)
0002 %
0003 % [node,elem,face]=cgalv2m(vol,opt,maxvol)
0004 %
0005 % wrapper for CGAL 3D mesher (CGAL 3.5 or up)
0006 % convert a binary (or multi-valued) volume to tetrahedral mesh
0007 %
0008 % http://www.cgal.org/Manual/3.5/doc_html/cgal_manual/Mesh_3/Chapter_main.html
0009 %
0010 % author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
0011 %
0012 % input:
0013 %     vol: a volumetric binary image
0014 %     ix,iy,iz: subvolume selection indices in x,y,z directions
0015 %     opt: parameters for CGAL mesher, if opt is a structure, then
0016 %         opt.radbound: defines the maximum surface element size
0017 %         opt.angbound: defines the miminum angle of a surface triangle
0018 %         opt.distbound: defines the maximum distance between the
0019 %         center of the surface bounding circle and center of the
0020 %         element bounding sphere
0021 %         opt.reratio:  maximum radius-edge ratio
0022 %         if opt is a scalar, it only specifies radbound.
0023 %     maxvol: target maximum tetrahedral elem volume
0024 %
0025 % output:
0026 %     node: output, node coordinates of the tetrahedral mesh
0027 %     elem: output, element list of the tetrahedral mesh, the last
0028 %          column is the region id
0029 %     face: output, mesh surface element list of the tetrahedral mesh
0030 %          the last column denotes the boundary ID
0031 %          note: each triangle will appear twice in the face list with each
0032 %            one attaches to each side of the interface. one can remove
0033 %            the redundant triangles by unique(face(:,1:3),'rows')
0034 %
0035 % -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
0036 %
0037 
0038 fprintf(1,'creating surface and tetrahedral mesh from a multi-domain volume ...\n');
0039 
0040 dtype=class(vol);
0041 if(~(islogical(vol) | strcmp(dtype,'uint8')))
0042     error('cgalmesher can only handle uint8 volumes, you have to convert your image to unit8 first.');
0043 end
0044 
0045 if(~any(vol))
0046         error('no labeled regions found in the input volume.');
0047 end
0048 
0049 exesuff=getexeext;
0050 exesuff=fallbackexeext(exesuff,'cgalmesh');
0051 
0052 ang=30;
0053 ssize=6;
0054 approx=0.5;
0055 reratio=3;
0056 
0057 if(~isstruct(opt))
0058     ssize=opt;
0059 end
0060 
0061 if(isstruct(opt) & length(opt)==1)  % does not support settings for multiple labels
0062     if(isfield(opt,'radbound'))   ssize=opt.radbound; end
0063     if(isfield(opt,'angbound'))   ang=opt.angbound; end
0064     if(isfield(opt,'distbound')) approx=opt.distbound; end
0065     if(isfield(opt,'reratio'))    reratio=opt.reratio; end
0066 end
0067 
0068 saveinr(vol,mwpath('pre_cgalmesh.inr'));
0069 deletemeshfile(mwpath('post_cgalmesh.mesh'));
0070 
0071 randseed=hex2dec('623F9A9E'); % "U+623F U+9A9E"
0072 
0073 if(~isempty(getvarfrom('base','ISO2MESH_RANDSEED')))
0074         randseed=getvarfrom('base','ISO2MESH_RANDSEED');
0075 end
0076 
0077 if(ischar(maxvol))
0078     format_maxvol='%s';
0079 else
0080     format_maxvol='%f';
0081 end
0082 cmd=sprintf(['"%s%s" "%s" "%s" %f %f %f %f ' format_maxvol ' %d'],mcpath('cgalmesh'),exesuff,...
0083     mwpath('pre_cgalmesh.inr'),mwpath('post_cgalmesh.mesh'),ang,ssize,...
0084     approx,reratio,maxvol,randseed);
0085 system(cmd);
0086 if(~exist(mwpath('post_cgalmesh.mesh'),'file'))
0087     error(['output file was not found, failure was encountered when running command: \n',cmd]);
0088 end
0089 [node,elem,face]=readmedit(mwpath('post_cgalmesh.mesh'));
0090 
0091 % if a transformation matrix/offset vector supplied, apply them
0092 if (isstruct(opt) & length(opt)==1)
0093     if(isfield(opt,'A') & isfield(opt,'B'))
0094         node(:,1:3)=(opt.A*node(:,1:3)'+repmat(opt.B(:),1,size(node,1)))';
0095     end
0096 end
0097 
0098 fprintf(1,'node number:\t%d\ntriangles:\t%d\ntetrahedra:\t%d\nregions:\t%d\n',...
0099     size(node,1),size(face,1),size(elem,1),length(unique(elem(:,end))));
0100 fprintf(1,'surface and volume meshes complete\n');
0101 
0102 if(size(node,1)>0)
0103     [node,elem,face]=sortmesh(node(1,:),node,elem,1:4,face,1:3);
0104 end
0105 
0106 node=node+0.5;