Merge pull request #349 from yig/master

seam_edges.cpp fixed duplicate boundary edges bug and removed helper

Former-commit-id: b29647e053ff8e81b657a6d7212f878b8ac792b3

include/igl/point_simplex_squared_distance.cpp

@@ -20,26 +20,29 @@ template <
   typename DerivedV,
   typename DerivedEle,
   typename Derivedsqr_d,
-  typename Derivedc>
+  typename Derivedc,
+  typename Derivedb>
 IGL_INLINE void igl::point_simplex_squared_distance(
   const Eigen::MatrixBase<Derivedp> & p,
   const Eigen::MatrixBase<DerivedV> & V,
   const Eigen::MatrixBase<DerivedEle> & Ele,
   const typename DerivedEle::Index primitive,
   Derivedsqr_d & sqr_d,
-  Eigen::PlainObjectBase<Derivedc> & c)
+  Eigen::PlainObjectBase<Derivedc> & c,
+  Eigen::PlainObjectBase<Derivedb> & bary)
 {
   typedef typename Derivedp::Scalar Scalar;
   typedef typename Eigen::Matrix<Scalar,1,DIM> Vector;
   typedef Vector Point;
+  typedef Eigen::PlainObjectBase<Derivedb> BaryPoint;
 
   const auto & Dot = [](const Point & a, const Point & b)->Scalar
   {
     return a.dot(b);
   };
   // Real-time collision detection, Ericson, Chapter 5
-  const auto & ClosestPtPointTriangle = 
-    [&Dot](Point p, Point a, Point b, Point c)->Point 
+  const auto & ClosestBaryPtPointTriangle = 
+    [&Dot](Point p, Point a, Point b, Point c, BaryPoint& bary_out )->Point 
   {
     // Check if P in vertex region outside A
     Vector ab = b - a;
@@ -47,42 +50,61 @@ IGL_INLINE void igl::point_simplex_squared_distance(
     Vector ap = p - a;
     Scalar d1 = Dot(ab, ap);
     Scalar d2 = Dot(ac, ap);
-    if (d1 <= 0.0 && d2 <= 0.0) return a; // barycentric coordinates (1,0,0)
+    if (d1 <= 0.0 && d2 <= 0.0) {
+      // barycentric coordinates (1,0,0)
+      bary_out << 1, 0, 0;
+      return a;
+    }
     // Check if P in vertex region outside B
     Vector bp = p - b;
     Scalar d3 = Dot(ab, bp);
     Scalar d4 = Dot(ac, bp);
-    if (d3 >= 0.0 && d4 <= d3) return b; // barycentric coordinates (0,1,0)
+    if (d3 >= 0.0 && d4 <= d3) {
+      // barycentric coordinates (0,1,0)
+      bary_out << 0, 1, 0;
+      return b;
+    }
     // Check if P in edge region of AB, if so return projection of P onto AB
     Scalar vc = d1*d4 - d3*d2;
     if( a != b)
     {
       if (vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0) {
         Scalar v = d1 / (d1 - d3);
-        return a + v * ab; // barycentric coordinates (1-v,v,0)
+        // barycentric coordinates (1-v,v,0)
+        bary_out << 1-v, v, 0;
+        return a + v * ab;
       }
     }
     // Check if P in vertex region outside C
     Vector cp = p - c;
     Scalar d5 = Dot(ab, cp);
     Scalar d6 = Dot(ac, cp);
-    if (d6 >= 0.0 && d5 <= d6) return c; // barycentric coordinates (0,0,1)
+    if (d6 >= 0.0 && d5 <= d6) {
+      // barycentric coordinates (0,0,1)
+      bary_out << 0, 0, 1;
+      return c;
+    }
     // Check if P in edge region of AC, if so return projection of P onto AC
     Scalar vb = d5*d2 - d1*d6;
     if (vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0) {
       Scalar w = d2 / (d2 - d6);
-      return a + w * ac; // barycentric coordinates (1-w,0,w)
+      // barycentric coordinates (1-w,0,w)
+      bary_out << 1-w, 0, w;
+      return a + w * ac;
     }
     // Check if P in edge region of BC, if so return projection of P onto BC
     Scalar va = d3*d6 - d5*d4;
     if (va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0) {
       Scalar w = (d4 - d3) / ((d4 - d3) + (d5 - d6));
-      return b + w * (c - b); // barycentric coordinates (0,1-w,w)
+      // barycentric coordinates (0,1-w,w)
+      bary_out << 0, 1-w, w;
+      return b + w * (c - b);
     }
     // P inside face region. Compute Q through its barycentric coordinates (u,v,w)
     Scalar denom = 1.0 / (va + vb + vc);
     Scalar v = vb * denom;
     Scalar w = vc * denom;
+    bary_out << 1.0-v-w, v, w;
     return a + ab * v + ac * w; // = u*a + v*b + w*c, u = va * denom = 1.0-v-w
   };
 
@@ -91,16 +113,38 @@ IGL_INLINE void igl::point_simplex_squared_distance(
   assert(Ele.cols() <= DIM+1);
   assert(Ele.cols() <= 3 && "Only simplices up to triangles are considered");
 
-  c = ClosestPtPointTriangle(
+  bary.resize( DIM, 1 );
+  c = ClosestBaryPtPointTriangle(
     p,
     V.row(Ele(primitive,0)),
     V.row(Ele(primitive,1%Ele.cols())),
-    V.row(Ele(primitive,2%Ele.cols())));
+    V.row(Ele(primitive,2%Ele.cols())),
+    bary);
   sqr_d = (p-c).squaredNorm();
 }
+template <
+  int DIM,
+  typename Derivedp,
+  typename DerivedV,
+  typename DerivedEle,
+  typename Derivedsqr_d,
+  typename Derivedc>
+IGL_INLINE void igl::point_simplex_squared_distance(
+  const Eigen::MatrixBase<Derivedp> & p,
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedEle> & Ele,
+  const typename DerivedEle::Index primitive,
+  Derivedsqr_d & sqr_d,
+  Eigen::PlainObjectBase<Derivedc> & c)
+{
+    Eigen::PlainObjectBase<Derivedc> b;
+    point_simplex_squared_distance<DIM>( p, V, Ele, primitive, sqr_d, c, b );
+}
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instanciation
 template void igl::point_simplex_squared_distance<3, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, double&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
 template void igl::point_simplex_squared_distance<2, Eigen::Matrix<double, 1, 2, 1, 1, 2>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<double, 1, 2, 1, 1, 2> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, double&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> >&);
+template void igl::point_simplex_squared_distance<3, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, double&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
+template void igl::point_simplex_squared_distance<2, Eigen::Matrix<double, 1, 2, 1, 1, 2>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<double, 1, 2, 1, 1, 2> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, double&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> >&);
 #endif

include/igl/point_simplex_squared_distance.h

@@ -36,6 +36,35 @@ namespace igl
     const typename DerivedEle::Index i,
     Derivedsqr_d & sqr_d,
     Eigen::PlainObjectBase<Derivedc> & c);
+  // Determine squared distance from a point to linear simplex.
+  // Also return barycentric coordinate of closest point. 
+  //
+  // Inputs:
+  //   p  d-long query point
+  //   V  #V by d list of vertices
+  //   Ele  #Ele by ss<=d+1 list of simplex indices into V
+  //   i  index into Ele of simplex
+  // Outputs:
+  //   sqr_d  squared distance of Ele(i) to p
+  //   c  closest point on Ele(i) 
+  //   b  barycentric coordinates of closest point on Ele(i) 
+  //
+  template <
+    int DIM,
+    typename Derivedp,
+    typename DerivedV,
+    typename DerivedEle,
+    typename Derivedsqr_d,
+    typename Derivedc,
+    typename Derivedb>
+  IGL_INLINE void point_simplex_squared_distance(
+    const Eigen::MatrixBase<Derivedp> & p,
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedEle> & Ele,
+    const typename DerivedEle::Index i,
+    Derivedsqr_d & sqr_d,
+    Eigen::PlainObjectBase<Derivedc> & c,
+    Eigen::PlainObjectBase<Derivedb> & b);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "point_simplex_squared_distance.cpp"

include/igl/seam_edges.cpp

@@ -11,23 +11,6 @@
 #include <unordered_set>
 #include <cassert>
 
-namespace {
-    // Computes the orientation of `c` relative to the line between `a` and `b`.
-    // Assumes 2D vector input.
-    // Based on: https://www.cs.cmu.edu/~quake/robust.html
-    template< typename scalar_t >
-    inline scalar_t orientation(
-        const Eigen::Matrix< scalar_t, 2, 1 >& a,
-        const Eigen::Matrix< scalar_t, 2, 1 >& b,
-        const Eigen::Matrix< scalar_t, 2, 1 >& c
-        ) {
-        typedef Eigen::Matrix< scalar_t, 2, 1 > Vector2S;
-        const Vector2S row0 = a - c;
-        const Vector2S row1 = b - c;
-        return row0(0)*row1(1) - row1(0)*row0(1);
-    }
-}
-
 // I have verified that this function produces the exact same output as
 // `find_seam_fast.py` for `cow_triangled.obj`.
 template <typename DerivedV, typename DerivedF, typename DerivedT>
@@ -49,6 +32,15 @@ IGL_INLINE void igl::seam_edges(
     // Assume 2D texture coordinates (foldovers tests).
     assert( TC.cols() == 2 );
     typedef Eigen::Matrix< typename DerivedT::Scalar, 2, 1 > Vector2S;
+    // Computes the orientation of `c` relative to the line between `a` and `b`.
+    // Assumes 2D vector input.
+    // Based on: https://www.cs.cmu.edu/~quake/robust.html
+    const auto& Orientation = []( const Vector2S& a, const Vector2S& b, const Vector2S& c ) -> typename DerivedT::Scalar
+    {
+        const Vector2S row0 = a - c;
+        const Vector2S row1 = b - c;
+        return row0(0)*row1(1) - row1(0)*row0(1);
+    };
     
     seams     .setZero( 3*F.rows(), 4 );
     boundaries.setZero( 3*F.rows(), 2 );
@@ -73,7 +65,7 @@ IGL_INLINE void igl::seam_edges(
     typedef std::pair< typename DerivedF::Scalar, typename DerivedF::Scalar > directed_edge;
 	const int numV = V.rows();
 	const int numF = F.rows();
-	auto edge_hasher = [numV]( directed_edge const& e ) { return e.first*numV + e.second; };
+	const auto& edge_hasher = [numV]( directed_edge const& e ) { return e.first*numV + e.second; };
 	// When we pass a hash function object, we also need to specify the number of buckets.
 	// The Euler characteristic says that the number of undirected edges is numV + numF -2*genus.
 	std::unordered_map< directed_edge, std::pair< int, int >, decltype( edge_hasher ) > directed_position_edge2face_position_index( 2*( numV + numF ), edge_hasher );
@@ -84,12 +76,13 @@ IGL_INLINE void igl::seam_edges(
         }
     }
     
-    // First find all undirected position edges (collect both orientations of
+    // First find all undirected position edges (collect a canonical orientation of
     // the directed edges).
     std::unordered_set< directed_edge, decltype( edge_hasher ) > undirected_position_edges( numV + numF, edge_hasher );
     for( const auto& el : directed_position_edge2face_position_index ) {
-        undirected_position_edges.insert( el.first );
-        undirected_position_edges.insert( std::make_pair( el.first.second, el.first.first ) );
+        // The canonical orientation is the one where the smaller of
+        // the two vertex indices is first.
+        undirected_position_edges.insert( std::make_pair( std::min( el.first.first, el.first.second ), std::max( el.first.first, el.first.second ) ) );
     }
     
     // Now we will iterate over all position edges.
@@ -97,6 +90,10 @@ IGL_INLINE void igl::seam_edges(
     // texcoord indices (or one doesn't exist at all in the case of a mesh
     // boundary).
     for( const auto& vp_edge : undirected_position_edges ) {
+        // We should only see canonical edges,
+        // where the first vertex index is smaller.
+        assert( vp_edge.first < vp_edge.second );
+        
         const auto vp_edge_reverse = std::make_pair( vp_edge.second, vp_edge.first );
         // If it and its opposite exist as directed edges, check if their
         // texture coordinate indices match.
@@ -108,12 +105,6 @@ IGL_INLINE void igl::seam_edges(
             // NOTE: They should never be equal.
             assert( forwards != backwards );
             
-            // NOTE: Non-matching seam edges will show up twice, once as
-            //       ( forwards, backwards ) and once as
-            //       ( backwards, forwards ). We don't need to examine both,
-            //       so continue only if forwards < backwards.
-            if( forwards < backwards ) continue;
-
             // If the texcoord indices match (are similarly flipped),
             // this edge is not a seam. It could be a foldover.
             if( std::make_pair( FTC( forwards.first, forwards.second ), FTC( forwards.first, ( forwards.second+1 ) % 3 ) )
@@ -129,8 +120,8 @@ IGL_INLINE void igl::seam_edges(
                 const Vector2S c_backwards = TC.row( FTC( backwards.first, (backwards.second+2) % 3 ) );
                 // If the opposite vertices' texture coordinates fall on the same side
                 // of the edge, we have a UV-space foldover.
-                const auto orientation_forwards = orientation( a, b, c_forwards );
-                const auto orientation_backwards = orientation( a, b, c_backwards );
+                const auto orientation_forwards = Orientation( a, b, c_forwards );
+                const auto orientation_backwards = Orientation( a, b, c_backwards );
                 if( ( orientation_forwards > 0 && orientation_backwards > 0 ) ||
                     ( orientation_forwards < 0 && orientation_backwards < 0 )
                     ) {