add mesh quality measure

git-svn-id: http://svn.code.sf.net/p/iso2mesh/code/trunk/iso2mesh@279 786e58fb-9377-0410-9ff7-e4ac0ac0635c

Author: fangq

README.txt

@@ -54,7 +54,6 @@ reference:
  (ISBI 2009), pp. 1142-1145, 2009
 
 
-== # List of functions ==
 === # Streamlined mesh generation - shortcuts ===
 
 ==== function [node,elem,face]=v2m(img,isovalues,opt,maxvol,method) ====
@@ -85,20 +84,21 @@ reference:
  output:
 	 img: a volumetric binary image at position of ndgrid(xi,yi,zi)
 
-==== function newnode=sms(node,face,iter,alpha) ====
- newnode=sms(node,face,iter,useralpha)
+==== function newnode=sms(node,face,iter,alpha,method) ====
+ newnode=sms(node,face,iter,useralpha,method)
  simplified version of surface mesh smoothing
  input:
     node:  node coordinates of a surface mesh
     face:  face element list of the surface mesh
     iter:  smoothing iteration number
     alpha: scaler, smoothing parameter, v(k+1)=alpha*v(k)+(1-alpha)*mean(neighbors)
+    method: same as in smoothsurf, default is 'laplacianhc'
  output:
     newnode: output, the smoothed node coordinates
 === # Streamlined mesh generation ===
 
-==== function [node,elem,face]=vol2mesh(img,ix,iy,iz,opt,maxvol,dofix,method,isovalues) ====
- [node,elem,face]=vol2mesh(img,ix,iy,iz,opt,maxvol,dofix,method,isovalues)
+==== 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
@@ -123,6 +123,8 @@ reference:
 	       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)
@@ -150,6 +152,8 @@ reference:
 	   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
@@ -224,7 +228,7 @@ reference:
 	 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.surfaceapprox: defines the maximum distance between the 
+	     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
@@ -242,7 +246,7 @@ reference:
 
 ==== function [node,elem,face]=cgals2m(v,f,opt,maxvol) ====
  [node,elem,face]=cgals2m(v,f,opt,maxvol)
- wrapper for CGAL 3D mesher (CGAL 3.5)
+ wrapper for CGAL 3D mesher (CGAL 3.5 and newer)
  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:
@@ -251,7 +255,7 @@ reference:
 	 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.surfaceapprox: defines the maximum distance between the 
+	     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
@@ -530,6 +534,65 @@ reference:
  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 distance (with signs) from p0. the
+       intersection coordinates can be expressed as p0+ts(i)*v0
+   idx: index to the face elements that intersect with the ray, order
+       match that of t
+
+==== function seeds=surfseeds(node,face) ====
+ seeds=surfseeds(node,face)
+ calculate a set of interior points with each enclosed by a closed
+ component of a surface
+ input:
+   node: a list of node coordinates (nn x 3)
+   face: a surface mesh triangle list (ne x 3)
+ output:
+   seeds: the interior points coordinates for each closed-surface
+          component
+
+==== function quality=meshquality(node,elem) ====
+ quality=meshquality(node,elem)
+ compute Joe-Liu mesh quality measure of a tetrahedral mesh
+ input:
+    node:  node coordinates of the mesh (nn x 3)
+    elem:  element table of a tetrahedral mesh (ne x 4)
+ output
+    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 edges=meshedge(elem) ====
+ edges=meshedge(elem)
+ return all edges in a surface or volumetric mesh
+ input:
+    elem:  element table of a mesh (support N-d space element)
+ 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. 
 === # Mesh processing and reparing ===
 
 ==== function [node,elem]=meshcheckrepair(node,elem,opt) ====
@@ -946,7 +1009,7 @@ reference:
 
 ==== function vol=smoothbinvol(vol,layer) ====
  vol=smoothbinvol(vol,layer)
- convolve a 3x3 gaussian kernel to a binary image multiple times
+ perform a memory-limited 3D image smoothing
  input:
      vol: a 3D volumetric image to be smoothed
      layer: number of iterations for the smoothing
@@ -1054,11 +1117,13 @@ reference:
            'FaceColor','interp');
 
 ==== function plottetview(session,method) ====
- runtetview(session,method)
+ plottetview(session,method)
  wrapper for tetview to plot the generated mesh
  input:
-	 session: a string to identify the output files for plotting
-        method:  method 
+	 session: a string to identify the output files for plotting, '' for
+	          the default session
+    method:  method can be 'cgalsurf' (default), 'simplify', 'cgalpoly'
+             'cgalmesh' and 'remesh'
 === # Miscellaneous functions ===
 
 ==== function valnew=surfdiffuse(node,tri,val,ddt,iter,type1,opt) ====
@@ -1076,15 +1141,16 @@ reference:
      valnew: nodal value vector after the smoothing
 
 ==== function [node,elem,face]=volmap2mesh(img,ix,iy,iz,elemnum,maxvol,thickness,Amat,Bvec) ====
- [node,elem,face]=vol2mesh(img,ix,iy,iz,thickness,elemnum,maxvol,A,B)
- convert a binary volume to tetrahedral mesh
+ [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
-         thickness: scale z-dimension of the mesh to specified thickness, if thickness==0, scaling is bypassed
+                   [x,y,z]_new=Amat*[x,y,z]_old+Bvec
 
 ==== function isoctave=isoctavemesh ====
  isoctave=isoctavemesh
@@ -1102,6 +1168,29 @@ reference:
     p: the value of the specified variable, if the variable does not
        exist, return empty array
 
+==== function [t,u,v]=raytrace(p,v,node,face) ====
+ [t,u,v]=raytrace(p,v,node,face)
+ perform a Havel-styled ray tracing for a triangular surface
+ input:
+   p: starting point coordinate of the ray
+   v: 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
+   u: bary-centric coordinate 1 of the intersection point
+   v: bary-centric coordinate 2 of the intersection point
+      the final bary-centric triplet is [u,v,1-u-v]
+  users can find the IDs of the elements intersecting with the ray by
+    idx=find(u>=0 & v>=0 & u+v<=1.0);
+ 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

TODO.txt

@@ -5,16 +5,13 @@
 == high priority ==
 
 * detect multi-value interfaces (i.e. voxels where more than 2 values meet) \
-and perturb the levelset surface mesh to avoid intersection [proposed on 02/11/09]
-
-* if one levelset contains multiple disjointed surfaces, separate them, \
-and do sub-region labeling for each component [proposed on 02/11/09]
+and perturb the levelset surface mesh to avoid intersection [proposed on 2009/02/11]
 
 == low priority ==
 
-* 2D image-based mesh generation support ? [proposed on 02/11/09]
+* 2D image-based mesh generation support ? [proposed on 2009/02/11]
 
-* adaptive refining [proposed on 02/11/09]
+* adaptive refining [proposed on 2009/02/11]
 
 * moving mesh
 
@@ -28,10 +25,14 @@ and do sub-region labeling for each component [proposed on 02/11/09]
 files were opened but not closed and can not be deleted \
 - done, by FangQ 2009/05/04
 
-* processing gray-scale images directly [proposed on 02/11/09] \
+* processing gray-scale images directly [proposed on 2009/02/11] \
 - done, by FangQ 2009/05/04
 
 * setting RNG seeds when calling CGAL tools to ensure repeatibility \
-of the meshing [proposed 07/17/2010] \
+of the meshing [proposed 2010/07/17] \
 - done, by FangQ 2010/11/08
 
+* if one levelset contains multiple disjointed surfaces, separate them, \
+and do sub-region labeling for each component [proposed on 2009/02/11]
+- done, by FangQ 2011/02/26
+

doc/Advanced_Features.txt

@@ -5,36 +5,36 @@ Advanced Features in iso2mesh
 Global variables
 
    ISO2MESH_TEMP
-          control the temporary file output directory
-
-   The default output directory is the system's temporary directory, i.e.
-   the output of mwpath(''). However, if one does not have permission to
-   write to this folder, or if this contains outputs from other users who
-   shares the same output folder (for example, /tmp on Linux), one can
-   avoid "write to file" error by setting his own output folder. To do
-   this, he need to define
+          controls the temporary file output directory
+
+   The default output directory is under the system's temporary directory,
+   i.e. the output of mwpath(''). However, if one does not have permission
+   to write to this folder, or if this contains outputs from other users
+   who shares the same output folder (for example, /tmp/iso2mesh-username
+   on Linux), one can avoid the permission error by setting his own output
+   folder. To do this, he need to define
    ISO2MESH_TEMP='/folder/path/you/can/write';
 
    in Matlab/Octave's "base workspace". The iso2mesh commands afterward
    will output the intermediate files under the specified temp directory.
 
    ISO2MESH_SESSION
-          control the output file name prefix
+          controls the output file name prefix
 
    If two users share a common output directory (for example /tmp for GNU
    Linux) on a single machine, some of the users may not able to write
    intermediate output files and encountered errors. If this happens, one
    can set "ISO2MESH_TEMP" to define a new output folder, alternatively,
    one can set ISO2MESH_SESSION to label all his output by a unique
    prefix. For example, if one define
-   ISO2MESH_SESSION='john_';
+   ISO2MESH_SESSION='foo_';
 
-   this will prepend 'john_' to all output file names, if this string is
+   this will prepend 'foo_' to all output file names, if this string is
    unique, users can produce all their output files without conflict with
    other users.
 
    ISO2MESH_BIN
-          specify the folder where to look for the external binaries
+          specifies the folder where to look for the external binaries
 
    There are several external tools are essential for iso2mesh. These
    pre-compiled binaries are saved under iso2mesh/bin folder. In our
@@ -54,14 +54,14 @@ Global variables
           sets the seed for the random number generators in the CGAL
           modules
 
-   Iso2mesh versions later than 1.0.0-pre can make reproducible meshes
-   across multiple runs. This was done by setting proper seeds for the
-   random number generators called inside the CGAL executables. By
-   default, iso2mesh will use 0x623F9A9E as the seed. If you prefer to set
-   your own seed, please define a global variable named ISO2MESH_RANDSEED
-   in the "base" workspace. The value of ISO2MESH_RANDSEED must be a
-   positive integer. If iso2mesh detects the definition of
-   ISO2MESH_RANDSEED, it will use it to seed the RNGs.
+   Iso2mesh versions later than 1.0.0 can make reproducible meshes across
+   multiple runs. This was done by setting proper seeds for the random
+   number generators called inside the CGAL executables. By default,
+   iso2mesh will use 0x623F9A9E as the seed. If you prefer to set your own
+   seed, please define a global variable named ISO2MESH_RANDSEED in the
+   "base" workspace. The value of ISO2MESH_RANDSEED must be a positive
+   integer. If iso2mesh detects the definition of ISO2MESH_RANDSEED, it
+   will use it to seed the RNGs.
 
 Intermediate outputs
 

doc/Get_Started.txt

@@ -18,15 +18,10 @@ Getting Started with iso2mesh
 
  load mydata.mat
  [node,elem,face]=v2m(mydata,0.5,5,100);
- % or you can use the verbose form:
- % [node,elem,face]=vol2mesh(mydata>0.5,1:size(mydata,1),1:size(mydata,2),...
-                           1:size(mydata,3), 5, 100, 1,'cgalsurf');
  plotmesh(node,face)
- % be careful, the last column of face is a label, should not be used for plotti
-ng, same for elem
 
    Needless to say, the first line loads the data to your current session.
-   The second line calls an iso2mesh function, 'v2m', shortcut for
+   The second line calls an iso2mesh function, 'v2m', the shortcut for
    vol2mesh, to construct a volumetric mesh from this data array. The
    first argument of v2m, "mydata", is the volumetric image you are about
    to mesh; the 2nd argument is the threshold value at which you define

doc/INSTALL.txt

@@ -17,8 +17,8 @@ System Requirements
    The recommended system configuration for running this toolbox includes
      * a computer running GNU/Linux, Windows, Mac OS (either 32bit/64bit)
      * standard installation of Matlab (v7 or above) or Octave (3.0 or
-       above) (for demo 3, you need to install image processing toolbox
-       for maltab/octave)
+       above) (for some demos, you need to install the image processing
+       toolbox for maltab/octave)
      * 20M disk space for the toolbox and the examples
      * a folder where you have write permission
 
@@ -54,17 +54,17 @@ For restricted users
    encounter is the "permission error" when saving the "pathtool" path
    list in order to add iso2mesh permanently. If this happens, you may
    work in a multi-user or network-based system. For Matlab users, you
-   typically need to create file named [5]startup.m under your home
-   directory (~/matlab/ for Linux/Unix), and put
+   typically need to create file named [5]startup.m under the [6]Matlab
+   startup directory (~/matlab/ for Linux/Unix), and put
    addpath('/path/to/iso2mesh/'); into this file. It will be automatically
-   executed when Matlab starts. For Octave, this file is [6].octaverc.
+   executed when Matlab starts. For Octave, this file is [7].octaverc.
 
    When using this toolbox under an extensively restricted mode, one may
-   encounter a "write to file" error, this may likely be caused by the
+   encounter a "fail to write" error, this may likely be caused by the
    default output folder is not writable from your account. If you do have
    another folder which you have permission to write, you need to
-   [7]define an variable ISO2MESH_TEMP in Matlab/Octave's "base workspace"
-   and set the value as the writable folder path, then rerun your meshing
+   [8]define an variable ISO2MESH_TEMP in Matlab/Octave's "base workspace"
+   and set the value as the writable folder path, then rer-un your meshing
    commands.
 
 References
@@ -74,5 +74,6 @@ References
    3. https://orbit.nmr.mgh.harvard.edu/plugins/scmcvs/cvsweb.php/?cvsroot=iso2mesh
    4. http://iso2mesh.sourceforge.net/cgi-bin/index.cgi?Doc/AddPath
    5. http://www.mathworks.com/access/helpdesk/help/techdoc/matlab_env/f8-4994.html#brlkmbe-1
-   6. http://en.wikibooks.org/wiki/MATLAB_Programming/Differences_between_Octave_and_MATLAB#startup.m
-   7. http://iso2mesh.sourceforge.net/cgi-bin/index.cgi?Advanced
+   6. http://www.mathworks.com/help/techdoc/matlab_env/f8-10506.html
+   7. http://en.wikibooks.org/wiki/MATLAB_Programming/Differences_between_Octave_and_MATLAB#startup.m
+   8. http://iso2mesh.sourceforge.net/cgi-bin/index.cgi?Advanced