Adjust igl::boundary_conditions(...) to support 3D cages (#2229)

* Calculate boundary conditions for cage faces

* Remove std:: prefix to conform to original libigl

* Adjust tests to new boundary_conditions(...)

* Adjusting tutorial 403 to new boundary_conditions(...)

* Removing unnecessary include

* Remove tolerances from boundary conditions for cages to prevent dups

---------

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

include/igl/boundary_conditions.cpp

@@ -22,6 +22,7 @@ IGL_INLINE bool igl::boundary_conditions(
   const Eigen::VectorXi & P  ,
   const Eigen::MatrixXi & BE ,
   const Eigen::MatrixXi & CE ,
+  const Eigen::MatrixXi & CF ,
   Eigen::VectorXi &       b  ,
   Eigen::MatrixXd &       bc )
 {
@@ -121,6 +122,58 @@ IGL_INLINE bool igl::boundary_conditions(
     }
   }
 
+  std::vector<uint8_t> vertices_marked(V.rows(), 0);
+  // loop over cage faces
+  for(int f = 0;f<CF.rows();f++)
+  {
+    Vector3d v_0 = C.row(P(CF(f, 0)));
+    Vector3d v_1 = C.row(P(CF(f, 1)));
+    Vector3d v_2 = C.row(P(CF(f, 2)));
+    Vector3d n = (v_1 - v_0).cross(v_2 - v_1);
+    n.normalize();
+    // loop over domain vertices
+    for (int i = 0;i<V.rows();i++)
+    {
+      // ensure each vertex is associated with only one face
+      if (vertices_marked[i])
+      {
+          continue;
+      }
+      Vector3d point = V.row(i);
+      Vector3d v = point - v_0;
+      double dist = abs(v.dot(n));
+      Vector3d projected_point = point - dist * n;
+      if (dist <= 1.e-1f)
+      {
+        //barycentric coordinates
+        Vector3d vec_0 = v_1 - v_0, vec_1 = v_2 - v_0, vec_2 = point - v_0;
+        double d00 = vec_0.dot(vec_0);
+        double d01 = vec_0.dot(vec_1);
+        double d11 = vec_1.dot(vec_1);
+        double d20 = vec_2.dot(vec_0);
+        double d21 = vec_2.dot(vec_1);
+        double denom = d00 * d11 - d01 * d01;
+        double v = (d11 * d20 - d01 * d21) / denom;
+        double w = (d00 * d21 - d01 * d20) / denom;
+        double u = 1.0 - v - w;
+
+        if (u>=0. && u<=1.0 && v>=0. && v<=1.0 && w >=0. && w<=1.0)
+        {
+          vertices_marked[i] = 1;
+          bci.push_back(i);
+          bcj.push_back(CF(f, 0));
+          bcv.push_back(u);
+          bci.push_back(i);
+          bcj.push_back(CF(f, 1));
+          bcv.push_back(v);
+          bci.push_back(i);
+          bcj.push_back(CF(f, 2));
+          bcv.push_back(w);
+        }
+      }
+    }
+  }
+
   // find unique boundary indices
   vector<int> vb = bci;
   sort(vb.begin(),vb.end());

include/igl/boundary_conditions.h

@@ -15,7 +15,7 @@ namespace igl
   /// Compute boundary conditions for automatic weights computation. This
   /// function expects that the given mesh (V,Ele) has sufficient samples
   /// (vertices) exactly at point handle locations and exactly along bone and
-  /// cage edges.
+  /// cage edges/faces.
   ///
   /// @param[in] V  #V by dim list of domain vertices
   /// @param[in] Ele  #Ele by simplex-size list of simplex indices
@@ -23,6 +23,7 @@ namespace igl
   /// @param[in] P  #P by 1 list of point handle indices into C
   /// @param[in] BE  #BE by 2 list of bone edge indices into C
   /// @param[in] CE  #CE by 2 list of cage edge indices into *P*
+  ///  @param[in] CF  #CF by 3 list of (triangular) cage face indices into *P*
   /// @param[out] b  #b list of boundary indices (indices into V of vertices which have
   ///     known, fixed values)
   /// @param[out] bc #b by #weights list of known/fixed values for boundary vertices
@@ -35,7 +36,6 @@ namespace igl
   ///   some column of bc doesn't have a 0 (assuming bc has >1 columns)
   ///   some column of bc doesn't have a 1 (assuming bc has >1 columns)
   ///
-  /// \note 3D cages are not yet supported.
   IGL_INLINE bool boundary_conditions(
     const Eigen::MatrixXd & V,
     const Eigen::MatrixXi & Ele,
@@ -43,6 +43,7 @@ namespace igl
     const Eigen::VectorXi & P,
     const Eigen::MatrixXi & BE,
     const Eigen::MatrixXi & CE,
+    const Eigen::MatrixXi & CF,
     Eigen::VectorXi & b,
     Eigen::MatrixXd & bc);
 }

tests/include/igl/bbw.cpp

@@ -16,7 +16,7 @@ TEST_CASE("bbw: decimated_knight", "[igl]" "[slow]")
     test_common::data_path("decimated-knight-matlab-active-set.dmat"),W_groundtruth);
   Eigen::VectorXi b;
   Eigen::MatrixXd bc;
-  igl::boundary_conditions(V,T,C,Eigen::VectorXi(),E,Eigen::MatrixXi(),b,bc);
+  igl::boundary_conditions(V,T,C,Eigen::VectorXi(),E,Eigen::MatrixXi(),Eigen::MatrixXi(),b,bc);
   igl::BBWData params;
   params.active_set_params.max_iter = 100;
   igl::bbw(V,T,b,bc,params,Was);

tests/include/igl/mosek/bbw.cpp

@@ -16,7 +16,7 @@ TEST_CASE("mosek_bbw: decimated_knight", "[igl/copyleft/mosek]")
     test_common::data_path("decimated-knight-matlab-active-set.dmat"),W_groundtruth);
   Eigen::VectorXi b;
   Eigen::MatrixXd bc;
-  igl::boundary_conditions(V,T,C,Eigen::VectorXi(),E,Eigen::MatrixXi(),b,bc);
+  igl::boundary_conditions(V,T,C,Eigen::VectorXi(),E,Eigen::MatrixXi(),Eigen::MatrixXi(),b,bc);
   igl::BBWData params;
   igl::mosek::MosekData mosek_params;
   igl::mosek::bbw(V,T,b,bc,params,mosek_params,Wmo);

tutorial/403_BoundedBiharmonicWeights/main.cpp

@@ -116,7 +116,7 @@ int main(int argc, char *argv[])
   VectorXi b;
   // List of boundary conditions of each weight function
   MatrixXd bc;
-  igl::boundary_conditions(V,T,C,VectorXi(),BE,MatrixXi(),b,bc);
+  igl::boundary_conditions(V,T,C,VectorXi(),BE,MatrixXi(),MatrixXi(),b,bc);
 
   // compute BBW weights matrix
   igl::BBWData bbw_data;