1. Iso2mesh Function List

1. Iso2mesh Function List
1.1. Streamlined mesh generation - shortcuts
1.2. Streamlined mesh generation
1.3. iso2mesh main function backend
1.4. iso2mesh primitive meshing functions
1.5. Mesh decomposition and query
1.6. Mesh processing and reparing
1.7. Mesh resampling and optimization
1.8. File I/O
1.9. Volumetric image pre-processing
1.10. Mesh plotting
1.11. Miscellaneous functions

1.1. 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,regions,holes)

 [node,elem,face]=s2m(v,f,keepratio,maxvol,method)
 [node,elem,face]=s2m(v,f,keepratio,maxvol,'tetgen',regions,holes)
 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)
 input default values:
       method: if ignored, iso2mesh uses surf2mesh ('tetgen') to do the
               tetrahedral mesh generation
       regions,holes: if ignored, iso2mesh assumes both are empty

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 varargout=m2v(varargin)

 vol=m2v(node,face,Nxyz)
  or
 vol=m2v(node,face,xi,yi,zi)
 shortcut for mesh2vol, rasterizing a teterahedral mesh to a volume using graphics
 input/output: please see details in the help for mesh2vol

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

1.2. 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)

function [mask weight]=mesh2vol(node,elem,xi,yi,zi)

 [mask weight]=mesh2vol(node,face,Nxyz)
 [mask weight]=mesh2vol(node,face,[Nx,Ny,Nz])
 [mask weight]=mesh2vol(node,face,xi,yi,zi,hf)
   or
 newval=mesh2vol(node_val,face,...)
 fast rasterization of a 3D mesh to a volume with tetrahedron index labels
 date for initial version: Feb 10,2014
 input:
      node: node coordinates, dimension N by 2 or N by 3 array
      nodeval: a 4-column array, the first 3 columns are the node coordinates, 
            the last column denotes the values associated with each node
      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: the handle of a pre-created figure window for faster rendering
 output:
      mask: a 3D 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
      newval: when the node has 4 columns, the last column represents the
            values (color) at each node, the output newval is the rasterized
            mesh value map over the specified grid.
 note: This function only works for matlab
 example:
   [no,el]=meshgrid6(0:5,0:5,0:3);
   mask=mesh2vol(no,el(:,1:4),0:0.1:5,0:0.1:5,0:0.1:4);
   imagesc(mask(:,:,3))

1.3. 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)

function [mask, weight]=mesh2mask(node,face,xi,yi,hf)

 [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
 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);

1.4. 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 

==== function [node,elem]=meshgrid5(varargin)

 ====
 [node,elem]=meshgrid5(v1,v2,v3,...)
 mesh an ND rectangular lattice by splitting 
 each hypercube into 5 tetrahedra
 inspired by John D'Errico
 URL: http://www.mathworks.com/matlabcentral/newsreader/view_thread/107191
 input:
    v1,v2,v3,... - numeric vectors defining the lattice in
                   each dimension.
                   Each vector must be of length >= 1
 output:
    node - factorial lattice created from (v1,v2,v3,...)
           Each row of this array is one node in the lattice
    elem - integer array defining simplexes as references to
           rows of "node".
 example:
     [node,elem]=meshgrid5(0:5,0:6,0:4);
     plotmesh(node,elem);

==== function [node,elem]=meshgrid6(varargin)

 ====
 [node,elem]=meshgrid6(v1,v2,v3,...)
 mesh an ND rectangular lattice by splitting 
 each hypercube into 6 tetrahedra
 author: John D'Errico
 URL: http://www.mathworks.com/matlabcentral/newsreader/view_thread/107191
 input:
    v1,v2,v3,... - numeric vectors defining the lattice in
                   each dimension.
                   Each vector must be of length >= 1
 output:
    node - factorial lattice created from (v1,v2,v3,...)
           Each row of this array is one node in the lattice
    elem - integer array defining simplexes as references to
           rows of "node".
 example:
     [node,elem]=meshgrid6(0:5,0:6,0:4);
     plotmesh(node,elem);

function [node,face,centroids]=latticegrid(varargin)

 [node,face,centroids]=latticegrid(xrange,yrange,zrange,...)
 generate a 3D lattice
 input: 
   xrange, yrange, zrange ...: 1D vectors specifying the range of each
         dimension of the lattice
 output:
   node: the vertices of the 3D lattice
   face: the list of cell faces of the lattice, including both internal
         and external facets. By default, face is in the form of a cell
         array, with each row representing a face. One can use
         cell2mat(face) to convert it to an array
   centroids: the centroids of each lattice cell
 example:
    % generate a 3D lattice
    [node,face,c0]=latticegrid([1 2 4],1:3,1:4);
    plotmesh(node,face)
    % mesh the 3D lattice based on the face info
    [no,el]=surf2mesh(node,face,[],[],1,0.01,c0);
    figure; plotmesh(no,el)
    % mesh a 2-layer structure using a simple lattice
    [node,face,c0]=latticegrid([0 10],[0 5],[0 3.5 4]);
    c0(:,4)=[0.01;0.001];
    [no,el]=surf2mesh(node,face,[],[],1,[],c0);
    figure; plotmesh(no,el)

1.5. 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,varargin)

 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 loop or polyline segment from a collection of edges
 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 3D polyline or loop consisted of node IDs
 example:
    edges=[1 2;2 3;1 4;3 4;7 3;1 9;5 6;6 7;10 9; 8 10;1 8;9 3;11 11;11 12];
    loops=extractloops(edges)

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,maxnode)

 quality=meshquality(node,elem)
 compute the Joe-Liu mesh quality measure of an N-D mesh (N<=3)
 input:
    node:  node coordinates of the mesh (nn x 3)
    elem:  element table of an N-D mesh (ne x (N+1))
 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=surfacenorm(node,face,varargin)

 snorm=surfacenorm(node,face)
    or
 snorm=surfacenorm(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 surfacenorm 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

1.6. 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

1.7. 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. the new nodes must be located on the 
          surface or inside the original mesh. external nodes are
          discarded, unless the opt.extcmdopt is specified.
        if opt is a vector with a length that equals to that of node,
          it will be used to specify the desired edge-length at each node;
          setting a node value to 0 will by-pass the refinement at this node
        if opt is a vector with a length that equals to that of elem,
          it will be used as the desired maximum element volume of each
          tetrahedron; setting to 0 will by-pass the refinement of that element.
        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
          opt.sizefield: a vector specifying either the desired edge-length
              at each node, or the maximum volume constraint within each 
              tetrahedron, see above for details.
          opt.extcmdopt: by default, meshrefine can only insert nodes
              that are inside the original mesh. if one prefers to insert
              nodes that are outside of the original mesh, one can define
              this parameter to specify the meshing option (for tetgen)
              for the extended domain, i.e. the convex hull including 
              both the original and the external nodes. If not defined, 
              '-Y' option is used by default (prevent tetgen from
              inserting new nodes on the surface).
          opt.extlabel: when external nodes are inserted, the new elements
              will be assigned with an element label to group them
              together, by default, this label is 0, unless opt.extlabel
              is given
          opt.extcorelabel: when external nodes are inserted, par of the 
              new elements share the polyhedra between the inserted nodes,
              these special elements will be marked by opt.extcorelabel, 
              otherwise the label will be set to -1 
 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
 examples:
     [node,face,elem]=meshasphere([0 0 0],24,1,2);
     elem(:,5)=1;
     % inserting nodes that are inside the original mesh
     innernodes=double([1 1 1; 2 2 2; 3 3 3]);
     [newno,newel]=meshrefine(node,elem,innernodes);
     all(ismember(round(innernodes*1e10)*1e-10,round(newno*1e10)*1e-10,'rows'))
     plotmesh(newno,[],newel,'x>-3')
     % inserting nodes that are external to the original mesh
     extnodes=double([-5 -5 25;-5 5 25;5 5 25;5 -5 25]);
     [newno,newel]=meshrefine(node,elem,struct('newnode',extnodes,'extcmdopt','-Y'));
     figure;
     plotmesh(newno,[],newel,'x>-3')

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');

1.8. 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 500 2015-09-19 17:07:03Z fangq $
 input:
      rootname: the name of the root-object, when set to '', the root name
        is ignored, however, when opt.ForceRootName is set to 1 (see below),
        the MATLAB variable name will be used as the root 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, ignore to use default values.
        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.JSONP='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.SaveBinary [0|1]: 1 - save the JSON file in binary mode; 0 - text mode.
        opt.Compact [0|1]: 1- out compact JSON format (remove all newlines and tabs)
        opt can be replaced by a list of ('param',value) pairs. The param 
        string is equivallent to a field in opt and is case sensitive.
 output:
      json: a string in the JSON format (see http://json.org)
 examples:
      jsonmesh=struct('MeshNode',[0 0 0;1 0 0;0 1 0;1 1 0;0 0 1;1 0 1;0 1 1;1 1 1],... 
               'MeshTetra',[1 2 4 8;1 3 4 8;1 2 6 8;1 5 6 8;1 5 7 8;1 3 7 8],...
               'MeshTri',[1 2 4;1 2 6;1 3 4;1 3 7;1 5 6;1 5 7;...
                          2 8 4;2 8 6;3 8 4;3 8 7;5 8 6;5 8 7],...
               'MeshCreator','FangQ','MeshTitle','T6 Cube',...
               'SpecialData',[nan, inf, -inf]);
      savejson('jmesh',jsonmesh)
      savejson('',jsonmesh,'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
 created on 2011/09/09, including previous works from 
         Nedialko Krouchev: http://www.mathworks.com/matlabcentral/fileexchange/25713
            created on 2009/11/02
         François Glineur: http://www.mathworks.com/matlabcentral/fileexchange/23393
            created on  2009/03/22
         Joel Feenstra:
         http://www.mathworks.com/matlabcentral/fileexchange/20565
            created on 2008/07/03
 $Id: loadjson.m 500 2015-09-19 17:07:03Z 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. opt can have the following 
           fields (first in [.|.] is the default)
           opt.SimplifyCell [0|1]: if set to 1, loadjson will call cell2mat
                         for each element of the JSON data, and group 
                         arrays based on the cell2mat rules.
           opt.FastArrayParser [1|0 or integer]: if set to 1, use a
                         speed-optimized array parser when loading an 
                         array object. The fast array parser may 
                         collapse block arrays into a single large
                         array similar to rules defined in cell2mat; 0 to 
                         use a legacy parser; if set to a larger-than-1
                         value, this option will specify the minimum
                         dimension to enable the fast array parser. For
                         example, if the input is a 3D array, setting
                         FastArrayParser to 1 will return a 3D array;
                         setting to 2 will return a cell array of 2D
                         arrays; setting to 3 will return to a 2D cell
                         array of 1D vectors; setting to 4 will return a
                         3D cell array.
           opt.ShowProgress [0|1]: if set to 1, loadjson displays a progress bar.
 output:
      dat: a cell array, where {...} blocks are converted into cell arrays,
           and [...] are converted to arrays
 examples:
      dat=loadjson('{"obj":{"string":"value","array":[1,2,3]}}')
      dat=loadjson(['examples' filesep 'example1.json'])
      dat=loadjson(['examples' filesep 'example1.json'],'SimplifyCell',1)
 license:
     BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details 

function json=saveubjson(rootname,obj,varargin)

 json=saveubjson(rootname,obj,filename)
    or
 json=saveubjson(rootname,obj,opt)
 json=saveubjson(rootname,obj,'param1',value1,'param2',value2,...)
 convert a MATLAB object (cell, struct or array) into a Universal 
 Binary JSON (UBJSON) binary string
 created on 2013/08/17
 $Id: saveubjson.m 465 2015-01-25 00:46:07Z fangq $
 input:
      rootname: the name of the root-object, when set to '', the root name
        is ignored, however, when opt.ForceRootName is set to 1 (see below),
        the MATLAB variable name will be used as the root name.
      obj: a MATLAB object (array, cell, cell array, struct, struct array)
      filename: a string for the file name to save the output UBJSON data
      opt: a struct for additional options, ignore to use default values.
        opt can have the following fields (first in [.|.] is the default)
        opt.FileName [''|string]: a file name to save the output JSON data
        opt.ArrayToStruct[0|1]: when set to 0, saveubjson 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 [1|0]: 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, saveubjson
                         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.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 and is case sensitive.
 output:
      json: a binary string in the UBJSON format (see http://ubjson.org)
 examples:
      jsonmesh=struct('MeshNode',[0 0 0;1 0 0;0 1 0;1 1 0;0 0 1;1 0 1;0 1 1;1 1 1],... 
               'MeshTetra',[1 2 4 8;1 3 4 8;1 2 6 8;1 5 6 8;1 5 7 8;1 3 7 8],...
               'MeshTri',[1 2 4;1 2 6;1 3 4;1 3 7;1 5 6;1 5 7;...
                          2 8 4;2 8 6;3 8 4;3 8 7;5 8 6;5 8 7],...
               'MeshCreator','FangQ','MeshTitle','T6 Cube',...
               'SpecialData',[nan, inf, -inf]);
      saveubjson('jsonmesh',jsonmesh)
      saveubjson('jsonmesh',jsonmesh,'meshdata.ubj')
 license:
     BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details

function data = loadubjson(fname,varargin)

 data=loadubjson(fname,opt)
    or
 data=loadubjson(fname,'param1',value1,'param2',value2,...)
 parse a JSON (JavaScript Object Notation) file or string
 created on 2013/08/01
 $Id: loadubjson.m 487 2015-05-06 18:19:07Z fangq $
 input:
      fname: input file name, if fname contains "{}" or "[]", fname
             will be interpreted as a UBJSON 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. opt can have the following 
           fields (first in [.|.] is the default)
           opt.SimplifyCell [0|1]: if set to 1, loadubjson will call cell2mat
                         for each element of the JSON data, and group 
                         arrays based on the cell2mat rules.
           opt.IntEndian [B|L]: specify the endianness of the integer fields
                         in the UBJSON input data. B - Big-Endian format for 
                         integers (as required in the UBJSON specification); 
                         L - input integer fields are in Little-Endian order.
           opt.NameIsString [0|1]: for UBJSON Specification Draft 8 or 
                         earlier versions (JSONLab 1.0 final or earlier), 
                         the "name" tag is treated as a string. To load 
                         these UBJSON data, you need to manually set this 
                         flag to 1.
 output:
      dat: a cell array, where {...} blocks are converted into cell arrays,
           and [...] are converted to arrays
 examples:
      obj=struct('string','value','array',[1 2 3]);
      ubjdata=saveubjson('obj',obj);
      dat=loadubjson(ubjdata)
      dat=loadubjson(['examples' filesep 'example1.ubj'])
      dat=loadubjson(['examples' filesep 'example1.ubj'],'SimplifyCell',1)
 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

1.9. 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,nobd)

 vol=thinbinvol(vol,layer,nobd)
 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 run thinning for 
 non-zero elements (and hopefully faster)
 input:
     vol: a volumetric binary image
     layer: number of iterations for the thickenining
     nobd: (optional) if set to 1, boundaries will not 
            erode. if not given, nobd=0.
 output:
     vol: the volume image after the thinning operations

1.10. 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 (2d manifold) or polylines (1d manifold)
 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');
   h=plotsurf(node,edges,'linestyle','-','linewidth',2,'color','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 hh=plotedges(node,edges,varargin)

 hm=plotedges(node,edges,opt)
   or
 hm=plotedges(node,loops,opt)
 plot a 3D polyline or close loop (1d manifold)
 input: 
      node: node coordinates, dimension (nn,3); if node has a 
            4th column, it will be used to set the color at each node.
      edges:edge list: a 2-column index array, with each row being
            an edge connecting the two indexed node
      loops:loops is an NaN separated integer array, with each segment
            denoting a 3D polyline or loop represented by a list of node
            indices
      opt:  additional options for the plotting, see plotmesh
 output:
      hm: handle or handles (vector) to the plotted surfaces
 example:
   h=plotedges(node,[1 2 3 4 5 nan 6 7 8 9]);
   h=plotedges(node,edges,'marker','o','linewidth',2,'color','r');

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');

1.11. 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 
 -- this function is part of JSONlab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
len=length(varargin); opt=struct; if(len==0) return; end i=1; while(i<=len)
    if(isstruct(varargin{i}))
        opt=mergestruct(opt,varargin{i});
    elseif(ischar(varargin{i}) && i<len)
        opt=setfield(opt,lower(varargin{i}),varargin{i+1});
        i=i+1;
    else
        error('input must be in the form of ...,name,value,... pairs or structs');
    end
    i=i+1;
end

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
 -- this function is part of JSONlab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
val=default; if(nargin<=2) return; end opt=varargin{1}; if(isstruct(opt))
    if(isfield(opt,key))
       val=getfield(opt,key);
    elseif(isfield(opt,lower(key)))
       val=getfield(opt,lower(key));
    end
end

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 
 -- this function is part of jsonlab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
if(~isstruct(s1) || ~isstruct(s2))
    error('input parameters contain non-struct');
end if(length(s1)>1 || length(s2)>1)
    error('can not merge struct arrays');
end fn=fieldnames(s2); s=s1; for i=1:length(fn)
    s=setfield(s,fn{i},getfield(s2,fn{i}));
end

function node=orthdisk(c0,c1,r,ndiv)

 node=orthdisk(c0,c1,r,ndiv)
 Defining a 3D disk that is orthogonal to the vector c1-c0 
 input:
     c0: a 1x3 vector for the origin
     c1: a 1x3 vector to define a direction vector c1-c0
     r: the radius of the disk that is orthogonal to c1-c0, passing through c0
     ndiv: division count to approximate a circle by a polygon, if ignored, ndiv=20
 output:
     node: the 3D vertices of the disk

function newdata=struct2jdata(data,varargin)

 newdata=struct2jdata(data,opt,...)
 convert a JData object (in the form of a struct array) into an array
 input:
      data: a struct array. If data contains JData keywords in the first
            level children, these fields are parsed and regrouped into a
            data object (arrays, trees, graphs etc) based on JData 
            specification. The JData keywords are
               "_ArrayType_", "_ArraySize_", "_ArrayData_"
               "_ArrayIsSparse_", "_ArrayIsComplex_"
      opt: (optional) a list of 'Param',value pairs for additional options 
           The supported options include
               'Recursive', if set to 1, will apply the conversion to 
                            every child; 0 to disable
 output:
      newdata: the covnerted data if the input data does contain a JData 
               structure; otherwise, the same as the input.
 examples:
      obj=struct('_ArrayType_','double','_ArraySize_',[2 3],
                 '_ArrayIsSparse_',1 ,'_ArrayData_',null);
      ubjdata=struct2jdata(obj);
 license:
     BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details 
Powered by Habitat