From b4eee2bc3d2b6db9e1a1354479aeee6cacfbbbf4 Mon Sep 17 00:00:00 2001
From: Alec Jacobson <alecjacobson@gmail.com>
Date: Tue, 9 Sep 2014 01:43:45 -0400
Subject: [PATCH] more robust weight computation in skeleton poser example

Former-commit-id: 8309c73bd9e2fd6139df98f6238091f7fa0060ec
---
 examples/skeleton-posing/example.cpp      | 30 +++++++++++++++++------
 include/igl/tetgen/mesh_with_skeleton.cpp |  2 +-
 include/igl/volume.cpp                    |  9 ++++---
 3 files changed, 28 insertions(+), 13 deletions(-)

diff --git a/examples/skeleton-posing/example.cpp b/examples/skeleton-posing/example.cpp
index 1445e5215..427e39920 100644
--- a/examples/skeleton-posing/example.cpp
+++ b/examples/skeleton-posing/example.cpp
@@ -42,6 +42,7 @@
 #include <igl/writeOFF.h>
 #include <igl/writeTGF.h>
 #include <igl/next_filename.h>
+#include <igl/volume.h>
 
 #include <Eigen/Core>
 #include <Eigen/Geometry>
@@ -710,7 +711,6 @@ bool clean(
     MatrixXd oldCV = CV;
     MatrixXi oldCF = CF;
     remove_unreferenced(oldCV,oldCF,CV,CF,IM);
-    writeOBJ("remesh.obj",CV,CF);
   }
   MatrixXd TV;
   MatrixXi TT;
@@ -727,7 +727,6 @@ bool clean(
       cout<<REDRUM("CDT failed.")<<endl;
       return false;
     }
-    writeMESH("tetrahedralize-1.mesh",TV,TT,MatrixXi());
   }
   {
     MatrixXd BC;
@@ -737,8 +736,6 @@ bool clean(
     cout<<"winding_number"<<endl;
 #endif
     winding_number(V,F,BC,W);
-    writeDMAT("W.dmat",W);
-    writeDMAT("BC.dmat",BC);
     W = W.array().abs();
     const double thresh = 0.5;
     const int count = (W.array()>thresh).cast<int>().sum();
@@ -753,7 +750,6 @@ bool clean(
     }
     assert(c==count);
     boundary_facets(CT,CF);
-    writeOBJ("boundary_facets.obj",CV,CF);
   }
   return true;
 }
@@ -783,14 +779,32 @@ bool robust_weights(
 #ifdef VERBOSE
     cout<<"mesh_with_skeleton"<<endl;
 #endif
-    writeOBJ("before.obj",CV,CF);
-    if(!mesh_with_skeleton(CV,CF,C,{},BE,{},10,TV,TT,_1))
+    if(!mesh_with_skeleton(CV,CF,C,{},BE,{},10,"pq1.5Y",TV,TT,_1))
     {
       cout<<REDRUM("tetgen failed.")<<endl;
       return false;
     }
-    writeMESH("after.mesh",TV,TT,MatrixXi());
   }
+  // Finally, tetgen may have still included some insanely small tets.
+  // Just ignore these during weight computation (and hope they don't isolate
+  // any vertices).
+  {
+    const MatrixXi oldTT = TT;
+    VectorXd vol;
+    volume(TV,TT,vol);
+    const double thresh = 1e-17;
+    const int count = (vol.array()>thresh).cast<int>().sum();
+    TT.resize(count,oldTT.cols());
+    int c = 0;
+    for(int t = 0;t<oldTT.rows();t++)
+    {
+      if(vol(t)>thresh)
+      {
+        TT.row(c++) = oldTT.row(t);
+      }
+    }
+  }
+
   // compute weights
   VectorXi b;
   MatrixXd bc;
diff --git a/include/igl/tetgen/mesh_with_skeleton.cpp b/include/igl/tetgen/mesh_with_skeleton.cpp
index 7fc831936..b6038b18d 100644
--- a/include/igl/tetgen/mesh_with_skeleton.cpp
+++ b/include/igl/tetgen/mesh_with_skeleton.cpp
@@ -50,7 +50,7 @@ IGL_INLINE bool igl::mesh_with_skeleton(
   //   * has consistent orientation
   //   * has no self-intersections
   //   * has no 0-volume pieces
-  writeOBJ("mesh_with_skeleton.obj",VS,F);
+  //writeOBJ("mesh_with_skeleton.obj",VS,F);
   cerr<<"tetgen begin()"<<endl;
   int status = tetrahedralize( VS,F,eff_tetgen_flags,VV,TT,FF);
   cerr<<"tetgen end()"<<endl;
diff --git a/include/igl/volume.cpp b/include/igl/volume.cpp
index fc6288225..133303bb4 100644
--- a/include/igl/volume.cpp
+++ b/include/igl/volume.cpp
@@ -22,10 +22,10 @@ IGL_INLINE void igl::volume(
   vol.resize(m,1);
   for(int t = 0;t<m;t++)
   {
-    const auto & a = V.row(T(t,0));
-    const auto & b = V.row(T(t,1));
-    const auto & c = V.row(T(t,2));
-    const auto & d = V.row(T(t,3));
+    const RowVector3d & a = V.row(T(t,0));
+    const RowVector3d & b = V.row(T(t,1));
+    const RowVector3d & c = V.row(T(t,2));
+    const RowVector3d & d = V.row(T(t,3));
     vol(t) = -(a-d).dot((b-d).cross(c-d))/6.;
   }
 }
@@ -111,4 +111,5 @@ template void igl::volume<Eigen::Matrix<double, -1, 6, 0, -1, 6>, Eigen::Matrix<
 template Eigen::Matrix<double, 1, 3, 1, 1, 3>::Scalar igl::volume_single<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, Eigen::Matrix<double, 1, 3, 1, 1, 3> const&);
 template Eigen::Matrix<double, 3, 1, 0, 3, 1>::Scalar igl::volume_single<Eigen::Matrix<double, 3, 1, 0, 3, 1>, Eigen::Matrix<double, 3, 1, 0, 3, 1>, Eigen::Matrix<double, 3, 1, 0, 3, 1>, Eigen::Matrix<double, 3, 1, 0, 3, 1> >(Eigen::Matrix<double, 3, 1, 0, 3, 1> const&, Eigen::Matrix<double, 3, 1, 0, 3, 1> const&, Eigen::Matrix<double, 3, 1, 0, 3, 1> const&, Eigen::Matrix<double, 3, 1, 0, 3, 1> const&);
 template void igl::volume<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::volume<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 #endif