Showing revision 1

iso2mesh: an image-based 3D surface and volumetric mesh generator

Author: Qianqian Fang <fangq at nmr.mgh.harvard.edu>

Martinos Center for Biomedical Imaging
Massachusetts General Hospital (Harvard Medical School)
Bldg. 149, 13th St., Charlestown, MA 02148
Version: 0.9.0 (Tofu alpha) License: GPL v2 or later (see COPYING)
(this license does not cover the binaries in the bin/
directory, see section III for more details)
URL: http://iso2mesh.sf.net

Table of Content

1. Introduction
2. List of functions
2.1. Streamlined Mesh Generation (wrapper)
2.2. Mesh simplification
2.3. Mesh repair
2.4. Mesh decomposition and query
2.5. Binary image pre-processing (hole-filling and island removing)
2.6. File IO
3. Acknowledgement

1. Introduction

"iso2mesh" is a Matlab/Octave-based mesh generation toolbox. It contains a number of optimized mesh processing scripts and interfaces to several free meshing tools, such as surface mesh simplifier and extraction utilities (based on CGAL), tetgen, and mesh validation&repairing utility (based on JMeshLib). Iso2mesh toolbox provides a simple yet powerful way to create quality tetrahedral volumetric mesh directly from 3D binary, segmented or gray-scale images such as MRI or CT scans.

To be able to create volumetric finite-element mesh from 3D volumetric images is of great interest for researchers in the field of image-based modeling and analysis. Unfortunately, very limited software packages for this purpose can be found in publications and market; they are either limited in functionalities or very expensive (such as Amira from Mercury Computer Systems). This toolbox was developed as a free alternative to the expensive commertial tools and provides researchers a complete workflow from 3D binary image preprocessing (hole-filling and island-removing), to surface mesh modeling (extraction, remeshing, validation, repairing, and smoothing) and volumetric mesh creation. Users can either choose the streamlined wrapper scripts to create 3D mesh with a fully automated process, or to call invidual subroutines to perform specific meshing tasks, such as closed surface extraction or boundary extraction.

The details of this toolbox can be found in the following reference:

  • Qianqian Fang and David Boas, "Tetrahedral mesh generation from volumetric binary and gray-scale images," Proceedings of IEEE International Symposium on Biomedical Imaging (ISBI 2009), pp. 1142-1145, 2009

2. List of functions

2.1. Streamlined Mesh Generation (wrapper)

vol2mesh.m

function [node,elem,bound]=vol2mesh(img,ix,iy,iz,opt,maxvol,dofix,method,isovalues) 
  or
function [node,elem,face]=v2m(img,isovalues,opt,maxvol,method) 

   convert a binary (or multi-valued) volume to tetrahedral mesh
   author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
   inputs: 
          img: a volumetric binary image 
          ix,iy,iz: subvolume selection indices in x,y,z directions
          opt: see vol2surf 
          maxvol: target maximum tetrahedral elem volume
          dofix: 1: perform mesh validation&repair, 0: skip repairing
          method: see vol2surf
          isovalues: extract level-sets at the specified isovalues and
                fill the volume in-between; if not specified, will use all
                distinct labels in the image.
   outputs:
          node: output, node coordinates of the tetrahedral mesh
          elem: output, element list of the tetrahedral mesh
          bound: output, mesh surface element list of the tetrahedral mesh
               the last column denotes the boundary ID

vol2surf.m

function [no,el,regions,holes]=vol2surf(img,ix,iy,iz,opt,dofix,method,isovalues) 
  or
function [no,el,regions,holes]=v2s(img,isovalues,opt,method) 

   converting a 3D volumetric image to surfaces
   author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
   inputs:
          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':
	      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(1,2,...).radbound: same as above, for each levelset
	    if method is 'simplify':
	      opt=a float number<1: compression rate for surf. simplification
	      opt.keeyratio=a float less than 1: same as above, same for all surf.
	      opt(1,2,..).keeyratio: 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,..).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
	  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: extract level-sets at the specified isovalues and
                fill the volume in-between; if not specified, will use all
                distinct labels in the image.
   outputs:
	  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]

surf2mesh.m

function [node,elem,bound]=surf2mesh(v,f,p0,p1,keepratio,maxvol,method) 
  or
function [node,elem,face]=s2m(v,f,keepratio,maxvol,method) 

   surf2mesh - create quality volumetric mesh from isosurface patches
   author: fangq (fangq<at> nmr.mgh.harvard.edu)
   parameters:
        v: input, isosurface node list, dimension (nn,3)
        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
        node: output, node coordinates of the tetrahedral mesh
        elem: output, element list of the tetrahedral mesh
        bound: output, mesh surface element list of the tetrahedral mesh 
               the last column denotes the boundary ID

2.2. Mesh simplification

meshresample.m

function [node,elem]=meshresample(v,f,elemnum) 
   meshresample: resample mesh using CGAL mesh simplification utility
   author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)

bbxflatsegment.m

function seg=bbxflatsegment(node,loop) 
   bbxflatsegment: decompose edge loops into flat segments
                    alone x/y/z planes of the bounding box
   author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
   parameters:   
      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 
      seg:   output, a single vector separated by NaN, each segment
               is a close-polygon on x/y/z plane

2.3. Mesh repair

meshcheckrepair.m

function [node,elem]=meshcheckrepair(node,elem,opt) 
   check and repair a surface mesh
   author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
   parameters:
        node: input/output, surface node list, dimension (nn,3)
        elem: input/output, surface face element list, dimension (be,3)
        opt: options, including
              'duplicated': remove duplicated elements
              'isolated': remove isolated nodes
              'deep': call external jmeshlib to remove non-manifold vertices

removedupelem.m

function elem=removedupelem(elem) 
   remove doubly duplicated elements

removeisolatednode.m

function [no,el]=removeisolatednode(node,elem) 
   remove isolated nodes: nodes that are not included in any element

delendelem.m

function elem=delendelem(elem,mask) 
   delendelem - delete elements whose nodes are all edge nodes
   author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
   parameters: 
        elem: input/output, surface/volumetric element list
        mask: node label, =0 for internal nodes, =1 for edge nodes

surfaceclean.m

function f=surfaceclean(f,v) 
   surfaceclean: remove surface patches that are located inside 
                 the bounding box faces
   author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
   parameters: 
        v: input, surface node list, dimension (nn,3)
        f: input, surface face element list, dimension (be,3)  
        f: output, faces free of those on the bounding box

smoothsurf.m

function nodenew=smoothsurf(node,mask,conn,iter,useralpha,usermethod,userbeta) or function nodenew=sms(node,face,iter,alpha) 

   smoothsurf: smooth a surface mesh by 3 smoothing algorithms
   author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
   input:
      node:  node coordinates of a surface mesh
      mask: of length of node number, =0 for internal nodes, =1 for edge nodes
      conn:  input, a cell structure of length size(node), conn{n}
             contains a list of all neighboring node ID for node n
      iter:  smoothing iteration number
      useralpha: scaler, smoothing parameter, v(k+1)=alpha*v(k)+(1-alpha)*mean(neighbors)
      usermethod: smoothing method, including 'laplacian','laplacianhc' and 'lowpass'
      userbeta: scaler, smoothing parameter, for 'laplacianhc'
   output:
      nodenew: output, the smoothed node coordinates

2.4. Mesh decomposition and query

finddisconnsurf.m

function facecell=finddisconnsurf(f) 
   subroutine to extract disconnected surfaces from a 
   cluster of surfaces
   author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
   parameters:
       facelist: input, node indices for all surface triangles
       facecell: separated disconnected surface node indices

maxsurf.m

function f=maxsurf(facecell) 
   return a surface with contains the most elements in a 
   collection of surfaces (stored in a cell array)

extractloops.m

function loops=extractloops(edges) 
   extractloops: extract individual loops from an edge table of a loop
                 collection
   author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
   parameters:   
      edges:  two column matrix recording the starting/ending 
               points of all edge segments
      loops:  output, a single vector separated by NaN, each segment
               is a close-polygon consisted by node IDs

flatsegment.m

function mask=flatsegment(node,edge) 
   flatsegment: decompose edge loops into flat segments
                alone arbitrary planes of the bounding box
   author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
    this code is fragile: it can not handle curves with many co-linear
    nodes near the corner point


   parameters:   
      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 
      mask:  output, a cell, each element is a close-polygon 
             on x/y/z plane

internalpoint.m

function p=internalpoint(v,aloop) 
   internalpoint: imperical function to find an internal point
                  of a planar polygon
   author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
   parameters:   
      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 
      p:   output, [x y z] of an internal point of aloop

meshconn.m

function [conn,connnum,count]=meshconn(elem,nn); 
   meshconn: create node neighbor list from a mesh
   author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
   parameters:
      elem:  element table of a mesh
      nn  :  total node number of the mesh
      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

orderloopedge.m

function newedge=orderloopedge(edge) 
   order the node list of a simple loop based on connection sequence
   author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
   inputs: 
          edge: a loop consisted by a sequence of edges, each row 
                is an edge with two integers: start/end node index
          newedge: reordered edge node list
   (this subroutine can not process bifercation)

surfedge.m

function tri=surfedge(f) 
   surfedge: find the edge of an open surface
   author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
   parameters:
        v: input, surface node list, dimension (nn,3)
        f: input, surface face element list, dimension (be,3)

2.5. Binary image pre-processing (hole-filling and island removing)

bwislands.m

function islands=bwislands(img) 
   decompose a 2D binary image into isolated regions
   and return the indices of the region pixels as a cell array
   author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
   parameters:
        img: a binary (0-1) 2D image
        islands:a cell array of length of region numbers, each 
                element stores all indices in img for that region

deislands2d.m

function cleanimg=deislands2d(img,sizelim) 
   remove all small islands smaller than specified size from  
   a binary image
   author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
   parameters:
       	img: a binary (0-1) 2D image
        sizelim: maximum pixel counts of the islands to be removed
        cleanimg: a binary image after removing the small islands

deislands3d.m

function cleanimg=deislands3d(img,sizelim) 
   remove all small islands smaller than specified size from
   a 3D binary image (it applies deislands2d for each image slice)
   author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
   parameters:
        img: a binary (0-1) 3D image
        sizelim: maximum pixel counts of the islands to	be removed
       	cleanimg: a binary image after removing	the small islands

2.6. File IO

readoff.m

function [node,elem]=readoff(fname) 
   readoff: read  Geomview Object File Format
   by FangQ, 2008/03/28

readsmf.m

function [node,elem]=readsmf(fname) 
   readsmf: read simple model format
   by FangQ, 2007/11/21

readtetgen.m

function [node,elem,bound]=readtetgen(fstub) 
   readtetgen: read tetgen output files
   author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
   parameters:
      fstub: file name stub

saveoff.m

function saveoff(v,f,fname) 
   saveoff: save a surface mesh to  Geomview Object File Format
   author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
   parameters:
        v: input, surface node list, dimension (nn,3)
        f: input, surface face element list, dimension (be,3)
        fname: output file name

savesmf.m

function savesmf(v,f,fname) 
   savesmf: save a surface mesh to smf format
   author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
   parameters:
        v: input, surface node list, dimension (nn,3)
        f: input, surface face element list, dimension (be,3)
        fname: output file name

savesurfpoly.m

function savesurfpoly(v,f,p0,p1,fname) 
   meshconn: create node neighbor list from a mesh
   author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
   parameters:
        v: input, surface node list, dimension (nn,3)
        f: input, surface 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]
        fname: output file name

mwpath.m

function fullname=mwpath(fname) 
    get full temp-file name by prepend working-directory and current session name.
    author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)


    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.

mcpath.m

function fullname=mcpath(cmdname) 
    get full executable path by prepending a command directory path


    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.

3. Acknowledgement

This toolbox interacts with a number external meshing tools to perform the essential functionalities. These tools are listed below:

bin/tetgen

  • Summary:tetgen is a compact and fast 3D mesh generator, it performs
  • License: (IMPORTANT) tetgen is free for academic research and non-commertial uses only.
  • URL: http://tetgen.berlios.de/
  • Author: Hang Si <si at wias-berlin.de>
Research Group: Numerical Mathematics and Scientific Computing
Weierstrass Institute for Applied Analysis and Stochastics
Mohrenstr. 39, 10117 Berlin, Germany

bin/cgalsurf

bin/cgalmesh and bin/cgalpoly

  • Summary: cgalmesh and cgalpoly are utilities to produce surface and volumetric meshes from a multi-valued volumetric image
  • Source: it is a slightly modified version from Mesh_3 example from CGAL 3.5
  • License: CGAL v3 core library is licensed under QPL (Q Public License) other modules are under the Lesser General Public License (LGPL)
  • URL: http://www.cgal.org/Manual/3.5/doc_html/cgal_manual/Mesh_3/Chapter_main.html

bin/cgalsimp2

bin/meshfix

  • Summary: meshfix is adapted from the sample code of JMeshLib
  • License: GPL (GNU General Public License) v2 or later
  • URL:http://jmeshlib.sourceforge.net/
  • Author:Marco Attene <attene at ge.imati.cnr.it>
Istituto di Matematica Applicata e Tecnologie Informatiche
Consiglio Nazionale delle Ricerche
Via De Marini, 6 (Torre di Francia)
16149 Genoa - ITALY

bin/tetview

Note: iso2mesh and the above meshing utilities are considered as an "aggregate" rather than "derived work", based on the definitions in GPL FAQ (http://www.gnu.org/licenses/gpl-faq.html#MereAggregation) Therefore, the license of iso2mesh and these utilities are independent. The iso2mesh license only applies to the scripts and documents/data in this package and exclude those programs stored in the bin/ directory. The source codes of the modified meshing utilities are available separatly at iso2mesh's website and retain their upstream licenses.

Your acknowledgement of iso2mesh in your publications or presentations would be greatly appreciated by the author of this toolbox. The citation information can be found in the Introduction section.

Powered by Habitat