=== # Streamlined mesh generation - shortcuts ===

==== function [node,elem,face]=v2m(img,isovalues,opt,maxvol,method) ====
 [node,elem,face]=v2m(img,isovalues,opt,maxvol,method)
 volumetric mesh generation from binary or gray-scale volumetric images
 shortcut for vol2mesh
 inputs and outputs are similar to those defined in vol2mesh

==== function [no,el,regions,holes]=v2s(img,isovalues,opt,method) ====
 [no,el,regions,holes]=v2s(img,isovalues,opt,method)
 surface mesh generation from binary or gray-scale volumetric images
 shortcut for vol2surf
 inputs and outputs are similar to those defined in vol2surf; In v2s, 
 method can be set to 'cgalmesh' in addition to those allowed by vol2surf.

==== function [node,elem,face]=s2m(v,f,keepratio,maxvol,method) ====
 [node,elem,face]=s2m(v,f,keepratio,maxvol)
 volumetric mesh generation from a closed surface, shortcut for surf2mesh
 inputs and outputs are similar to those defined in surf2mesh
 if method='cgalpoly', s2m will call cgals2m and keepratio should be a 
 structure (as the 'opt' input in cgals2m)

==== function img=s2v(node,face,div) ====
 img=s2v(node,face,div)
 shortcut for surf2vol, coverting a surface to a volumetric image
 input:
	 node: node list of the triangular surface, 3 columns for x/y/z
	 face: triangle node indices, each row is a triangle
	 div:  division number along the shortest edge of the mesh (resolution)
              if not given, div=50
 output:
	 img: a volumetric binary image at position of ndgrid(xi,yi,zi)

==== function newnode=sms(node,face,iter,alpha,method) ====
 newnode=sms(node,face,iter,useralpha,method)
 simplified version of surface mesh smoothing
 input:
    node:  node coordinates of a surface mesh
    face:  face element list of the surface mesh
    iter:  smoothing iteration number
    alpha: scaler, smoothing parameter, v(k+1)=alpha*v(k)+(1-alpha)*mean(neighbors)
    method: same as in smoothsurf, default is 'laplacianhc'
 output:
    newnode: output, the smoothed node coordinates
=== # Streamlined mesh generation ===

==== function [node,elem,face,regions]=vol2mesh(img,ix,iy,iz,opt,maxvol,dofix,method,isovalues) ====
 [node,elem,face,regions]=vol2mesh(img,ix,iy,iz,opt,maxvol,dofix,method,isovalues)
 convert a binary (or multi-valued) volume to tetrahedral mesh
 input:
	 img: a volumetric binary image
	 ix,iy,iz: subvolume selection indices in x,y,z directions
	 opt: as defined in vol2surf.m
	 maxvol: target maximum tetrahedral elem volume
                when method='cgalmesh', maxvol can specify the target
                for each label (subregion index) by the following syntax
                'label1=size1:label2=size2:...'
	 dofix: 1: perform mesh validation&repair, 0: skip repairing
	 method: 'cgalsurf' or omit: use CGAL surface mesher
		 'simplify': use binsurface and then simplify
		 'cgalmesh': use CGAL 3.5 3D mesher for direct mesh generation [new]
		 generally speaking, 'cgalmesh' is the most robust path
		 if you want to product meshes from binary or multi-region
		 volumes, however, its limitations include 1) only accept 
		 uint8 volume, and 2) can not extract meshes from gray-scale
		 volumes. If ones goal is to process a gray-scale volume,
		 he/she should use the 'cgalsurf' option. 'simplify' approach
		 is not recommended unless other options has failed.
	 isovalues: a list of isovalues where the levelset is defined
 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
    region: optional output. if opt.autoregion is set to 1, region
          saves the interior points for each closed surface component

==== function [no,el,regions,holes]=vol2surf(img,ix,iy,iz,opt,dofix,method,isovalues) ====
 [no,el,regions,holes]=vol2surf(img,ix,iy,iz,opt,dofix,method,isovalues)
 converting a 3D volumetric image to surfaces
 input:
	 img: a volumetric binary image; if img is empty, vol2surf will
	      return user defined surfaces via opt.surf if it exists
	 ix,iy,iz: subvolume selection indices in x,y,z directions
	 opt: function parameters
	   if method is 'cgalsurf' or 'cgalpoly':
	     opt=a float number>1: max radius of the Delaunay sphere(element size) 
	     opt.radbound: same as above, max radius of the Delaunay sphere
	     opt.distbound: maximum deviation from the specified isosurfaces
	     opt(1,2,...).radbound: same as above, for each levelset
	   if method is 'simplify':
	     opt=a float number<1: compression rate for surf. simplification
	     opt.keepratio=a float less than 1: same as above, same for all surf.
	     opt(1,2,..).keepratio: setting compression rate for each levelset
	   opt(1,2,..).maxsurf: 1 - only use the largest disjointed surface
				0 - use all surfaces for that levelset
          opt(1,2,..).side: - 'upper': threshold at upper interface
                              'lower': threshold at lower interface
	   opt(1,2,..).maxnode: - the maximum number of surface node per levelset
	   opt(1,2,..).holes: user specified holes interior pt list
	   opt(1,2,..).regions: user specified regions interior pt list
	   opt(1,2,..).surf.{node,elem}: add additional surfaces
	   opt(1,2,..).{A,B}: linear transformation for each surface
	   opt.autoregion: if set to 1, vol2surf will try to determine 
              the interior points for each closed surface automatically
	 dofix: 1: perform mesh validation&repair, 0: skip repairing
	 method: - if method is 'simplify', iso2mesh will first call
		   binsurface to generate a voxel-based surface mesh and then
		   use meshresample/meshcheckrepair to create a coarser mesh;
		 - if method is 'cgalsurf', iso2mesh will call the surface
		   extraction program from CGAL to make surface mesh
		 - if method is not specified, 'cgalsurf' is assumed by default
	 isovalues: a list of isovalues where the levelset is defined
 output: 
	 no:  list of nodes on the resulting suface mesh, 3 columns for x,y,z
	 el:  list of trianglular elements on the surface, [n1,n2,n3,region_id]
	 regions: list of interior points for all sub-region, [x,y,z]
	 holes:   list of interior points for all holes, [x,y,z]

==== function [node,elem,face]=surf2mesh(v,f,p0,p1,keepratio,maxvol,regions,holes,forcebox) ====
 [node,elem,face]=surf2mesh(v,f,p0,p1,keepratio,maxvol,regions,holes,forcebox)
 create quality volumetric mesh from isosurface patches
 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

==== function img=surf2vol(node,face,xi,yi,zi) ====
 img=surf2vol(node,face,xi,yi,zi)
 convert a triangular surface to a shell of voxels in a 3D image
 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)
=== # iso2mesh main function backend ===

==== function [node,elem]=binsurface(img,nface) ====
 [node,elem]=binsurface(img,nface)
 fast isosurface extraction from 3D binary images
 input: 
   img:  a 3D binary image
   nface: nface=3 or ignored - for triangular faces, 
          nface=4 - square faces
          nface=0 - return a boundary mask image via node
 output:
   elem: integer array with dimensions of NE x nface, each row represents
         a surface mesh face element 
   node: node coordinates, 3 columns for x, y and z respectively
 the outputs of this subroutine can be easily plotted using 
     patch('Vertices',node,'faces',elem,'FaceVertexCData',node(:,3),
           'FaceColor','interp');
 if the surface mesh has triangular faces, one can plot it with
     trisurf(elem,node(:,1),node(:,2),node(:,3))

==== function [node,elem,face]=cgalv2m(vol,opt,maxvol) ====
 [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
 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')

==== function [node,elem,face]=cgals2m(v,f,opt,maxvol,varargin) ====
 [node,elem,face]=cgals2m(v,f,opt,maxvol)
 wrapper for CGAL 3D mesher (CGAL 3.5 and newer)
 convert a triangular surface to tetrahedral mesh
 http://www.cgal.org/Manual/3.5/doc_html/cgal_manual/Mesh_3/Chapter_main.html
 input:
	 v: the node coordinate list of a surface mesh (nn x 3)
	 f: the face element list of a surface mesh (be x 3)
	 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

==== function [node,elem]=vol2restrictedtri(vol,thres,cent,brad,ang,radbound,distbound,maxnode) ====
 [node,elem]=vol2restrictedtri(vol,thres,cent,brad,ang,radbound,distbound,maxnode)
 surface mesh extraction using CGAL mesher
 input:
       vol: a 3D volumetric image
       thres: a scalar as the threshold of of the extraction
       cent: a 3d position (x,y,z) which locates inside the resulting
             mesh, this is automatically computed from vol2surf
       brad: maximum bounding sphere squared of the resulting mesh
       ang: minimum angular constrains of the resulting tranglar elements
            (in degrees)
       radbound: maximum triangle delaunay circle radius
       distbound: maximum delaunay sphere distances
       maxnode: maximum number of surface nodes (even radbound is not reached)
 output:
       node: the list of 3d nodes in the resulting surface (x,y,z)
       elem: the element list of the resulting mesh (3 columns of integers)

==== function img=surf2volz(node,face,xi,yi,zi) ====
 img=surf2volz(node,face,xi,yi,zi)
 convert a triangular surface to a shell of voxels in a 3D image
 along the z-axis
 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)
=== # iso2mesh primitive meshing functions ===

==== function [node,face,elem]=meshabox(p0,p1,opt,nodesize) ====
 [node,face,elem]=meshabox(p0,p1,opt,maxvol)
 create the surface and tetrahedral mesh of a box geometry
 input: 
   p0:  coordinates (x,y,z) for one end of the box diagnoal
   p1:  coordinates (x,y,z) for the other end of the box diagnoal
   opt: maximum volume of the tetrahedral elements
   nodesize: 1 or a 8x1 array, size of the element near each vertex
 output:
   node: node coordinates, 3 columns for x, y and z respectively
   face: integer array with dimensions of NB x 3, each row represents
         a surface mesh face element 
   elem: integer array with dimensions of NE x 4, each row represents
         a tetrahedron 
 example:
   [node,face,elem]=meshabox([2 3 2],[6 12 15],0.1,1);
   plotmesh(node,elem,'x>4');

==== function [node,face,elem]=meshasphere(c0,r,tsize,maxvol) ====
 [node,face,elem]=meshasphere(c0,r,tsize,maxvol)
 create the surface and tetrahedral mesh of a sphere
 input: 
   c0:  center coordinates (x0,y0,z0) of the sphere
   r:   radius of the sphere
   tsize: maximum surface triangle size on the sphere
   maxvol: maximu volume of the tetrahedral elements
 output:
   node: node coordinates, 3 columns for x, y and z respectively
   face: integer array with dimensions of NB x 3, each row represents
         a surface mesh face element 
   elem: integer array with dimensions of NE x 4, each row represents
         a tetrahedron 

==== function [node,face,elem]=meshanellip(c0,rr,tsize,maxvol) ====
 [node,face,elem]=meshanellip(c0,rr,opt)
 create the surface and tetrahedral mesh of an ellipsoid
 input: 
   c0:  center coordinates (x0,y0,z0) of the ellipsoid
   rr:  radii of an ellipsoid, 
        if rr is a scalar, this is a sphere with radius rr
        if rr is a 1x3 or 3x1 vector, it specifies the ellipsoid radii [a,b,c]
        if rr is a 1x5 or 5x1 vector, it specifies [a,b,c,theta,phi]
           where theta and phi are the rotation angles along z and x 
           axes, respectively. Rotation is applied before translation.
   tsize: maximum surface triangle size on the sphere
   maxvol: maximu volume of the tetrahedral elements
 output:
   node: node coordinates, 3 columns for x, y and z respectively
   face: integer array with dimensions of NB x 3, each row represents
         a surface mesh face element 
   elem: integer array with dimensions of NE x 4, each row represents
         a tetrahedron; if ignored, only produces the surface
 example:
   [node,face,elem]=meshanellip([10,10,-10],[30,20,10,pi/4,pi/4],0.5,0.4);
   plotmesh(node,elem,'x>10');axis equal;

==== function [node,face,elem]=meshunitsphere(tsize,maxvol) ====
 [node,face,elem]=meshunitsphere(tsize,maxvol)
 create the surface and/or volumetric mesh of a unit sphere 
 centered at [0 0 0] and radius 1
 input: 
   tsize: maximum size of the surface triangles (from 0 to 1)
   maxvol: maximum volume of the tetrahedron; if one wants to return
           elem without specifying maxvol, maxvol=tsize^3
 output:
   node: node coordinates, 3 columns for x, y and z respectively
   face: integer array with dimensions of NB x 3, each row represents
         a surface mesh face element 
   elem: integer array with dimensions of NE x 4, each row represents
         a tetrahedron. If ignored, this function only produces the surface
 example:
   [node,face]=meshunitsphere(0.05);
   [node,face,elem]=meshunitsphere(0.05,0.01);
   plotmesh(node,elem,'x>0'); axis equal;

==== function [node,face,elem]=meshacylinder(c0,c1,r,tsize,maxvol,ndiv) ====
 [node,face]=meshacylinder(c0,c1,r,tsize,maxvol,ndiv)
    or
 [node,face,elem]=meshacylinder(c0,c1,r,tsize,maxvol,ndiv)
 create the surface and (optionally) tetrahedral mesh of a 3D cylinder
 input: 
   c0, c1:  cylinder axis end points
   r:   radius of the cylinder
   tsize: maximum surface triangle size on the sphere
   maxvol: maximu volume of the tetrahedral elements
   ndiv: approximate the cylinder surface into ndiv flat pieces, if 
         ignored, ndiv is set to 20
 output:
   node: node coordinates, 3 columns for x, y and z respectively
   face: integer array with dimensions of NB x 3, each row represents
         a surface mesh triangle 
   elem: (optional) integer array with dimensions of NE x 4, each row 
         represents a tetrahedron 
=== # Mesh decomposition and query ===

==== function facecell=finddisconnsurf(f) ====
 facecell=finddisconnsurf(f)
 subroutine to extract disconnected surfaces from a 
 cluster of surfaces
 Date: 2008/03/06
 input: 
     f: faces defined by node indices for all surface triangles
 output:
     facecell: separated disconnected surface node indices

==== function openedge=surfedge(f) ====
 openedge=surfedge(f)
 find the edge of an open surface or surface of a volume
 input:
      f: input, surface face element list, dimension (be,3)
 output:
      openedge: list of edges of the specified surface

==== function openface=volface(t) ====
 openface=volface(t)
 find the surface patches of a volume
 input:
      t: input, volumetric element list, dimension (ne,4)
 output:
      openface: list of faces of the specified volume

==== function loops=extractloops(edges) ====
 loops=extractloops(edges)
 extract individual loops from an edge table of a loop
 collection
 input:   
    edges:  two column matrix recording the starting/ending 
             points of all edge segments
 output:
    loops:  output, a single vector separated by NaN, each segment
             is a close-polygon consisted by node IDs

==== function [conn,connnum,count]=meshconn(elem,nn) ====
 [conn,connnum,count]=meshconn(elem,nn)
 create node neighbor list from a mesh
 input:
    elem:  element table of a mesh
    nn  :  total node number of the mesh
 output:
    conn:  output, a cell structure of length nn, conn{n}
           contains a list of all neighboring node ID for node n
    connnum: vector of length nn, denotes the neighbor number of each node
    count: total neighbor numbers

==== function centroid=meshcentroid(v,f) ====
 centroid=meshcentroid(v,f)
 compute the centroids of a mesh defined by nodes and elements
 (surface or tetrahedra) in R^n space
 input:
      v: surface node list, dimension (nn,3)
      f: surface face element list, dimension (be,3)
 output:
      centroid: centroid positions, one row for each element

==== function nodevol=nodevolume(node,elem) ====
 nodevol=nodevolume(node,elem)
 calculate the Voronoi volume of each node in a simplex mesh
 input:
    node:  node coordinates
    elem:  element table of a mesh
 output:
    nodevol:   volume values for all nodes

==== function vol=elemvolume(node,elem,option) ====
 vol=elemvolume(node,elem,option)
 calculate the volume for a list of simplexes
 input:
    node:  node coordinates
    elem:  element table of a mesh
    option: if option='signed', the volume is the raw determinant,
            else, the results will be the absolute values
 output:
    vol:   volume values for all elements

==== function [conn,connnum,count]=neighborelem(elem,nn); ====
 [conn,connnum,count]=neighborelem(elem,nn)
 create node neighbor list from a mesh
 input:
    elem:  element table of a mesh
    nn  :  total node number of the mesh
 output:
    conn:  output, a cell structure of length nn, conn{n}
           contains a list of all neighboring elem ID for node n
    connnum: vector of length nn, denotes the neighbor number of each node
    count: total neighbor numbers

==== function facenb=faceneighbors(t,opt) ====
 facenb=faceneighbors(t,opt)
 to find 4 face-neighboring elements of a tetrahedron
 input:
     t: tetrahedron element list, 4 columns of integers
     opt: if opt='surface', return boundary triangle list 
          (should be the same as the face output from v2m)
          otherwise, return the element list for each element:
          each row contains 4 numbers, representing the element
          indices sharing triangular faces [1 2 3],[1 2 4],[1 3 4]
          and [2 3 4] in order, where 1~4 is the node local index.
          if the index is 0, indicating the face has no neighbor
          (i.e. a boundary face)
 output:
     facenb: see opt

==== function edgenb=edgeneighbors(t,opt) ====
 edgenb=edgeneighbors(t,opt)
 to find neighboring triangular elements in a triangule surface
 input:
     t: a triangular surface element list, 3 columns of integers
     opt: if opt='general', return the edge neighbors for a general
          triangular surface: each edge can be shared by more than 2
          triangles; if ignored, we assume all triangles are shared by no
          more than 2 triangles.
 output:
     edgenb: if opt is not supplied, edgenb is a size(t,1) by 3 array with
     each element being the triangle ID of the edge neighbor of that
     triangle. For each row, the order of the neighbors is listed as those
     sharing edges [1 2], [2 3] and [3 1] between the triangle nodes.
     when opt='general', edgenb is a cell array with a length of size(t).
     each member of the cell array is a list of edge neighbors (the order 
     is not defined).

==== function [f maxsize]=maxsurf(facecell,node) ====
 [f maxsize]=maxsurf(facecell,node)
 return the surface with the maximum element number or  
 total area from a cell arry of surfaces
 input:
    facecell: a cell array, each element is a face array
    node: optional, node list, if given, the output is the
          surface with the largest surface area.
 output:
    f: the surface data (node indices) for the surface with the 
       most elements (or largest area when node is given)
    maxsize: if node is not given, maxisize is row number of f;
             otherwise, maxsize is the total area of f

==== function mask=flatsegment(node,edge) ====
 mask=flatsegment(node,edge)
 decompose edge loops into flat segments alone arbitrary planes of the bounding box
 this code is fragile: it can not handle curves with many co-linear
 nodes near the corner point
 input:   
    node:  x,y,z coordinates of each node of the mesh
    edge:  input, a single vector separated by NaN, each segment
           is a close-polygon consisted by node IDs 
 output:
    mask:  output, a cell, each element is a close-polygon 
           on x/y/z plane 

==== function newedge=orderloopedge(edge) ====
 [newedge]=orderloopedge(edge)
 order the node list of a simple loop based on connection sequence
 input: 
        edge: a loop consisted by a sequence of edges, each row 
              is an edge with two integers: start/end node index
 output:
        newedge: reordered edge node list

==== function [X,V,E,F]=mesheuler(face) ====
 [X,V,E,F]=mesheuler(face)
 Euler's charastistics of a mesh
 input: 
   face: a closed surface mesh
 output:
   X: Euler's charastistics
   V: number of vertices 
   E: number of edges
   F: number of faces

==== function seg=bbxflatsegment(node,loop) ====
 seg=bbxflatsegment(node,loop)
 decompose edge loops into flat segments along the x/y/z 
 planes of the bounding box
 input:   
    node:  x,y,z coordinates of each node of the mesh
    loop:  input, a single vector separated by NaN, each segment
             is a close-polygon consisted by node IDs 
 output:
    seg:   output, a single vector separated by NaN, each segment
             is a close-polygon on x/y/z plane 

==== function plane=surfplane(node,face) ====
 plane=surfplane(node,face)
 plane equation coefficients for each face in a surface
 input:
   node: a list of node coordinates (nn x 3)
   face: a surface mesh triangle list (ne x 3)
 output:
   plane: a (ne x 4) array, in each row, it has [a b c d]
        to denote the plane equation as "a*x+b*y+c*z+d=0"

==== function [pt,p0,v0,t,idx]=surfinterior(node,face) ====
 [pt,p0,v0,t,idx]=surfinterior(node,face)
 identify a point that is enclosed by the (closed) surface
 input:
   node: a list of node coordinates (nn x 3)
   face: a surface mesh triangle list (ne x 3)
 output:
   pt: the interior point coordinates [x y z]
   p0: ray origin used to determine the interior point
   v0: the vector used to determine the interior point
   t : ray-tracing intersection distances (with signs) from p0. the
       intersection coordinates can be expressed as p0+t(i)*v0
   idx: index to the face elements that intersect with the ray, order
       match that of t

==== function elist=surfpart(f,loopedge) ====
 elist=surfpart(f,loopedge)
 partition a triangular surface using a closed loop defined by existing edges
 input:
      f: input, surface face element list, dimension (be,3)
      loopedge: a 2-column array, specifying a closed loop in CCW order
 output:
      elist: list of triangles that is enclosed by the loop

==== function seeds=surfseeds(node,face) ====
 seeds=surfseeds(node,face)
 calculate a set of interior points with each enclosed by a closed
 component of a surface
 input:
   node: a list of node coordinates (nn x 3)
   face: a surface mesh triangle list (ne x 3)
 output:
   seeds: the interior points coordinates for each closed-surface
          component

==== function quality=meshquality(node,elem) ====
 quality=meshquality(node,elem)
 compute the Joe-Liu mesh quality measure of a tetrahedral mesh
 input:
    node:  node coordinates of the mesh (nn x 3)
    elem:  element table of a tetrahedral mesh (ne x 4)
 output:
    quality: a vector of the same length as size(elem,1), with 
           each element being the Joe-Liu mesh quality metric (0-1) of 
           the corresponding element. A value close to 1 represents
           higher mesh quality (1 means equilateral tetrahedron); 
           a value close to 0 means nearly degenerated element.
 reference:
    A. Liu, B. Joe, Relationship between tetrahedron shape measures, 
                    BIT 34 (2) (1994) 268-287.

==== function edges=meshedge(elem,varargin) ====
 edges=meshedge(elem,opt)
 return all edges in a surface or volumetric mesh
 input:
    elem:  element table of a mesh (support N-d space element)
    opt: optional input, giving the additional options. If opt
         is a struct, it can have the following field:
       opt.nodeorder: if 1, assuming the elem node indices is in CCW 
                      orientation; 0 use nchoosek() output to order edges
         you can replace opt by a series of ('param', value) pairs.
 output:
    edge:  edge list; each row is an edge, specified by the starting and
           ending node indices, the total edge number is
           size(elem,1) x nchoosek(size(elem,2),2). All edges are ordered
           by looping through each element first. 

==== function snorm=surfnorm(node,face,varargin) ====
 snorm=surfnorm(node,face)
    or
 snorm=surfnorm(node,face,'Normalize',0)
 compute the normal vectors for a triangular surface
 input:
   node: a list of node coordinates (nn x 3)
   face: a surface mesh triangle list (ne x 3)
   opt: a list of optional parameters, currently surfnorm supports:
        'Normalize': [1|0] if set to 1, the normal vectors will be 
                           unitary (default)
 output:
   snorm: output surface normal vector at each face

==== function [edges,idx,edgemap]=uniqedges(elem) ====
 [edges,idx,edgemap]=uniqedges(elem)
 return the unique edge list from a surface or tetrahedral mesh
 input:
     elem: a list of elements, each row is a list of nodes for an element.
           elem can have 2, 3 or 4 columns
 output:
     edge: unique edges in the mesh, denoted by a pair of node indices
     idx:  index of the output in the raw edge list (returned by meshedge)
     edgemap: index of the raw edges in the output list (for triangular mesh)

==== function [elist,nextfront]=advancefront(edges,loop,elen) ====
 [elist,nextfront]=advancefront(edges,loop,elen)
 advance an edge-front on an oriented surface to the next separated by
 one-element width
 input:
      edges: edge list of an oriented surface mesh, must be in CCW order
      loop: a 2-column array, specifying a closed loop in CCW order
      elen: node number inside each element, if ignored, elen is set to 3
 output:
      elist: list of triangles that is enclosed between the two
             edge-fronts
      nextfront: a new edge loop list representing the next edge-front
=== # Mesh processing and reparing ===

==== function [node,elem]=meshcheckrepair(node,elem,opt,varargin) ====
 [node,elem]=meshcheckrepair(node,elem,opt)
 check and repair a surface mesh
 input/output:
      node: input/output, surface node list, dimension (nn,3)
      elem: input/output, surface face element list, dimension (be,3)
      opt: options, including
            'dupnode': remove duplicated nodes
            'dupelem' or 'duplicated': remove duplicated elements
            'dup': both above
            'isolated': remove isolated nodes
            'open': abort when open surface is found
            'deep': call external jmeshlib to remove non-manifold vertices
            'meshfix': repair a closed surface by the meshfix utility (new)
                       it can remove self-intersecting elements and fill holes
            'intersect': test a surface for self-intersecting elements

==== function newelem=meshreorient(node,elem) ====
 newelem=meshreorient(node,elem)
 reorder nodes in a surface or tetrahedral mesh to ensure all
 elements are oriented consistently
 input:
    node: list of nodes
    elem: list of elements (each row are indices of nodes of each element)
 output:
    newelem: the element list with consistent ordering

==== function elem=removedupelem(elem) ====
 elem=removedupelem(elem)
 remove doubly duplicated (folded) elements
 input:
    elem: list of elements (node indices)
 output:
    elem: element list after removing the duplicated elements

==== function [newnode,newelem]=removedupnodes(node,elem) ====
 [newnode,newelem]=removedupnodes(node,elem)
 removing the duplicated nodes from a mesh
 input:
   elem: integer array with dimensions of NE x 4, each row contains
         the indices of all the nodes for each tetrahedron
   node: node coordinates, 3 columns for x, y and z respectively
 output:
   newnode: nodes without duplicates
   newelem: elements with only the unique nodes

==== function [no,el]=removeisolatednode(node,elem) ====
 [no,el]=removeisolatednode(node,elem)
 remove isolated nodes: nodes that are not included in any element
 input:
     node: list of node coordinates
     elem: list of elements of the mesh
 output:
     no: node coordinates after removing the isolated nodes
     el: element list of the resulting mesh

==== function fnew=removeisolatedsurf(v,f,maxdiameter) ====
 fnew=removeisolatedsurf(v,f,maxdiameter)
 remove disjointed surface fragment filtered by using mesh diameter
 input:
    v: list of nodes of the input surface
    f: list of triangles of the input surface
    maxdiameter: maximum bounding box size for surface removal
 ouput:
    fnew: new face list after removing the components smaller than 
          maxdiameter

==== function f=surfaceclean(f,v) ====
 f=surfaceclean(f,v)
 remove surface patches that are located inside 
               the bounding box faces
 input: 
      v: surface node list, dimension (nn,3)
      f: surface face element list, dimension (be,3)  
 output:
      f: faces free of those on the bounding box

==== function eid=getintersecttri(tmppath) ====
 eid=getintersecttri(tmppath)
 get the IDs of self-intersecting elements from tetgen
 call this when tetgen complains about self-intersection
 input: 
   tmppath: working dir, use mwpath('') in most cases
 output:
   eid: an array of all intersecting surface elements, 
     one can read the corresponding node/elem by
     [no,el]=readoff(mwpath('post_vmesh.off'));

==== function elem=delendelem(elem,mask) ====
 elem=delendelem(elem,mask)
 delete elements whose nodes are all edge nodes
 input/output: 
      elem: input/output, surface/volumetric element list
      mask: of length of node number, =0 for internal nodes, =1 for edge nodes

==== function [newnode,newface]=surfreorient(node,face) ====
 [newnode,newface]=surfreorient(node,elem)
 reorder nodes in a single closed surface to ensure the norms of all
 triangles are pointing outward
 input:
    node: list of nodes
    face: list of surface triangles (each row are indices of nodes of each triangle)
 output:
    newnode: the output node list, in most cases it equals node
    newface: the face list with consistent ordering
=== # Mesh resampling and optimization ===

==== function [node,elem]=meshresample(v,f,keepratio) ====
 [node,elem]=meshresample(v,f,keepratio)
 resample mesh using CGAL mesh simplification utility
 input:
    v: list of nodes
    f: list of surface elements (each row for each triangle)
    keepratio: decimation rate, a number less than 1, as the percentage
               of the elements after the sampling
 output:
    node: the node coordinates of the sampled surface mesh
    elem: the element list of the sampled surface mesh

==== function [newno,newfc]=remeshsurf(node,face,opt) ====
 [newno,newfc]=remeshsurf(node,face,opt)
 remesh a triangular surface and the output is guaranteed to be
 free of self-intersecting element. This function is similar to 
 meshresample, but it can both downsample or upsample a mesh
 input:
	 node: list of nodes on the input suface mesh, 3 columns for x,y,z
	 face: list of trianglular elements on the surface, [n1,n2,n3,region_id]
	 opt: function parameters
	   opt.gridsize:  resolution for the voxelization of the mesh
	   opt.closesize: if there are openings, set the closing diameter
	   opt.elemsize:  the size of the element of the output surface
	   if opt is a scalar, it defines the elemsize and gridsize=opt/4
 output:
	 newno:  list of nodes on the resulting suface mesh, 3 columns for x,y,z
	 newfc:  list of trianglular elements on the surface, [n1,n2,n3,region_id]

==== function p=smoothsurf(node,mask,conn,iter,useralpha,usermethod,userbeta) ====
 p=smoothsurf(node,mask,conn,iter,useralpha,usermethod,userbeta)
 smoothing a surface mesh
 input:
    node:  node coordinates of a surface mesh
    mask:  flag whether a node is movable: 0 movable, 1 non-movable
           if mask=[], it assumes all nodes are movable
    conn:  input, a cell structure of length size(node), conn{n}
           contains a list of all neighboring node ID for node n,
           this can be computed from meshconn function
    iter:  smoothing iteration number
    useralpha: scaler, smoothing parameter, v(k+1)=(1-alpha)*v(k)+alpha*mean(neighbors)
    usermethod: smoothing method, including 'laplacian','laplacianhc' and 'lowpass'
    userbeta: scaler, smoothing parameter, for 'laplacianhc'
 output:
    p: output, the smoothed node coordinates
 recommendations
    Based on [Bade2006], 'Lowpass' method outperforms 'Laplacian-HC' in volume
    preserving and both are significantly better than the standard Laplacian method
    [Bade2006]  R. Bade, H. Haase, B. Preim, "Comparison of Fundamental Mesh 
                Smoothing Algorithms for Medical Surface Models," 
                Simulation and Visualization, pp. 289-304, 2006. 

==== function [no,el,fc,nodemap]=sortmesh(origin,node,elem,ecol,face,fcol) ====
 [no,el,fc]=sortmesh(origin,node,elem,face)
 sort nodes and elements in a mesh so that the indexed
 nodes and elements are closer to each order
 (this may reduce cache-miss in a calculation)
 input:
    origin: sorting all nodes and elements with the distance and
            angles wrt this location, if origin=[], it will be 
            node(1,:)
    node: list of nodes
    elem: list of elements (each row are indices of nodes of each element)
    ecol: list of columns in elem to participate sorting
    face: list of surface triangles (this can be omitted)
    fcol: list of columns in face to participate sorting
 output:
    no: node coordinates in the sorted order
    el: the element list in the sorted order
    fc: the surface triangle list in the sorted order (can be ignored)
    nodemap: the new node mapping order, no=node(nodemap,:)

==== function [newnode,newelem]=mergemesh(node,elem,varargin) ====
 [newnode,newelem]=mergemesh(node,elem,varargin)
 concatenate two or more tetrahedral meshes or triangular surfaces
 input: 
      node: node coordinates, dimension (nn,3)
      elem: tetrahedral element or triangle surface (nn,3) to (nn,5)
 output:
      newnode: the node coordinates after merging, dimension (nn,3)
      newelem: tetrahedral element or surfaces after merging (nn,4) or (nhn,5)
 note: you can call meshcheckrepair for the output newnode and
 newelem to remove the duplicated nodes or elements. mergemesh does
 detect self-intersecting elements when merging; to remove self-intersecting
 elements, you need to use mergesurf().
 example:
   [node1,face1,elem1]=meshabox([0 0 0],[10 10 10],1,1);
   [node2,face2,elem2]=meshasphere([5 5 13.1],3,0.3,3);
   [newnode,newelem]=mergemesh(node1,elem1,node2,elem2);
   plotmesh(newnode,newelem);
   figure;
   [newnode,newface]=mergemesh(node1,face1,node2,face2);
   plotmesh(newnode,newface,'x>5');

==== function [newnode,newelem,newface]=meshrefine(node,elem,varargin) ====
 [newnode,newelem,newface]=meshrefine(node,elem,face,opt)
 refine a tetrahedral mesh by adding new nodes or constraints
 input parameters:
      node: existing tetrahedral mesh node list
      elem: existing tetrahedral element list
      face: (optional) existing tetrahedral mesh surface triangle list
      opt:  options for mesh refinement:
        if opt is a Nx3 array, opt is treated as a list of new nodes to
          be inserted into the mesh (must be located on the surface or inside)
        if opt is a struct, it can have the following fields:
          opt.newnode: same as setting opt to an Nx3 array
          opt.reratio: radius-edge ratio, by default, iso2mesh uses 1.414
          opt.maxvol: maximum element volume
 outputs:
      newnode: node coordinates of the tetrahedral mesh
      newelem: element list of the tetrahedral mesh
      newface: mesh surface element list of the tetrahedral mesh 
             the last column denotes the boundary ID

==== function [newnode,newelem]=mergesurf(node,elem,varargin) ====
 [newnode,newelem]=mergesurf(node1,elem1,node2,elem2,...)
 merge two or more triangular meshes and split intersecting elements
 input:
      node: node coordinates, dimension (nn,3)
      elem: tetrahedral element or triangle surface (nn,3)
 output:
      newnode: the node coordinates after merging, dimension (nn,3)
      newelem: tetrahedral element or surfaces after merging (nn,4) or (nhn,5)
 note: you can call meshcheckrepair for the output newnode and
 newelem to remove the duplicated nodes or elements
 example:
   [node1,face1,elem1]=meshabox([0 0 0],[10 10 10],1,1);
   [node2,face2,elem2]=meshasphere([5 5 10],3,0.3,3);
   [newnode,newface]=mergemesh(node1,face1,node2,face2);
   plotmesh(newnode,newface,'x>5');

==== function [newnode,newelem,newelem0]=surfboolean(node,elem,varargin) ====
 [newnode,newelem,newelem0]=surfboolean(node1,elem1,op2,node2,elem2,op3,node3,elem3,...)
 merge two or more triangular meshes and resolve intersecting elements
 input:
      node: node coordinates, dimension (nn,3)
      elem: tetrahedral element or triangle surface (ne,3)
      op:  a string of a boolean operator, possible op values include
           'union' or 'or': the outter surface of the union of the enclosed space
           'inter' or 'and': the surface of the domain contained by both meshes
           'diff' or '-': the surface of the domain in mesh 1 excluding that of
                   mesh 2
           'all' or 'xor' or '+': the output contains 4 subsurfaces, identified by the 4th
                  column of newelem:
                    1: mesh 1 outside of mesh 2
                    2: mesh 2 outside of mesh 1
                    3: mesh 1 inside of mesh 2
                    4: mesh 2 inside of mesh 1
                  you can use newelem(find(mod(newelem(:,4),2)==1),:) to
                  get mesh 1 cut by mesh 2, or newelem(find(mod(newelem(:,4),2)==0),:) 
                  to get mesh 2 cut by mesh 1;
           'first': combine 1 and 3 from the output of 'all'
           'second': combine 2 and 4 from the output of 'all'
           'self': test for self-intersections; only the first mesh is
                   tested; other inputs are ignored.
           'decouple': separate two shells and make sure there is no intersection;
                   the input surfaces must be closed and ordered from outer to inner
 output:
      newnode: the node coordinates after boolean operations, dimension (nn,3)
      newelem: tetrahedral element or surfaces after boolean operations (nn,4) or (nhn,5)
      newelem0: when the operator is 'self', return the intersecting
               element list in terms of the input node list (experimental)
 example:
   [node1,face1,elem1]=meshabox([0 0 0],[10 10 10],1,1);
   [node2,face2,elem2]=meshabox([0 0 0]+5,[10 10 10]+5,1,1);
   [newnode,newface]=surfboolean(node1,face1,'union',node2,face2);
   plotmesh(newnode,newface);
   figure;
   [newnode,newface]=surfboolean(node1,face1,'diff',node2,face2);
   plotmesh(newnode,newface,'x>5');
=== # File I/O ===

==== function saveasc(v,f,fname) ====
 saveasc(v,f,fname)
 save a surface mesh to FreeSurfer ASC mesh format
 input:
      v: input, surface node list, dimension (nn,3)
      f: input, surface face element list, dimension (be,3)
      fname: output file name

==== function savedxf(node,face,elem,fname) ====
 savedxf(node,face,elem,fname)
 save a surface mesh to DXF format
 input:
      node: input, surface node list, dimension (nn,3)
      face: input, surface face element list, dimension (be,3)
      elem: input, tetrahedral element list, dimension (ne,4)
      fname: output file name

==== function savestl(node,elem,fname,solidname) ====
 savestl(node,elem,fname,solidname)
 save a tetrahedral mesh to an STL (Standard Tessellation Language) file
 input:
      node: input, surface node list, dimension Nx3
      elem: input, tetrahedral element list; if size(elem,2)==3, it is a surface
      fname: output file name
      solidname: an optional string for the name of the object

==== function savebinstl(node,elem,fname,solidname) ====
 savebinstl(node,elem,fname,solidname)
 save a tetrahedral mesh to a binary STL (Standard Tessellation Language) file
 input:
      node: input, surface node list, dimension Nx3
      elem: input, tetrahedral element list; if size(elem,2)==3, it is a surface
      fname: output file name
      solidname: an optional string for the name of the object

==== function saveinr(vol,fname) ====
 saveinr(vol,fname)
 save a surface mesh to INR Format
 input:
      vol: input, a binary volume
      fname: output file name

==== function saveoff(v,f,fname) ====
 saveoff(v,f,fname)
 save a surface mesh to Geomview Object File Format (OFF)
 input:
      v: input, surface node list, dimension (nn,3)
      f: input, surface face element list, dimension (be,3)
      fname: output file name

==== function savesmf(v,f,fname) ====
 savesmf(v,f,fname)
 save a surface mesh to smf format
 input:
      v: input, surface node list, dimension (nn,3)
      f: input, surface face element list, dimension (be,3)
      fname: output file name

==== function savesurfpoly(v,f,holelist,regionlist,p0,p1,fname,forcebox) ====
 savesurfpoly(v,f,holelist,regionlist,p0,p1,fname)
 save a set of surfaces into poly format (for tetgen)
 input:
      v: input, surface node list, dimension (nn,3)
         if v has 4 columns, the last column specifies mesh density near each node
      f: input, surface face element list, dimension (be,3)
      holelist: list of holes, each hole is represented by an internal point
      regionlist: list of regions, similar to holelist
      p0: coordinate of one of the end of the bounding box
      p1: coordinate for the other end of the bounding box
      fname: output file name
      forcebox: non-empty: add bounding box, []: automatic
                if forcebox is a 8x1 vector, it will be used to 
                specify max-edge size near the bounding box corners

==== function nedge=savegts(v,f,fname,edges) ====
 nedge=savegts(v,f,fname,edges)
 save a surface mesh to GNU Triangulated Surface Format (GTS)
 input:
      v: input, surface node list, dimension (nn,3)
      f: input, surface face element list, dimension (be,3)
      fname: output file name
      edges: edge list, if ignored, savegts will compute
 output:
      nedge: the number of unique edges in the mesh

==== function [node,elem,edges,edgemap]=readgts(fname) ====
 [node,elem,edges,edgemap]=readgts(fname)
 read GNU Triangulated Surface files (GTS)
 input:
    fname: name of the OFF data file
 output:
    node: node coordinates of the mesh
    elem: list of elements of the surface mesh
    edges: the edge list section in the GTS file (optional)
    edgemap: the face section (in terms of edge indices) in the GTS file
             (optional)

==== function savevrml(node,face,elem,fname) ====
 savevrml(node,face,elem,fname)
 save a surface mesh to VRML 1.0 format
 input:
      node: input, surface node list, dimension (nn,3)
      face: input, surface face element list, dimension (be,3)
      elem: input, tetrahedral element list, dimension (ne,4)
      fname: output file name

==== function [node,elem]=readasc(fname) ====
 [node,elem]=readasc(fname)
 read FreeSurfer ASC mesh format
 input:
      fname: name of the asc file
 output:
      node: node positions of the mesh
      elem: element list of the mesh

==== function dat=readinr(fname) ====
 vol=readinr(fname)
 load a volume from an INR file
 input:
      fname: input file name
 output:
      dat: output, data read from the inr file

==== function [node,elem,face]=readmedit(filename) ====
 [node,elem,face]=readmedit(filename)
 read Medit mesh format
 input:
    fname: name of the medit data file
 output:
    node: node coordinates of the mesh
    elem: list of elements of the mesh	    
    face: list of surface triangles of the mesh	    

==== function [node,elem]=readoff(fname) ====
 [node,elem]=readoff(fname)
 read Geomview Object File Format (OFF)
 input:
    fname: name of the OFF data file
 output:
    node: node coordinates of the mesh
    elem: list of elements of the mesh	    

==== function [node,elem]=readsmf(fname) ====
 [node,elem]=readsmf(fname)
 read simple model format (SMF)
 input: 
    fname: name of the	SMF data file
 output:
    node: node coordinates of the mesh
    elem: list of elements of the mesh

==== function [node,elem,face]=readtetgen(fstub) ====
 [node,elem,face]=readtetgen(fstub)
 read tetgen output files
 input:
    fstub: file name stub
 output:
    node: node coordinates of the tetgen mesh
    elem: tetrahedra element list of the tetgen mesh
    face: surface triangles of the tetgen mesh

==== function flag=deletemeshfile(fname) ====
 flag=deletemeshfile(fname)
 delete a given work mesh file under the working directory
 input: 
     fname: specified file name (without path)
 output:
     flag: not used

==== function binname=mcpath(fname) ====
 binname=mcpath(fname)
 get full executable path by prepending a command directory path
 parameters:
 input:
    fname: input, a file name string
 output:
    binname: output, full file name located in the bin directory
    if global variable ISO2MESH_BIN is set in 'base', it will
    use [ISO2MESH_BIN filesep cmdname] as the command full path,
    otherwise, let matlab pass the cmdname to the shell, which
    will search command in the directories listed in system
    $PATH variable.

==== function tempname=mwpath(fname) ====
 tempname=meshtemppath(fname)
 get full temp-file name by prepend working-directory and current session name
 input:
    fname: input, a file name string
 output:
    tempname: output, full file name located in the working directory
    if global variable ISO2MESH_TEMP is set in 'base', it will use it
    as the working directory; otherwise, will use matlab function tempdir
    to return a working directory.
    if global variable ISO2MESH_SESSION is set in 'base', it will be
    prepended for each file name, otherwise, use supplied file name.

==== function savemedit(node,face,elem,fname) ====
 savemedit(node,face,elem,fname)
 save a surface or tetrahedral mesh to Medit format
 input:
      node: input, surface node list, dimension (nn,3 or 4)
      face: input, surface face element list, dimension (be,3 or 4)
      elem: input, tetrahedral element list, dimension (ne,4 or 5)
      fname: output file name

==== function json=savejson(rootname,obj,varargin) ====
 json=savejson(rootname,obj,filename)
    or
 json=savejson(rootname,obj,opt)
 json=savejson(rootname,obj,'param1',value1,'param2',value2,...)
 convert a MATLAB object (cell, struct or array) into a JSON (JavaScript
 Object Notation) string
            created on 2011/09/09
 $Id: savejson.m 394 2012-12-18 17:58:11Z fangq $
 input:
      rootname: name of the root-object, if set to '', will use variable name
      obj: a MATLAB object (array, cell, cell array, struct, struct array)
      filename: a string for the file name to save the output JSON data
      opt: a struct for additional options, use [] if all use default
        opt can have the following fields (first in [.|.] is the default)
        opt.FileName [''|string]: a file name to save the output JSON data
        opt.FloatFormat ['%.10g'|string]: format to show each numeric element
                         of a 1D/2D array;
        opt.ArrayIndent [1|0]: if 1, output explicit data array with
                         precedent indentation; if 0, no indentation
        opt.ArrayToStruct[0|1]: when set to 0, savejson outputs 1D/2D
                         array in JSON array format; if sets to 1, an
                         array will be shown as a struct with fields
                         "_ArrayType_", "_ArraySize_" and "_ArrayData_"; for
                         sparse arrays, the non-zero elements will be
                         saved to _ArrayData_ field in triplet-format i.e.
                         (ix,iy,val) and "_ArrayIsSparse_" will be added
                         with a value of 1; for a complex array, the 
                         _ArrayData_ array will include two columns 
                         (4 for sparse) to record the real and imaginary 
                         parts, and also "_ArrayIsComplex_":1 is added. 
        opt.ParseLogical [0|1]: if this is set to 1, logical array elem
                         will use true/false rather than 1/0.
        opt.NoRowBracket [1|0]: if this is set to 1, arrays with a single
                         numerical element will be shown without a square
                         bracket, unless it is the root object; if 0, square
                         brackets are forced for any numerical arrays.
        opt.ForceRootName [0|1]: when set to 1 and rootname is empty, savejson
                         will use the name of the passed obj variable as the 
                         root object name; if obj is an expression and 
                         does not have a name, 'root' will be used; if this 
                         is set to 0 and rootname is empty, the root level 
                         will be merged down to the lower level.
        opt.Inf ['"$1_Inf_"'|string]: a customized regular expression pattern
                         to represent +/-Inf. The matched pattern is '([-+]*)Inf'
                         and $1 represents the sign. For those who want to use
                         1e999 to represent Inf, they can set opt.Inf to '$11e999'
        opt.NaN ['"_NaN_"'|string]: a customized regular expression pattern
                         to represent NaN
        opt.JSONP [''|string]: to generate a JSONP output (JSON with padding),
                         for example, if opt.JSON='foo', the JSON data is
                         wrapped inside a function call as 'foo(...);'
        opt.UnpackHex [1|0]: conver the 0x[hex code] output by loadjson 
                         back to the string form
        opt can be replaced by a list of ('param',value) pairs. The param 
        string is equivallent to a field in opt.
 output:
      json: a string in the JSON format (see http://json.org)
 examples:
      a=struct('node',[1  9  10; 2 1 1.2], 'elem',[9 1;1 2;2 3],...
           'face',[9 01 2; 1 2 3; NaN,Inf,-Inf], 'author','FangQ');
      savejson('mesh',a)
      savejson('',a,'ArrayIndent',0,'FloatFormat','\t%.5g')
 license:
     BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details

==== function data = loadjson(fname,varargin) ====
 data=loadjson(fname,opt)
    or
 data=loadjson(fname,'param1',value1,'param2',value2,...)
 parse a JSON (JavaScript Object Notation) file or string
         Nedialko Krouchev: http://www.mathworks.com/matlabcentral/fileexchange/25713
         François Glineur: http://www.mathworks.com/matlabcentral/fileexchange/23393
         Joel Feenstra:
         http://www.mathworks.com/matlabcentral/fileexchange/20565
 $Id: loadjson.m 394 2012-12-18 17:58:11Z fangq $
 input:
      fname: input file name, if fname contains "{}" or "[]", fname
             will be interpreted as a JSON string
      opt: a struct to store parsing options, opt can be replaced by 
           a list of ('param',value) pairs. The param string is equivallent
           to a field in opt.
 output:
      dat: a cell array, where {...} blocks are converted into cell arrays,
           and [...] are converted to arrays
 license:
     BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details 

==== function savejmesh(node,face,elem,fname,varargin) ====
 savejmesh(node,face,elem,fname,opt)
 export a mesh to the JMesh format
 input:
      node: input, node list, dimension (nn,3)
      face: input, optional, surface face element list, dimension (be,3)
      elem: input, tetrahedral element list, dimension (ne,4)
      fname: output file name
      opt: additional parameters in the form of 'parameter',value pairs
           valid parameters include:
           'Dimension': 0 - a user defined mesh, 2- a 2D mesh, 3- a 3D mesh
           'Author': a string to set the author of the mesh
           'MeshTitle': a string to set the title of the mesh
           'MeshTag': a value as the tag of the mesh data
           'MeshGroup': a value to group this mesh with other mesh data
           'Comment': a string as the additional note for the mesh data

==== function savemphtxt(node, face, elem, filename)
 ====
 savemphtxt(node, face, elem, filename)
 save tetrahedron mesh to comsol file (.mphtxt)
 author: Donghyeon Kim (danielkim<at> gist.ac.kr)
 input:
      node: input, node list, dimension (nn,3)
      face: input, surface face element list with label, dimension (be,4)
      elem: input, tetrahedron element list with label, dimension (ne,5)
      filename: input, output file name

==== function savetetgenele(elem,fname) ====
 savetetgenele(elem,fname)
 save a mesh tetrahedral element list to tetgen .ele format
 input:
      elem: tetrahedral element list, dimension (ne,4)
            columns beyound the 4rd column are treated as 
            markers and attributes associated with the element
      fname: output file name

==== function savetetgennode(node,fname) ====
 savetetgennode(node,fname)
 save a mesh node list to tetgen .node format
 input:
      node: node coordinates, dimension (nn,3)
            columns beyound the 3rd column are treated as 
            markers and attributes associated with the node
      fname: output file name

==== function saveabaqus(node,face,elem,fname,heading) ====
 saveabaqus(node,fname)
 saveabaqus(node,face,fname)
 saveabaqus(node,face,elem,fname)
 save a tetrahedral and/or surface mesh as an ABAQUS input file
 input:
      node: input, surface node list, dimension (nn,3)
      face: input, surface face element list, dimension (be,3)
      elem: input, tetrahedral element list, dimension (ne,4)
      fname: output file name
      heading: optional, a descriptive string for the mesh
=== # Volumetric image pre-processing ===

==== function islands=bwislands(img) ====
 islands=bwislands(img)
 return the indices of non-zero elements in a 2D or 3D image
 grouped by connected regions in a cell array
 input:
	 img: a 2D or 3D array
 output:
	 islands: a cell array, each cell records the indices 
		  of the non-zero elements in img for a connected 
		  region (or an island)

==== function resimg=fillholes3d(img,maxgap) ====
 resimg=fillholes3d(img,maxgap)
 close a 3D image with the speicified gap size and then fill the holes
 input:
    img: a 3D binary image
    maxgap: maximum gap size for image closing
 output:
    resimg: the image free of holes
 this function requires the image processing toolbox for matlab/octave

==== function cleanimg=deislands2d(img,sizelim) ====
 cleanimg=deislands2d(img,sizelim)
 remove isolated islands on a 2D image below speicified size limit
 input:
 	img: a 2D binary image
 	sizelim: a integer as the maximum pixel size of a isolated region
 output:
 	cleanimg: a binary image after removing islands below sizelim

==== function cleanimg=deislands3d(img,sizelim) ====
 cleanimg=deislands3d(img,sizelim)
 remove isolated islands for 3D image (for each slice)
 input:
      img: a 3D volumetric image
      sizelim: maximum island size (in pixels) for each x/y/z slice
 output:
      cleanimg: 3D image after removing the islands

==== function imgdiff=imedge3d(binimg,isdiff) ====
 imgdiff=imedge3d(binimg,isdiff)
 Extract the boundary voxels from a binary image
 author: Aslak Grinsted <ag at glaciology.net>
 input: 
   binimg: a 3D binary image
   isdiff: if isdiff=1, output will be all voxels which 
         is different from its neighbors; if isdiff=0 or 
         ignored, output will be the edge voxels of the 
         non-zero regions in binimg
 output:
   imgdiff: a 3D logical array with the same size as binimg
            with 1 for voxels on the boundary and 0 otherwise 

==== function p=internalpoint(v,aloop) ====
 p=internalpoint(v,aloop)
 imperical function to find an internal point
 of a planar polygon
 input:   
    v:     x,y,z coordinates of each node of the mesh
    aloop:  input, a single vector separated by NaN, each segment
             is a close-polygon consisted by node IDs 
 output:
    p:   output, [x y z] of an internal point of aloop

==== function vol=smoothbinvol(vol,layer) ====
 vol=smoothbinvol(vol,layer)
 perform a memory-limited 3D image smoothing
 input:
     vol: a 3D volumetric image to be smoothed
     layer: number of iterations for the smoothing
 output:
     vol: the volumetric image after smoothing

==== function vol=thickenbinvol(vol,layer) ====
 vol=thickenbinvol(vol,layer)
 thickening a binary volume by a given pixel width
 this is similar to bwmorph(vol,'thicken',3) except 
 this does it in 3d and only does thickening for 
 non-zero elements (and hopefully faster)
 input:
     vol: a volumetric binary image
     layer: number of iterations for the thickenining
 output:
     vol: the volume image after the thickening

==== function vol=thinbinvol(vol,layer) ====
 vol=thinbinvol(vol,layer)
 thinning a binary volume by a given pixel width
 this is similar to bwmorph(vol,'thin',n) except 
 this does it in 3d and only does thinning for 
 non-zero elements (and hopefully faster)
 input:
     vol: a volumetric binary image
     layer: number of iterations for the thickenining
 output:
     vol: the volume image after the thickening
=== # Mesh plotting ===

==== function hm=plotmesh(node,varargin) ====
 hm=plotmesh(node,face,elem,opt)
 plot surface and volumetric meshes
 input: 
      node: a node coordinate list, 3 columns for x/y/z; if node has a 
            4th column, it will be used to set the color at each node.
      face: a triangular surface face list; if face has a 4th column,
            it will be used to separate the surface into 
            sub-surfaces and display them in different colors;
            face can be a cell array, each element of the array represents
            a polyhedral facet of the mesh, if an element is an array with
            two array subelements, the first one is the node index, the
            second one is a scalar as the group id of the facet.
      elem: a tetrahedral element list; if elem has a 5th column,
            it will be used to separate the mesh into 
            sub-domains and display them in different colors.
      opt:  additional options for the plotting
            for simple point plotting, opt can be markers
            or color options, such as 'r.', or opt can be 
            a logic statement to select a subset of the mesh,
            such as 'x>0 & y+z<1', or an equation defining
            a plane at which a mesh cross-section is plotted, for
            example 'y=2*x'; opt can have more than one
            items to combine these options, for example: 
            plotmesh(...,'x>0','r.'); the range selector must
            appear before the color/marker specifier
 in the event where all of the above inputs have extra settings related to 
 the color of the plot, the priorities are given in the following order:
          opt > node(:,4) > elem(:,5) > face(:,4)
 output:
   hm: handle or handles (vector) to the plotted surfaces
 example:
   h=plotmesh(node,'r.');
   h=plotmesh(node,'x<20','r.');
   h=plotmesh(node,face);
   h=plotmesh(node,face,'y>10');
   h=plotmesh(node,face,'facecolor','r');
   h=plotmesh(node,elem,'x<20');
   h=plotmesh(node,elem,'x<20 & y>0');
   h=plotmesh(node,face,elem);
   h=plotmesh(node,face,elem,'linestyle','--');
   h=plotmesh(node,elem,'z=20');

==== function hm=plotsurf(node,face,varargin) ====
 hm=plotsurf(node,face,opt)
 plot 3D surface meshes
 input: 
      node: node coordinates, dimension (nn,3); if node has a 
            4th column, it will be used to set the color at each node.
      face: triangular surface face list; if face has a 4th column,
            it will be used to separate the surface into 
            sub-surfaces and display them in different colors;
            face can be a cell array, each element of the array represents
            a polyhedral facet of the mesh, if an element is an array with
            two array subelements, the first one is the node index, the
            second one is a scalar as the group id of the facet.
      opt:  additional options for the plotting, see plotmesh
 output:
   hm: handle or handles (vector) to the plotted surfaces
 example:
   h=plotsurf(node,face);
   h=plotsurf(node,face,'facecolor','r');

==== function hm=plottetra(node,elem,varargin) ====
 hm=plottetra(node,elem,opt)
 plot 3D surface meshes
 input: 
      node: a node coordinate list, 3 columns for x/y/z; if node has a 
            4th column, it will be used to set the color at each node.
      elem: a tetrahedral element list; if elem has a 5th column,
            it will be used to separate the mesh into 
            sub-domains and display them in different colors.
      opt:  additional options for a patch object, see plotmesh
 output:
   hm: handle or handles (vector) to the plotted surfaces
 example:
   h=plottetra(node,elem);
   h=plottetra(node,elem,'facealpha',0.5);

==== function [cutpos,cutvalue,facedata,elemid]=qmeshcut(elem,node,value,cutat,varargin) ====
 [cutpos,cutvalue,facedata,elemid]=qmeshcut(elem,node,value,cutat)
 fast tetrahedral mesh slicer
 input: 
   elem: integer array with dimensions of NE x 4, each row contains
         the indices of all the nodes for each tetrahedron
   node: node coordinates, 3 columns for x, y and z respectively
   value: a scalar array with the length of node numbers, can have 
          multiple columns 
   cutat: cutat can have different forms:
          if cutat is a 3x3 matrix, it defines a plane by 3 points: 
                 cutat=[x1 y1 z1;x2 y2 z2;x3 y3 z3]
          if cutat is a vector of 4 element, it defines a plane by
                 a*x+b*y+c*z+d=0  and cutat=[a b c d]
          if cutat is a single scalar, it defines an isosurface 
                 inside the mesh at value=cutat
          if cutat is a string, it defines an implicit surface
                 at which the cut is made. it must has form expr1=expr2
                 where expr1 expr2 are expressions made of x,y,z,v and
                 constants
          if cutat is a cell in the form of {'expression',scalar}, 
                 the expression will be evaluated at each node to 
                 produce a new value, then an isosurface is produced 
                 at the levelset where new value=scalar; the 
                 expression can contain constants and x,y,z,v
 output:
   cutpos: all the intersections of mesh edges by the cutat
           cutpos is similar to node, containing 3 columns for x/y/z
   cutvalue: interpolated values at the intersections, with row number
           being the num. of the intersections, column number being the 
           same as "value".
   facedata: define the intersection polygons in the form of patch "Faces"
   elemid: the index of the elem in which each intersection polygon locates
   without any output, qmeshcut generates a cross-section plot
 the outputs of this subroutine can be easily plotted using 
  % if value has a length of node:
     patch('Vertices',cutpos,'Faces',facedata,'FaceVertexCData',cutvalue,'FaceColor','interp');
  % if value has a length of elem:
     patch('Vertices',cutpos,'Faces',facedata,'CData',cutvalue,'FaceColor','flat');
=== # Miscellaneous functions ===

==== function valnew=surfdiffuse(node,tri,val,ddt,iter,type1,opt) ====
 valnew=surfdiffuse(node,tri,val,ddt,iter,type1,opt)
 apply a smoothing/diffusion process on a surface
 input:
     node: list of nodes of the surface mesh
     tri: triangular element list of the surface
     val: vector, scalar value for each node
     ddt: diffusion coefficient multiplied by delta t
     iter: iterations for applying the smoothing
     type1: indices of the nodes which will not be updated
     opt: method, 'grad' for gradient based, and 'simple' for simple average
 output:
     valnew: nodal value vector after the smoothing

==== function [node,elem,face]=volmap2mesh(img,ix,iy,iz,elemnum,maxvol,thickness,Amat,Bvec) ====
 [node,elem,face]=volmap2mesh(img,ix,iy,iz,thickness,elemnum,maxvol,A,B)
 convert a binary volume to tetrahedral mesh followed by an Affine transform
 input: 
        img, ix,iy,iz, elemnum and  maxvol: see vol2mesh.m
        thickness: scale z-dimension of the mesh to specified thickness, 
                   if thickness==0, scaling is bypassed
        Amat: a 3x3 transformation matrix
        Bvec: a 3x1 vector
        Amat and Bvec maps the image index space to real world coordnate system by
                   [x,y,z]_new=Amat*[x,y,z]_old+Bvec

==== function [isoctave verinfo]=isoctavemesh ====
 [isoctave verinfo]=isoctavemesh
 determine whether the code is running in octave
 output:
   isoctave: 1 if in octave, otherwise 0
   verinfo: a string, showing the version of octave (OCTAVE_VERSION)

==== function p=getvarfrom(ws,name) ====
 p=getvarfrom(ws,name)
 get variable value by name from specified work-space
 input:
    ws: name of the work-space, for example, 'base'
    name: name string of the variable
 output:
    p: the value of the specified variable, if the variable does not
       exist, return empty array

==== function [t,u,v,idx]=raytrace(p0,v0,node,face) ====
 [t,u,v,idx]=raytrace(p0,v0,node,face)
 perform a Havel-styled ray tracing for a triangular surface
 input:
   p0: starting point coordinate of the ray
   v0: directional vector of the ray
   node: a list of node coordinates (nn x 3)
   face: a surface mesh triangle list (ne x 3)
 output:
   t: signed distance from p to the intersection point for each surface
      triangle, if ray is parallel to the triangle, t is set to Inf
   u: bary-centric coordinate 1 of all intersection points
   v: bary-centric coordinate 2 of all intersection points
      the final bary-centric triplet is [u,v,1-u-v]
   idx: optional output, if requested, idx lists the IDs of the face
      elements that intersects the ray; users can manually calc idx by
      idx=find(u>=0 & v>=0 & u+v<=1.0 & ~isinf(t));
 Reference: 
  [1] J. Havel and A. Herout, "Yet faster ray-triangle intersection (using 
          SSE4)," IEEE Trans. on Visualization and Computer Graphics,
          16(3):434-438 (2010)
  [2] Q. Fang, "Comment on 'A study on tetrahedron-based inhomogeneous 
          Monte-Carlo optical simulation'," Biomed. Opt. Express, (in
          press)

==== function [a,b,c,d]=getplanefrom3pt(plane) ====
 [a,b,c,d]=getplanefrom3pt(plane)
 define a plane equation ax+by+cz+d=0 from three 3D points
 input: 
    plane: a 3x3 matrix with each row specifying a 3D point (x,y,z)
 output:
    a,b,c,d: the coefficient for plane equation ax+by+cz+d=0

==== function exesuff=getexeext() ====
 exesuff=getexeext()
 get meshing external tool extension names for the current platform
 output:
     exesuff: file extension for iso2mesh tool binaries

==== function exesuff=fallbackexeext(exesuffix, exename) ====
 exesuff=fallbackexeext(exesuffix, exename)
 get the fallback external tool extension names for the current platform
 input:
     exesuffix: the output executable suffix from getexeext
     exename: the executable name
 output:
     exesuff: file extension for iso2mesh tool binaries

==== function [major,minor,patchnum,extra]=iso2meshver ====
 [major,minor,patchnum,extra]=iso2meshver
      or
 v=iso2meshver
 get the version number of iso2mesh toolbox
 output:
    if you ask for a single output:
      v: a string denotes the current version number; the string is 
       typically in the following format: "major.minor.patch-extra"
       where major/minor/patch are typically integers, and extra can
       be an arbitrary string and is optional
    if you ask for 4 outputs:
     [major,minor,patchnum,extra] are each field of the version string

==== function [t,u,v,idx,xnode]=raysurf(p0,v0,node,face) ====
 [t,u,v,idx,xnode]=raysurf(p,v,node,face)
 perform a Havel-styled ray tracing for a triangular surface
 input:
   p0: list of starting points of the rays
   v0: directional vector of the rays, 
   node: a list of node coordinates (nn x 3)
   face: a surface mesh triangle list (ne x 3)
 output:
   t: distance from p0 to the intersection point for each surface
      triangle, if t(i)=NaN, no intersection was found for that ray
   u: bary-centric coordinate 1 of all intersection points
   v: bary-centric coordinate 2 of all intersection points
      the final bary-centric triplet is [u,v,1-u-v]
   idx: idx lists the IDs of the face elements that intersects 
      each ray
   xnode: optional output, if requested, xnode gives the intersection
      point coordinates; to compute manually, xnode=p0+repmat(t,1,3).*v0
 Reference: 
  [1] J. Havel and A. Herout, "Yet faster ray-triangle intersection (using 
          SSE4)," IEEE Trans. on Visualization and Computer Graphics,
          16(3):434-438 (2010)
  [2] Q. Fang, "Comment on 'A study on tetrahedron-based inhomogeneous 
          Monte-Carlo optical simulation'," Biomed. Opt. Express, (in
          press)

==== function val=getoptkey(key,default,varargin) ====
 val=getoptkey(key,default,opt)
    or
 val=getoptkey(key,default,'key1',val1,'key2',val2, ...)
 query the value of a key from a structure or a list of key/value pairs
 input:
   key: a string name for the target struct field name
   default: the default value of the key is not found
   opt: a struct object; the field names will be searched to match the 
        key input, opt can be a list of 'keyname'/value pairs
 output:
   val: val=opt.key if found, otherwise val=default

==== function newpt=rotatevec3d(pt,v1,u1,p0) ====
 newpt=rotatevec3d(pt,v1,u1,p0)
 rotate 3D points from one Cartesian coordinate system to another
 input: 
   pt: 3D points defined in a standard Cartesian system where a unitary
       z-vector is (0,0,1), 3 columns for x, y and z 
   v1: the unitary z-vector for the target coordinate system
   u1: the unitary z-vector for the source coordinate system, if ignored,
       u1=(0,0,1) 
   p0: offset of the new coordinate system, if ignored, p0=(0,0,0)
 output:
   newpt: the transformed 3D points

==== function [R,s]=rotmat2vec(u,v) ====
 [R,s]=rotmat2vec(u,v)
 the rotation matrix from vector u to v, satisfying R*u*s=v
 author: Bruno Luong
 URL:http://www.mathworks.com/matlabcentral/newsreader/view_original/827969
 input: 
   u: a 3D vector in the source coordinate system;
   v: a 3D vector in the target coordinate system;
 output:
   R: a rotation matrix to transform normalized u to normalized v
   s: a scaling factor, so that R*u(:)*s=v(:)

==== function opt=varargin2struct(varargin) ====
 opt=varargin2struct('param1',value1,'param2',value2,...)
   or
 opt=varargin2struct(...,optstruct,...)
 convert a series of input parameters into a structure
 input:
      'param', value: the input parameters should be pairs of a string and a value
       optstruct: if a parameter is a struct, the fields will be merged to the output struct
 output:
      opt: a struct where opt.param1=value1, opt.param2=value2 ...
 license:
     BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details 

==== function val=jsonopt(key,default,varargin) ====
 val=jsonopt(key,default,optstruct)
 setting options based on a struct. The struct can be produced
 by varargin2struct from a list of 'param','value' pairs
 $Id: loadjson.m 371 2012-06-20 12:43:06Z fangq $
 input:
      key: a string with which one look up a value from a struct
      default: if the key does not exist, return default
      optstruct: a struct where each sub-field is a key 
 output:
      val: if key exists, val=optstruct.key; otherwise val=default
 license:
     BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details

==== function s=mergestruct(s1,s2) ====
 s=mergestruct(s1,s2)
 merge two struct objects into one
 input:
      s1,s2: a struct object, s1 and s2 can not be arrays
 output:
      s: the merged struct object. fields in s1 and s2 will be combined in s.
 license:
     BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details 
