Boundary facets orientation (#2362)

* fix boundary_facets orientation + test

* rm print in tests

---------

Co-authored-by: Alec Jacobson <[email protected]>

include/igl/boundary_facets.cpp

@@ -47,21 +47,21 @@ IGL_INLINE void igl::boundary_facets(
       for(int i = 0; i< (int)T.rows();i++)
       {
         // get face in correct order
-        allF(i*simplex_size+0,0) = T(i,1);
+        allF(i*simplex_size+0,2) = T(i,1);
         allF(i*simplex_size+0,1) = T(i,3);
-        allF(i*simplex_size+0,2) = T(i,2);
+        allF(i*simplex_size+0,0) = T(i,2);
         // get face in correct order
-        allF(i*simplex_size+1,0) = T(i,0);
+        allF(i*simplex_size+1,2) = T(i,0);
         allF(i*simplex_size+1,1) = T(i,2);
-        allF(i*simplex_size+1,2) = T(i,3);
+        allF(i*simplex_size+1,0) = T(i,3);
         // get face in correct order
-        allF(i*simplex_size+2,0) = T(i,0);
+        allF(i*simplex_size+2,2) = T(i,0);
         allF(i*simplex_size+2,1) = T(i,3);
-        allF(i*simplex_size+2,2) = T(i,1);
+        allF(i*simplex_size+2,0) = T(i,1);
         // get face in correct order
-        allF(i*simplex_size+3,0) = T(i,0);
+        allF(i*simplex_size+3,2) = T(i,0);
         allF(i*simplex_size+3,1) = T(i,1);
-        allF(i*simplex_size+3,2) = T(i,2);
+        allF(i*simplex_size+3,0) = T(i,2);
       }
       break;
     case 3:
@@ -124,101 +124,12 @@ Ret igl::boundary_facets(
   return F;
 }
 
-template <typename IntegerT, typename IntegerF>
-IGL_INLINE void igl::boundary_facets(
-  const std::vector<std::vector<IntegerT> > & T,
-  std::vector<std::vector<IntegerF> > & F)
-{
-  // Kept for legacy reasons. Could probably just delete.
-  using namespace std;
-
-  if(T.size() == 0)
-  {
-    F.clear();
-    return;
-  }
-
-  int simplex_size = T[0].size();
-  // Get a list of all faces
-  vector<vector<IntegerF> > allF(
-    T.size()*simplex_size,
-    vector<IntegerF>(simplex_size-1));
-
-  // Gather faces, loop over tets
-  for(int i = 0; i< (int)T.size();i++)
-  {
-    assert((int)T[i].size() == simplex_size);
-    switch(simplex_size)
-    {
-      case 4:
-        // get face in correct order
-        allF[i*simplex_size+0][0] = T[i][1];
-        allF[i*simplex_size+0][1] = T[i][3];
-        allF[i*simplex_size+0][2] = T[i][2];
-        // get face in correct order
-        allF[i*simplex_size+1][0] = T[i][0];
-        allF[i*simplex_size+1][1] = T[i][2];
-        allF[i*simplex_size+1][2] = T[i][3];
-        // get face in correct order
-        allF[i*simplex_size+2][0] = T[i][0];
-        allF[i*simplex_size+2][1] = T[i][3];
-        allF[i*simplex_size+2][2] = T[i][1];
-        // get face in correct order
-        allF[i*simplex_size+3][0] = T[i][0];
-        allF[i*simplex_size+3][1] = T[i][1];
-        allF[i*simplex_size+3][2] = T[i][2];
-        break;
-      case 3:
-        allF[i*simplex_size+0][0] = T[i][1];
-        allF[i*simplex_size+0][1] = T[i][2];
-        allF[i*simplex_size+1][0] = T[i][2];
-        allF[i*simplex_size+1][1] = T[i][0];
-        allF[i*simplex_size+2][0] = T[i][0];
-        allF[i*simplex_size+2][1] = T[i][1];
-        break;
-    }
-  }
-
-  // Counts
-  vector<int> C;
-  face_occurrences(allF,C);
-
-  // Q: Why not just count the number of ones?
-  // A: because we are including non-manifold edges as boundary edges
-  int twos = (int) count(C.begin(),C.end(),2);
-  //int ones = (int) count(C.begin(),C.end(),1);
-  // Resize output to fit number of ones
-  F.resize(allF.size() - twos);
-  //F.resize(ones);
-  int k = 0;
-  for(int i = 0;i< (int)allF.size();i++)
-  {
-    if(C[i] != 2)
-    {
-      assert(k<(int)F.size());
-      F[k] = allF[i];
-      k++;
-    }
-  }
-  assert(k==(int)F.size());
-  //if(k != F.size())
-  //{
-  //  printf("%d =? %d\n",k,F.size());
-  //}
-
-}
-
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-// generated by autoexplicit.sh
 template void igl::boundary_facets<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
 template void igl::boundary_facets<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >&);
-// generated by autoexplicit.sh
 template void igl::boundary_facets<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
-template void igl::boundary_facets<int, int>(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);
-//template Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > igl::boundary_facets(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
 template Eigen::Matrix<int, -1, -1, 0, -1, -1> igl::boundary_facets<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
 template void igl::boundary_facets<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&);
 #endif

include/igl/boundary_facets.h

@@ -19,7 +19,9 @@ namespace igl
   /// (analogous to qptoolbox's `outline` and `boundary_faces`).
   ///
   /// @param[in] T  tetrahedron (triangle) index list, m by 4 (3), where m is the number of tetrahedra
-  /// @param[out] F  list of boundary faces, n by 3 (2), where n is the number of boundary faces
+  /// @param[out] F  list of boundary faces, n by 3 (2), where n is the number
+  ///   of boundary faces. Faces are oriented so that igl::centroid(V,F,…)
+  /// computes the same sign volume as igl::volume(V,T)
   /// @param[out] J  list of indices into T, n by 1
   /// @param[out] K  list of indices revealing across from which vertex is this facet
   ///
@@ -48,15 +50,6 @@ namespace igl
   template <typename DerivedT, typename Ret>
   Ret boundary_facets(
     const Eigen::MatrixBase<DerivedT>& T);
-  /// Determine boundary faces (edges) of tetrahedra (triangles) stored in T;
-  /// inputs and outputs lists.
-  ///
-  /// @param[in] T  tetrahedron (triangle) index list, m by 4 (3), where m is the number of tetrahedra
-  /// @param[out] F  list of boundary faces, n by 3 (2), where n is the number of boundary faces
-  template <typename IntegerT, typename IntegerF>
-  IGL_INLINE void boundary_facets(
-    const std::vector<std::vector<IntegerT> > & T,
-    std::vector<std::vector<IntegerF> > & F);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/read_triangle_mesh.cpp

@@ -99,8 +99,6 @@ IGL_INLINE bool igl::read_triangle_mesh(
     if(mF.rows() == 0 && T.rows() > 0)
     {
       boundary_facets(T,F);
-      // outward facing
-      F = F.rowwise().reverse().eval();
     }else
     {
       F = mF.template cast<typename DerivedF::Scalar>();
@@ -144,8 +142,6 @@ IGL_INLINE bool igl::read_triangle_mesh(
     //if(F.size() > T.size() || F.size() == 0)
     {
       boundary_facets(T,F);
-      // outward facing
-      F = F.rowwise().reverse().eval();
     }
   }else if(ext == "obj")
   {

include/igl/vector_area_matrix.cpp

@@ -24,6 +24,7 @@ IGL_INLINE void igl::vector_area_matrix(
   // number of vertices
   const int n = F.maxCoeff()+1;
 
+  assert(F.cols() == 3);
   MatrixXi E;
   boundary_facets(F,E);
 

tests/include/igl/boundary_facets.cpp

@@ -2,11 +2,55 @@
 #include <igl/boundary_facets.h>
 #include <igl/sort.h>
 #include <igl/sortrows.h>
-#include <igl/setxor.h>
+#include <igl/centroid.h>
+#include <igl/volume.h>
 
 #include <igl/matlab_format.h>
 #include <iostream>
 
+TEST_CASE("boundary_facets: single_tet_volume", "[igl]")
+{
+  Eigen::MatrixXd V(4,3);
+  V<<
+    0,0,0,
+    1,0,0,
+    0,1,0,
+    0,0,1;
+  Eigen::MatrixXi T(1,4);
+  T<<0,1,2,3;
+  Eigen::MatrixXi F;
+  igl::boundary_facets(T,F);
+  REQUIRE( F.rows () == 4 );
+  double total_volume;
+  Eigen::RowVector3d centroid;
+  igl::centroid(V,F,centroid,total_volume);
+  Eigen::VectorXd volumes;
+  igl::volume(V,T,volumes);
+  REQUIRE( total_volume == Approx(volumes(0)) );
+}
+
+TEST_CASE("boundary_facets: single_tri", "[igl]")
+{
+  Eigen::MatrixXd V(4,3);
+  V<<
+    0,0,
+    1,0,
+    0,1;
+  Eigen::MatrixXi F(1,3);
+  F<<0,1,2;
+  Eigen::MatrixXi E;
+  igl::boundary_facets(F,E);
+  REQUIRE( E.rows () == 3 );
+  // orientation should match triangle edges
+  Eigen::MatrixXi FE(3,2);
+  FE<<
+    F(0,1),F(0,2),
+    F(0,2),F(0,0),
+    F(0,0),F(0,1);
+  igl::sortrows(Eigen::MatrixXi(E),true,E);
+  igl::sortrows(Eigen::MatrixXi(FE),true,FE);
+  test_common::assert_eq(FE,E);
+}
 
 TEST_CASE("boundary_facets: single_tet", "[igl]")
 {