bug fix in winding number

Former-commit-id: bb9df666147576166352007f80d959f302e2651e

RELEASE_HISTORY.txt

@@ -1,3 +1,5 @@
+1.0.2  Bug fix in winding number code
+1.0.1  Bug fixes and more CGAL support
 1.0.0  Major beta release: many renames, tutorial, triangle wrapper, org. build
 0.4.6  Generalized Winding Numbers
 0.4.5  CGAL extra: mesh selfintersection

VERSION.txt

@@ -3,4 +3,4 @@
 # Anyone may increment Minor to indicate a small change.
 # Major indicates a large change or large number of changes (upload to website)
 # World indicates a substantial change or release
-1.0.1
+1.0.2

examples/intersections/example.cpp

@@ -7,6 +7,7 @@
 #include <igl/quat_to_mat.h>
 #include <igl/report_gl_error.h>
 #include <igl/readOBJ.h>
+#include <igl/writeOBJ.h>
 #include <igl/readDMAT.h>
 #include <igl/readOFF.h>
 #include <igl/readMESH.h>
@@ -273,6 +274,7 @@ void display()
     const MatrixXd & N,
     const MatrixXd & C)
   {
+    glEnable(GL_COLOR_MATERIAL);
     glEnable(GL_POLYGON_OFFSET_FILL); // Avoid Stitching!
     glPolygonOffset(1.0,1);
     glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
@@ -476,12 +478,14 @@ void color_selfintersections(
 {
   using namespace igl;
   using namespace Eigen;
+  using namespace std;
   MatrixXd SV;
   MatrixXi SF,IF;
   VectorXi J,IM;
   RemeshSelfIntersectionsParam params;
-  params.detect_only = true;
+  params.detect_only = false;
   remesh_self_intersections(V,F,params,SV,SF,IF,J,IM);
+  writeOBJ("FUCK.obj",SV,SF);
   C.resize(F.rows(),3);
   C.col(0).setConstant(0.4);
   C.col(1).setConstant(0.8);
@@ -639,6 +643,7 @@ int main(int argc, char * argv[])
   }else
   {
     VU = V;
+    color_selfintersections(V,F,C);
   }
   mid = 0.5*(VU.colwise().maxCoeff() + VU.colwise().minCoeff());
   bbd = (VU.colwise().maxCoeff() - VU.colwise().minCoeff()).maxCoeff();

examples/skeleton-posing/example.cpp

@@ -38,6 +38,7 @@
 #include <igl/winding_number.h>
 #include <igl/writeDMAT.h>
 #include <igl/writeOBJ.h>
+#include <igl/writeMESH.h>
 #include <igl/writeOFF.h>
 #include <igl/writeTGF.h>
 #include <igl/next_filename.h>
@@ -59,6 +60,8 @@
 #include <iostream>
 #include <iomanip>
 
+#define VERBOSE
+
 enum SkelStyleType
 {
   SKEL_STYLE_TYPE_3D = 0,
@@ -698,8 +701,16 @@ bool clean(
   using namespace std;
   {
     MatrixXi _1;
-    VectorXi _2,_3;
-    remesh_self_intersections(V,F,{},CV,CF,_1,_2,_3);
+    VectorXi _2,IM;
+#ifdef VERBOSE
+    cout<<"remesh_self_intersections"<<endl;
+#endif
+    remesh_self_intersections(V,F,{},CV,CF,_1,_2,IM);
+    for_each(CF.data(),CF.data()+CF.size(),[&IM](int & a){a=IM(a);});
+    MatrixXd oldCV = CV;
+    MatrixXi oldCF = CF;
+    remove_unreferenced(oldCV,oldCF,CV,CF,IM);
+    writeOBJ("remesh.obj",CV,CF);
   }
   MatrixXd TV;
   MatrixXi TT;
@@ -708,17 +719,26 @@ bool clean(
     // c  convex hull
     // Y  no boundary steiners
     // p  polygon input
-    if(tetrahedralize(V,F,"cYp",TV,TT,_1) != 0)
+#ifdef VERBOSE
+    cout<<"tetrahedralize"<<endl;
+#endif
+    if(tetrahedralize(CV,CF,"cYpC",TV,TT,_1) != 0)
     {
       cout<<REDRUM("CDT failed.")<<endl;
       return false;
     }
+    writeMESH("tetrahedralize-1.mesh",TV,TT,MatrixXi());
   }
   {
     MatrixXd BC;
     barycenter(TV,TT,BC);
     VectorXd W;
+#ifdef VERBOSE
+    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();
@@ -733,6 +753,7 @@ bool clean(
     }
     assert(c==count);
     boundary_facets(CT,CF);
+    writeOBJ("boundary_facets.obj",CV,CF);
   }
   return true;
 }
@@ -759,11 +780,16 @@ bool robust_weights(
   // compute tet-mesh
   {
     MatrixXi _1;
-    if(!mesh_with_skeleton(V,F,C,{},BE,{},10,TV,TT,_1))
+#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))
     {
       cout<<REDRUM("tetgen failed.")<<endl;
       return false;
     }
+    writeMESH("after.mesh",TV,TT,MatrixXi());
   }
   // compute weights
   VectorXi b;

include/igl/WindingNumberAABB.h

@@ -214,7 +214,10 @@ inline bool igl::WindingNumberAABB<Point>::inside(const Point & p) const
   assert(p.size() == min_corner.size());
   for(int i = 0;i<p.size();i++)
   {
-    if( p(i) < min_corner(i) || p(i) >= max_corner(i))
+    //// Perfect matching is **not** robust
+    //if( p(i) < min_corner(i) || p(i) >= max_corner(i))
+    // **MUST** be conservative!!
+    if( p(i) < min_corner(i) || p(i) > max_corner(i))
     {
       return false;
     }

include/igl/avg_edge_length.cpp

@@ -19,10 +19,10 @@ IGL_INLINE double igl::avg_edge_length(
 
   for (unsigned i=0;i<F.rows();++i)
   {
-    for (unsigned j=0;j<3;++j)
+    for (unsigned j=0;j<F.cols();++j)
     {
       ++count;
-      avg += (V.row(F(i,j)) - V.row(F(i,(j+1)%3))).norm();
+      avg += (V.row(F(i,j)) - V.row(F(i,(j+1)%F.cols()))).norm();
     }
   }
 

include/igl/avg_edge_length.h

@@ -22,7 +22,7 @@ namespace igl
   //   DerivedL derived from edge lengths matrix type: i.e. MatrixXd
   // Inputs:
   //   V  eigen matrix #V by 3
-  //   F  #F by 3 list of mesh faces (must be triangles)
+  //   F  #F by simplex-size list of mesh faces (must be simplex)
   // Outputs:
   //   l  average edge length
   //

include/igl/edge_lengths.cpp

@@ -17,6 +17,15 @@ IGL_INLINE void igl::edge_lengths(
   using namespace std;
   switch(F.cols())
   {
+    case 2:
+    {
+      L.resize(F.rows(),1);
+      for(int i = 0;i<F.rows();i++)
+      {
+        L(i,0) = (V.row(F(i,1))-V.row(F(i,0))).norm();
+      }
+      break;
+    }
     case 3:
     {
       L.resize(F.rows(),3);

include/igl/edge_lengths.h

@@ -22,11 +22,14 @@ namespace igl
   //   DerivedL derived from edge lengths matrix type: i.e. MatrixXd
   // Inputs:
   //   V  eigen matrix #V by 3
+  //   F  #F by 2 list of mesh edges
+  //    or
   //   F  #F by 3 list of mesh faces (must be triangles)
   //    or
   //   T  #T by 4 list of mesh elements (must be tets)
   // Outputs:
-  //   L  #F by {3|6} list of edge lengths 
+  //   L  #F by {1|3|6} list of edge lengths 
+  //     for edges, column of lengths
   //     for triangles, columns correspond to edges [1,2],[2,0],[0,1]
   //     for tets, columns correspond to edges
   //     [1,2],[2,0],[0,1],[3,0],[3,1],[3,2]

include/igl/matlab/mexStream.h

@@ -30,7 +30,6 @@ namespace igl
       inline virtual int overflow(int c = EOF);
   }; 
 }
-#ifndef IGL_STATIC_LIBRARY
-#  include "MexStream.cpp"
-#endif
+// ALWAYS INCLUDE
+#include "MexStream.cpp"
 #endif

include/igl/tetgen/mesh_with_skeleton.cpp

@@ -11,6 +11,7 @@
 #include <igl/cat.h>
 #include <igl/tetgen/tetrahedralize.h>
 #include <igl/writeOFF.h>
+#include <igl/writeOBJ.h>
 
 #include <iostream>
 // Default settings pq2Y tell tetgen to mesh interior of triangle mesh and
@@ -43,15 +44,13 @@ IGL_INLINE bool igl::mesh_with_skeleton(
   sample_edges(C,BECE,samples_per_bone,S);
   // Vertices we'll constrain tet mesh to meet
   MatrixXd VS = cat(1,V,S);
-  // Boundary faces
-  MatrixXi BF;
   // Use tetgen to mesh the interior of surface, this assumes surface:
   //   * has no holes
   //   * has no non-manifold edges or vertices
   //   * has consistent orientation
   //   * has no self-intersections
   //   * has no 0-volume pieces
-  //writeOFF("mesh_with_skeleton.off",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;

include/igl/triangle_fan.cpp

@@ -20,7 +20,8 @@ IGL_INLINE void igl::triangle_fan(
   // outgoing (left or right). Thus this will not work.
   assert(E.cols() == 2);
   // Arbitrary starting vertex
-  int s = E(int(((double)rand() / RAND_MAX)*E.rows()),0);
+  //int s = E(int(((double)rand() / RAND_MAX)*E.rows()),0);
+  int s = E(rand()%E.rows(),0);
   vector<vector<int> >  lcap;
   for(int i = 0;i<E.rows();i++)
   {

include/igl/writeMESH.cpp

@@ -95,7 +95,7 @@ IGL_INLINE bool igl::writeMESH(
   for(int i = 0;i<number_of_tet_vertices;i++)
   {
     // print position of ith tet vertex
-    fprintf(mesh_file,"%lg %lg %lg 1\n",
+    fprintf(mesh_file,"%.17lg %.17lg %.17lg 1\n",
       (double)V(i,0),
       (double)V(i,1),
       (double)V(i,2));

include/igl/writeOBJ.cpp

@@ -8,6 +8,8 @@
 #include "writeOBJ.h"
 
 #include <iostream>
+#include <limits>
+#include <iomanip>
 #include <fstream>
 #include <cstdio>
 
@@ -18,6 +20,7 @@ IGL_INLINE bool igl::writeOBJ(
   const Eigen::PlainObjectBase<DerivedF>& F)
 {
   std::ofstream s(str.c_str());
+  s.precision(std::numeric_limits<double>::digits10 + 1);
 
   if(!s.is_open())
   {