Adds interface for the multi-zone mesh generation capabilities of TetGen (#923)

* Update tetrahedralize.h

* Update tetrahedralize.h

* WIP: adding overload for allowing holes and regions in tetgen mesh

* Update tetrahedralize.h

* Update tetrahedralize.h

* WIP: adding overload for allowing holes and regions in tetgen mesh

* WIP: skeleton code for tetrahedralize supporting holes and regions

* removes comments

* call to overloaded mesh_to_tetgenio()

* overload to allow for holes and regions

* WIP: overload to allow for holes and regions

* fixes typos

* fixed function signature

* adds comments

* adds overload for returning regionID information

* added overload for holes and regions

* added Eigen wrapper

* WIP: add support for neighborlist and point2tet

* WIP: adding point2tet and neighborlist info

* adds support for face to tet list

* WIP: bug fixes, regions, holes

* support for holes, regions, point2tet, face2tet, tet neighbors

* Update tetgen to 1.5.1.


Former-commit-id: b4eaa0dac7af02d88e358bad4033d1afbb128283

include/igl/copyleft/tetgen/mesh_to_tetgenio.cpp

@@ -13,6 +13,94 @@
 // STL includes
 #include <cassert>
 
+
+IGL_INLINE bool igl::copyleft::tetgen::mesh_to_tetgenio(
+  const std::vector<std::vector<REAL > > & V,
+  const std::vector<std::vector<int> > & F,
+  const std::vector<std::vector<REAL> > & H, 
+  const std::vector<std::vector<REAL> > & R, 
+  tetgenio & in)
+{
+  using namespace std;
+  in.firstnumber = 0;
+  in.numberofpoints = V.size();
+  in.pointlist = new REAL[in.numberofpoints * 3];
+  //loop over points
+  for(size_t i = 0; i < (size_t)V.size(); i++)
+  {
+    assert(V[i].size() == 3);
+    in.pointlist[i*3+0] = V[i][0];
+    in.pointlist[i*3+1] = V[i][1];
+    in.pointlist[i*3+2] = V[i][2];    
+  }
+  in.numberoffacets = F.size();
+  in.facetlist = new tetgenio::facet[in.numberoffacets];
+  in.facetmarkerlist = new int[in.numberoffacets];
+
+  // loop over face
+  for(size_t i = 0;i < (size_t)F.size(); i++)
+  {
+    in.facetmarkerlist[i] = i;
+    tetgenio::facet * f = &in.facetlist[i];
+    f->numberofpolygons = 1;
+    f->polygonlist = new tetgenio::polygon[f->numberofpolygons];
+    f->numberofholes = 0;
+    f->holelist = NULL;
+    tetgenio::polygon * p = &f->polygonlist[0];
+    p->numberofvertices = F[i].size();
+    p->vertexlist = new int[p->numberofvertices];
+    // loop around face
+    for(int j = 0;j < (int)F[i].size(); j++)
+    {
+      p->vertexlist[j] = F[i][j];
+    }
+  }
+  
+  in.numberofholes = H.size(); 
+  in.holelist = new double[3 * in.numberofholes];
+  // loop over holes
+  for(size_t holeID = 0, nHoles = H.size(); holeID < nHoles; holeID++)
+  {
+    in.holelist[holeID * 3 + 0] = H[holeID][0]; 
+    in.holelist[holeID * 3 + 1] = H[holeID][1];
+    in.holelist[holeID * 3 + 2] = H[holeID][2];
+  }	  
+
+  in.numberofregions = R.size();
+  in.regionlist = new REAL[ 5 * in.numberofregions];
+  // loop over regions
+  for(size_t regionID = 0, nRegions = R.size(); regionID < nRegions; regionID++)
+  {
+    in.regionlist[regionID * 5 + 0] = R[regionID][0]; 
+    in.regionlist[regionID * 5 + 1] = R[regionID][1];
+    in.regionlist[regionID * 5 + 2] = R[regionID][2];
+    in.regionlist[regionID * 5 + 3] = R[regionID][3];
+    in.regionlist[regionID * 5 + 4] = R[regionID][4];
+  }	  
+
+  return true;
+    
+}	
+
+template <typename DerivedV, typename DerivedF, typename DerivedH, typename DerivedR>
+IGL_INLINE bool igl::copyleft::tetgen::mesh_to_tetgenio(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  const Eigen::PlainObjectBase<DerivedH>& H,
+  const Eigen::PlainObjectBase<DerivedR>& R,
+  tetgenio & in)
+{
+  using namespace std;
+  vector<vector<REAL> > vV, vH, vR;
+  vector<vector<int> > vF;
+  matrix_to_list(V,vV);
+  matrix_to_list(F,vF);
+  matrix_to_list(H, vH);
+  matrix_to_list(R, vR);
+  return mesh_to_tetgenio(vV,vF,vH,vR,in);
+}
+
+
 IGL_INLINE bool igl::copyleft::tetgen::mesh_to_tetgenio(
   const std::vector<std::vector<REAL > > & V, 
   const std::vector<std::vector<int> > & F, 

include/igl/copyleft/tetgen/mesh_to_tetgenio.h

@@ -43,6 +43,38 @@ namespace igl
         const Eigen::PlainObjectBase<DerivedV>& V,
         const Eigen::PlainObjectBase<DerivedF>& F,
         tetgenio & in);
+    
+
+      // Load a vertex list and face list into a tetgenio object
+      // Inputs:
+      //   V  #V by 3 vertex position list
+      //   F  #F list of polygon face indices into V (0-indexed)
+      //   H  #H list of seed point inside each hole
+      //   R  #R list of seed point inside each region	
+      // Outputs:
+      //   in  tetgenio input object
+      // Returns true on success, false on error
+      IGL_INLINE bool mesh_to_tetgenio(
+        const std::vector<std::vector<REAL> > & V,
+	const std::vector<std::vector<int> > & F,
+	const std::vector<std::vector<REAL > > & H, 
+	const std::vector<std::vector<REAL > > & R, 
+	tetgenio & in);	
+
+
+      // Wrapper with Eigen types
+      // Templates:
+      //   DerivedV  real-value: i.e. from MatrixXd
+      //   DerivedF  integer-value: i.e. from MatrixXi
+      //   DerivedH  real-value
+      //   DerivedR  real-value		      
+      template <typename DerivedV, typename DerivedF, typename DerivedH, typename DerivedR>
+      IGL_INLINE bool mesh_to_tetgenio(
+	const Eigen::PlainObjectBase<DerivedV>& V,
+	const Eigen::PlainObjectBase<DerivedF>& F,
+	const Eigen::PlainObjectBase<DerivedH>& H, 
+	const Eigen::PlainObjectBase<DerivedR>& R, 
+	tetgenio& in);	
     }
   }
 }

include/igl/copyleft/tetgen/tetgenio_to_tetmesh.cpp

@@ -13,6 +13,120 @@
 // STL includes
 #include <iostream>
 
+IGL_INLINE bool igl::copyleft::tetgen::tetgenio_to_tetmesh(
+  const tetgenio & out,
+  std::vector<std::vector<REAL > > & V,
+  std::vector<std::vector<int> > & T,
+  std::vector<std::vector<int > > & F,
+  std::vector<std::vector<REAL > >&  R, 
+  std::vector<std::vector<int > >& N,
+  std::vector<std::vector<int > >& PT,
+  std::vector<std::vector<int > >& FT,
+  size_t & nR ) 
+{
+   using namespace std;
+   // process points
+   if(out.pointlist == NULL)
+   {
+     cerr<<"^tetgenio_to_tetmesh Error: point list is NULL\n"<<endl;
+     return false;
+   }
+   V.resize(out.numberofpoints,vector<REAL>(3));
+   // loop over points
+   for(int i = 0;i < out.numberofpoints; i++)
+   {
+     V[i][0] = out.pointlist[i*3+0];
+     V[i][1] = out.pointlist[i*3+1];
+     V[i][2] = out.pointlist[i*3+2];
+   }
+
+   // process tets
+   if(out.tetrahedronlist == NULL)
+   {
+     cerr<<"^tetgenio_to_tetmesh Error: tet list is NULL\n"<<endl;
+     return false;
+   }
+
+   // When would this not be 4?
+   assert(out.numberofcorners == 4);
+   T.resize(out.numberoftetrahedra,vector<int>(out.numberofcorners));
+   int min_index = 1e7;
+   int max_index = -1e7;
+   // loop over tetrahedra
+   for(int i = 0; i < out.numberoftetrahedra; i++)
+   {
+     for(int j = 0; j<out.numberofcorners; j++)
+     {
+       int index = out.tetrahedronlist[i * out.numberofcorners + j];
+       T[i][j] = index;
+       min_index = (min_index > index ? index : min_index);
+       max_index = (max_index < index ? index : max_index);
+     }
+   }
+  
+   assert(min_index >= 0);
+   assert(max_index >= 0);
+   assert(max_index < (int)V.size());
+ 
+   cout<<out.numberoftrifaces<<endl;
+
+   // When would this not be 4?
+   F.clear();
+   // loop over tetrahedra
+   for(int i = 0; i < out.numberoftrifaces; i++)
+   {
+     if(out.trifacemarkerlist[i]>=0)
+     {
+       vector<int> face(3);
+       for(int j = 0; j<3; j++)
+       {
+         face[j] = out.trifacelist[i * 3 + j];
+       }
+       F.push_back(face);
+     }
+   }
+
+   R.resize(out.numberoftetrahedra, vector<REAL>(1));
+   unordered_map<REAL, REAL> hashUniqueRegions;
+   for(size_t i = 0; i < out.numberoftetrahedra; i++)
+   {
+	R[i][0] = out.tetrahedronattributelist[i];
+	hashUniqueRegions[R[i][0]] = i;
+   }
+   // extract region marks
+   nR = hashUniqueRegions.size();
+
+   // extract neighbor list 
+   N.resize(out.numberoftetrahedra, vector<int>(4));   
+   for (size_t i = 0; i < out.numberoftetrahedra; i++)
+   {
+     for (size_t j = 0; j < 4; j++)
+       N[i][j] = out.neighborlist[i * 4 + j];
+   } 
+  
+   // extract point 2 tetrahedron list 
+   PT.resize(out.numberofpoints, vector<int>(1));
+   for (size_t i = 0; i < out.numberofpoints; i++)
+   {
+     PT[i][0] = out.point2tetlist[i]; 
+   }	
+ 
+   //extract face to tetrahedron list
+   FT.resize(out.numberoftrifaces, vector<int>(2)); 
+   int triface;
+   
+   for (size_t i = 0; i < out.numberoftrifaces; i++)
+   {
+     for (size_t j = 0; j < 2; j++)
+     {
+       FT[i][j] = out.face2tetlist[0]; 
+     }
+   }
+
+   return true;
+}
+
+
 IGL_INLINE bool igl::copyleft::tetgen::tetgenio_to_tetmesh(
   const tetgenio & out,
   std::vector<std::vector<REAL > > & V, 

include/igl/copyleft/tetgen/tetgenio_to_tetmesh.h

@@ -14,6 +14,7 @@
 #endif
 #include "tetgen.h" // Defined tetgenio, REAL
 #include <vector>
+#include <unordered_map>
 #include <Eigen/Core>
 namespace igl
 {
@@ -54,7 +55,32 @@ namespace igl
         const tetgenio & out,
         Eigen::PlainObjectBase<DerivedV>& V,
         Eigen::PlainObjectBase<DerivedT>& T);
-    }
+    
+      // Extract a tetrahedral mesh from a tetgenio object
+      // Inputs:
+      //   out tetgenio output object
+      // Outputs:
+      //   V  #V by 3 vertex position list
+      //   T  #T by 4 list of tetrahedra indices into V
+      //   F  #F by 3 list of marked facets
+      //   R  #T list of region IDs for tetrahedra
+      //   N  #T by 2 list of neighbors for each tetrahedron
+      //   PT #V list of incident tetrahedron for each vertex
+      //   FT #F by 2 list of tetrahedra sharing each face 
+      //   nR number of regions in output mesh
+      // Returns true on success, false on error
+      IGL_INLINE bool tetgenio_to_tetmesh(
+        const tetgenio & out,
+	std::vector<std::vector<REAL > > & V,
+	std::vector<std::vector<int> > & T,
+        std::vector<std::vector<int> > & F, 
+	std::vector<std::vector<REAL> > & R,// region marks for tetrahedrons
+	std::vector<std::vector<int > > &N, // neighborlist per tet
+	std::vector<std::vector<int > >	&PT, // Point to tet list per point
+	std::vector<std::vector<int > > &FT, // face to tet list
+	size_t & nR); // number of regions    
+
+  }
   }
 }
 

include/igl/copyleft/tetgen/tetrahedralize.cpp

@@ -18,6 +18,131 @@
 #include <cassert>
 #include <iostream>
 
+
+IGL_INLINE int igl::copyleft::tetgen::tetrahedralize(
+  const std::vector<std::vector<REAL > > & V,
+  const std::vector<std::vector<int> > & F,
+  const std::vector<std::vector<REAL > > & H,
+  const std::vector<std::vector<REAL > > & R,
+  const std::string switches,
+  std::vector<std::vector<REAL > > & TV,
+  std::vector<std::vector<int > > & TT,
+  std::vector<std::vector<int > > & TF,
+  std::vector<std::vector<REAL > > &TR,
+  std::vector<std::vector<int > > & TN, 
+  std::vector<std::vector<int > > & PT,
+  std::vector<std::vector<int > > & FT,
+  size_t & numRegions)
+{
+  using namespace std;
+  tetgenio in,out;
+  bool success;
+  success = mesh_to_tetgenio(V, F, H, R, in);
+  if(!success)
+  {
+    return -1;
+  }
+  try 
+  {
+    char * cswitches = new char[switches.size() + 1];
+    strcpy(cswitches, switches.c_str());
+    
+    ::tetrahedralize(cswitches, &in, &out); 
+    delete[] cswitches;
+}catch(int e)
+  {
+    cerr <<"^"<<__FUNCTION__<<": TETGEN CRASHED...KABOOM!!"<<endl;
+    return 1;
+  }
+  if(out.numberoftetrahedra == 0)
+  {
+    cerr<<"^"<<__FUNCTION__<<": Tetgen failed to create tets"<<endl;
+    return 2;	  
+  }
+   success = tetgenio_to_tetmesh(out, TV, TT, TF, TR, TN, PT, FT, numRegions
+);
+  if(!success)
+  {
+    return -1;
+  }
+    return 0;
+}
+
+template <
+  typename DerivedV, 
+  typename DerivedF, 
+  typename DerivedH,
+  typename DerivedR,
+  typename DerivedTV, 
+  typename DerivedTT, 
+  typename DerivedTF,
+  typename DerivedTR>
+IGL_INLINE int igl::copyleft::tetgen::tetrahedralize(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  const Eigen::PlainObjectBase<DerivedH>& H,
+  const Eigen::PlainObjectBase<DerivedR>& R,
+  const std::string switches,
+  Eigen::PlainObjectBase<DerivedTV>& TV,
+  Eigen::PlainObjectBase<DerivedTT>& TT,
+  Eigen::PlainObjectBase<DerivedTF>& TF,
+  Eigen::PlainObjectBase<DerivedTR>& TR,
+  Eigen::PlainObjectBase<DerivedTT>& TN,
+  Eigen::PlainObjectBase<DerivedTT>& PT,
+  Eigen::PlainObjectBase<DerivedTT>& FT,
+  size_t & numRegions)
+{
+  using namespace std;
+  vector<vector<REAL> > vV, vH, vR, vTV, vTR;
+  vector<vector<int> > vF,vTT,vTF, vTN, vPT, vFT;
+  matrix_to_list(V,vV);
+  matrix_to_list(F,vF);	
+  matrix_to_list(H, vH);
+  matrix_to_list(R, vR);
+  
+  int e = tetrahedralize(vV,vF,vH,vR,switches,vTV,vTT,vTF,vTR,vTN,vPT,vFT, numRegions);
+  
+  if(e == 0)
+  {
+    bool TV_rect = list_to_matrix(vTV,TV);
+    if(!TV_rect)
+    {
+      return 3;
+    }
+    bool TT_rect = list_to_matrix(vTT,TT);
+    if(!TT_rect)
+    {
+      return 3;
+    }
+    bool TF_rect = list_to_matrix(vTF,TF);
+    if(!TF_rect)
+    {
+      return 3;
+    }
+    bool TR_rect = list_to_matrix(vTR, TR);
+    if(!TR_rect)
+    {
+      return 3;	    
+    }
+    bool TN_rect = list_to_matrix(vTN, TN);
+    if(!TN_rect)
+    {
+      return 3;	
+    }	    
+    bool PT_rect = list_to_matrix(vPT, PT);
+    if(!PT_rect)
+    {
+      return 3;	    
+    }	    
+    bool FT_rect = list_to_matrix(vFT, FT);
+    if(!FT_rect)
+    {
+      return 3;	    
+    }	    
+  }
+  return e;
+}
+
 IGL_INLINE int igl::copyleft::tetgen::tetrahedralize(
   const std::vector<std::vector<REAL > > & V,
   const std::vector<std::vector<int> > & F,

include/igl/copyleft/tetgen/tetrahedralize.h

@@ -123,7 +123,83 @@ namespace igl
         Eigen::PlainObjectBase<DerivedTT>& TT,
         Eigen::PlainObjectBase<DerivedTF>& TF,
         Eigen::PlainObjectBase<DerivedTM>& TM);
-    }
+	
+
+      // Mesh the interior of a surface mesh (V,F) using tetgen
+      //
+      // Inputs:
+      //   V  #V by 3 vertex position list
+      //   F  #F list of polygon face indices into V (0-indexed)
+      //   H  #H by 3 list of seed points inside holes
+      //   R  #R by 5 list of region attributes		
+
+      //   switches  string of tetgen options (See tetgen documentation) e.g.
+      //     "pq1.414a0.01" tries to mesh the interior of a given surface with
+      //       quality and area constraints
+      //     "" will mesh the convex hull constrained to pass through V (ignores F)
+      // Outputs:
+      //   TV  #V by 3 vertex position list
+      //   TT  #T by 4 list of tet face indices
+      //   TF  #F by 3 list of triangle face indices
+      //   TR  #T list of region ID for each tetrahedron	
+      //   TN  #T by 4 list of indices neighbors for each tetrahedron
+      //   PT  #V list of incident tetrahedron for a vertex
+      //   FT  #F by 2 list of tetrahedrons sharing a triface	
+      //   numRegions Number of regions in output mesh
+
+      // Returns status:
+      //   0 success
+      //   1 tetgen threw exception
+      //   2 tetgen did not crash but could not create any tets (probably there are
+      //     holes, duplicate faces etc.)
+      //   -1 other error
+	IGL_INLINE int tetrahedralize(
+	  const std::vector<std::vector<REAL> > &V, 
+	  const std::vector<std::vector<int> >  &F, 
+	  const std::vector<std::vector<REAL> > &H, 
+	  const std::vector<std::vector<REAL> > &R, 
+	  	  
+	  const std::string switches, 
+	  
+	  std::vector<std::vector<REAL > > & TV,
+	  std::vector<std::vector<int > >  & TT,
+	  std::vector<std::vector<int > >  & TF,
+	  std::vector<std::vector<REAL > > &TR,  
+	  std::vector<std::vector<int > > &TN, 
+	  std::vector<std::vector<int > > &PT, 
+	  std::vector<std::vector<int > > &FT, 
+	  size_t & numRegions);	     
+
+
+      // Wrapper with Eigen types
+      // Templates:
+      //   DerivedV  real-value: i.e. from MatrixXd
+      //   DerivedF  integer-value: i.e. from MatrixXi	
+      template <
+        typename DerivedV,
+   	typename DerivedF,
+	typename DerivedH,
+	typename DerivedR,
+	typename DerivedTV,
+	typename DerivedTT,
+	typename DerivedTF,
+	typename DerivedTR>	
+      IGL_INLINE int tetrahedralize(
+	const Eigen::PlainObjectBase<DerivedV>& V,
+	const Eigen::PlainObjectBase<DerivedF>& F,
+	const Eigen::PlainObjectBase<DerivedH>& H,
+	const Eigen::PlainObjectBase<DerivedR>& R,
+	const std::string switches,
+	Eigen::PlainObjectBase<DerivedTV>& TV,
+        Eigen::PlainObjectBase<DerivedTT>& TT,
+	Eigen::PlainObjectBase<DerivedTF>& TF,
+	Eigen::PlainObjectBase<DerivedTR>& TR, 
+	Eigen::PlainObjectBase<DerivedTT>& TN, 
+	Eigen::PlainObjectBase<DerivedTT>& PT, 
+	Eigen::PlainObjectBase<DerivedTT>& FT, 
+	size_t & numRegions);	      
+
+   }
   }
 }