SCAF: Expose linear system to be solved by user code (#1553)

* Add const qualifiers

* Add scaf_system() and expose scaf::compute_energy()

* Add test "scaf_system"

Co-authored-by: Patrick Schmidt <[email protected]>

include/igl/flip_avoiding_line_search.cpp

@@ -303,7 +303,7 @@ namespace igl
 IGL_INLINE double igl::flip_avoiding_line_search(
   const Eigen::MatrixXi F,
   Eigen::MatrixXd& cur_v,
-  Eigen::MatrixXd& dst_v,
+  const Eigen::MatrixXd& dst_v,
   std::function<double(Eigen::MatrixXd&)> energy,
   double cur_energy)
 {

include/igl/flip_avoiding_line_search.h

@@ -36,7 +36,7 @@ namespace igl
   IGL_INLINE double flip_avoiding_line_search(
     const Eigen::MatrixXi F,
     Eigen::MatrixXd& cur_v,
-    Eigen::MatrixXd& dst_v,
+    const Eigen::MatrixXd& dst_v,
     std::function<double(Eigen::MatrixXd&)> energy,
     double cur_energy = -1);
 

include/igl/scaf.cpp

@@ -332,7 +332,7 @@ IGL_INLINE double compute_soft_constraint_energy(const SCAFData &s)
   return e;
 }
 
-IGL_INLINE double compute_energy(SCAFData &s, Eigen::MatrixXd &w_uv, bool whole)
+IGL_INLINE double compute_energy(SCAFData &s, const Eigen::MatrixXd &w_uv, bool whole)
 {
   if (w_uv.rows() != s.v_num)
     assert(!whole);
@@ -529,6 +529,30 @@ IGL_INLINE void build_scaffold_linear_system(const SCAFData &s, Eigen::SparseMat
   rhs = Aut * (frhs - Ae * known_pos);
 }
 
+IGL_INLINE void build_weighted_arap_system(SCAFData &s, Eigen::SparseMatrix<double> &L, Eigen::VectorXd &rhs)
+{
+  // fixed frame solving:
+  // x_e as the fixed frame, x_u for unknowns (mesh + unknown scaffold)
+  // min ||(A_u*x_u + A_e*x_e) - b||^2
+  // => A_u'*A_u*x_u  = Au'* (b - A_e*x_e) := Au'* b_u
+  //
+  // separate matrix build:
+  // min ||A_m x_m - b_m||^2 + ||A_s x_all - b_s||^2 + soft + proximal
+  // First change dimension of A_m to fit for x_all
+  // (Not just at the end, since x_all is flattened along dimensions)
+  // L = A_m'*A_m + A_s'*A_s + soft + proximal
+  // rhs = A_m'* b_m + A_s' * b_s + soft + proximal
+  //
+  Eigen::SparseMatrix<double> L_m, L_s;
+  Eigen::VectorXd rhs_m, rhs_s;
+  build_surface_linear_system(s, L_m, rhs_m);  // complete Am, with soft
+  build_scaffold_linear_system(s, L_s, rhs_s); // complete As, without proximal
+
+  L = L_m + L_s;
+  rhs = rhs_m + rhs_s;
+  L.makeCompressed();
+}
+
 IGL_INLINE void solve_weighted_arap(SCAFData &s, Eigen::MatrixXd &uv)
 {
   using namespace Eigen;
@@ -562,27 +586,7 @@ IGL_INLINE void solve_weighted_arap(SCAFData &s, Eigen::MatrixXd &uv)
 
   Eigen::SparseMatrix<double> L;
   Eigen::VectorXd rhs;
-
-  // fixed frame solving:
-  // x_e as the fixed frame, x_u for unknowns (mesh + unknown scaffold)
-  // min ||(A_u*x_u + A_e*x_e) - b||^2
-  // => A_u'*A_u*x_u  = Au'* (b - A_e*x_e) := Au'* b_u
-  //
-  // separate matrix build:
-  // min ||A_m x_m - b_m||^2 + ||A_s x_all - b_s||^2 + soft + proximal
-  // First change dimension of A_m to fit for x_all
-  // (Not just at the end, since x_all is flattened along dimensions)
-  // L = A_m'*A_m + A_s'*A_s + soft + proximal
-  // rhs = A_m'* b_m + A_s' * b_s + soft + proximal
-  //
-  Eigen::SparseMatrix<double> L_m, L_s;
-  Eigen::VectorXd rhs_m, rhs_s;
-  build_surface_linear_system(s, L_m, rhs_m);  // complete Am, with soft
-  build_scaffold_linear_system(s, L_s, rhs_s); // complete As, without proximal
-
-  L = L_m + L_s;
-  rhs = rhs_m + rhs_s;
-  L.makeCompressed();
+  build_weighted_arap_system(s, L, rhs);
 
   Eigen::VectorXd unknown_Uc((v_n - s.frame_ids.size() - s.fixed_ids.size()) * dim), Uc(dim * v_n);
 
@@ -612,6 +616,7 @@ IGL_INLINE double perform_iteration(SCAFData &s)
                                         whole_E, -1) /
          s.mesh_measure;
 }
+
 }
 }
 
@@ -702,5 +707,24 @@ IGL_INLINE Eigen::MatrixXd igl::scaf_solve(SCAFData &s, int iter_num)
   return s.w_uv.topRows(s.mv_num);
 }
 
+IGL_INLINE void igl::scaf_system(SCAFData &s, Eigen::SparseMatrix<double> &L, Eigen::VectorXd &rhs)
+{
+    s.energy = igl::scaf::compute_energy(s, s.w_uv, false) / s.mesh_measure;
+
+    s.total_energy = igl::scaf::compute_energy(s, s.w_uv, true) / s.mesh_measure;
+    s.rect_frame_V = Eigen::MatrixXd();
+    igl::scaf::mesh_improve(s);
+
+    double new_weight = s.mesh_measure * s.energy / (s.sf_num * 100);
+    s.scaffold_factor = new_weight;
+    igl::scaf::update_scaffold(s);
+
+    igl::scaf::compute_jacobians(s, s.w_uv, true);
+    igl::slim_update_weights_and_closest_rotations_with_jacobians(s.Ji_m, s.slim_energy, 0, s.W_m, s.Ri_m);
+    igl::slim_update_weights_and_closest_rotations_with_jacobians(s.Ji_s, s.scaf_energy, 0, s.W_s, s.Ri_s);
+
+    igl::scaf::build_weighted_arap_system(s, L, rhs);
+}
+
 #ifdef IGL_STATIC_LIBRARY
 #endif

include/igl/scaf.h

@@ -68,15 +68,15 @@ namespace igl
     };
 
 
-// Compute necessary information to start using SCAF
-// Inputs:
-//		V           #V by 3 list of mesh vertex positions
-//		F           #F by 3/3 list of mesh faces (triangles/tets)
-//    data          igl::SCAFData
-//    slim_energy Energy type to minimize
-//    b           list of boundary indices into V (soft constraint)
-//    bc          #b by dim list of boundary conditions (soft constraint)
-//    soft_p      Soft penalty factor (can be zero)
+    // Compute necessary information to start using SCAF
+    // Inputs:
+    //		V           #V by 3 list of mesh vertex positions
+    //		F           #F by 3/3 list of mesh faces (triangles/tets)
+    //    data          igl::SCAFData
+    //    slim_energy Energy type to minimize
+    //    b           list of boundary indices into V (soft constraint)
+    //    bc          #b by dim list of boundary conditions (soft constraint)
+    //    soft_p      Soft penalty factor (can be zero)
     IGL_INLINE void scaf_precompute(
         const Eigen::MatrixXd &V,
         const Eigen::MatrixXi &F,
@@ -87,12 +87,30 @@ namespace igl
         Eigen::MatrixXd& bc,
         double soft_p);
 
+    // Run iter_num iterations of SCAF, with precomputed data
+    // Outputs:
+    //    V_o (in SLIMData): #V by dim list of mesh vertex positions
+    IGL_INLINE Eigen::MatrixXd scaf_solve(SCAFData &data, int iter_num);
 
-// Run iter_num iterations of SCAF, with precomputed data
-// Outputs:
-//    V_o (in SLIMData): #V by dim list of mesh vertex positions
-  IGL_INLINE Eigen::MatrixXd scaf_solve(SCAFData &data, int iter_num);
-  }
+    // Set up the SCAF system L * uv = rhs, without solving it.
+    // Inputs:
+    //    s:   igl::SCAFData. Will be modified by energy and Jacobian computation.
+    // Outputs:
+    //    L:   m by m matrix
+    //    rhs: m by 1 vector
+    //         with m = dim * (#V_mesh + #V_scaf - #V_frame)
+    IGL_INLINE void scaf_system(SCAFData &s, Eigen::SparseMatrix<double> &L, Eigen::VectorXd &rhs);
+
+    namespace scaf
+    {
+        // Compute SCAF energy
+        // Inputs:
+        //    s:     igl::SCAFData
+        //    w_uv:  (#V_mesh + #V_scaf) by dim matrix
+        //    whole: Include scaffold if true
+        IGL_INLINE double compute_energy(SCAFData &s, const Eigen::MatrixXd &w_uv, bool whole);
+    }
+}
 
 #ifndef IGL_STATIC_LIBRARY
 #  include "scaf.cpp"

include/igl/slim.cpp

@@ -51,11 +51,11 @@ namespace igl
     // Definitions of internal functions
     IGL_INLINE void buildRhs(igl::SLIMData& s, const Eigen::SparseMatrix<double> &A);
     IGL_INLINE void add_soft_constraints(igl::SLIMData& s, Eigen::SparseMatrix<double> &L);
-    IGL_INLINE double compute_energy(igl::SLIMData& s, Eigen::MatrixXd &V_new);
+    IGL_INLINE double compute_energy(igl::SLIMData& s, const Eigen::MatrixXd &V_new);
     IGL_INLINE double compute_soft_const_energy(igl::SLIMData& s,
                                                 const Eigen::MatrixXd &V,
                                                 const Eigen::MatrixXi &F,
-                                                Eigen::MatrixXd &V_o);
+                                                const Eigen::MatrixXd &V_o);
 
     IGL_INLINE void solve_weighted_arap(igl::SLIMData& s,
                                         const Eigen::MatrixXd &V,
@@ -287,7 +287,7 @@ namespace igl
       }
     }
 
-    IGL_INLINE double compute_energy(igl::SLIMData& s, Eigen::MatrixXd &V_new)
+    IGL_INLINE double compute_energy(igl::SLIMData& s, const Eigen::MatrixXd &V_new)
     {
       compute_jacobians(s,V_new);
       return mapping_energy_with_jacobians(s.Ji, s.M, s.slim_energy, s.exp_factor) +
@@ -297,7 +297,7 @@ namespace igl
     IGL_INLINE double compute_soft_const_energy(igl::SLIMData& s,
                                                 const Eigen::MatrixXd &V,
                                                 const Eigen::MatrixXi &F,
-                                                Eigen::MatrixXd &V_o)
+                                                const Eigen::MatrixXd &V_o)
     {
       double e = 0;
       for (int i = 0; i < s.b.rows(); i++)

tests/include/igl/scaf.cpp

@@ -0,0 +1,80 @@
+#include <test_common.h>
+#include <igl/cat.h>
+#include <igl/scaf.h>
+#include <igl/slice.h>
+#include <igl/harmonic.h>
+#include <igl/boundary_loop.h>
+#include <igl/map_vertices_to_circle.h>
+#include <igl/flip_avoiding_line_search.h>
+
+// Assert that building the SCAF equation system via scaf_system() and
+// solving it manually yields the same result as calling scaf_solve() directly.
+TEST_CASE("scaf_system: Test scaf_system() vs scaf_solve()", "[igl]")
+{
+  // Read cube mesh
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  igl::read_triangle_mesh(test_common::data_path("cube.obj"), V, F);
+
+  // Cut to disk
+  F = F.topRows(F.rows() - 1).eval();
+
+  // Compute Tutte's embedding
+  Eigen::MatrixXd uv_init;
+  {
+      Eigen::MatrixXd bnd_uv;
+      std::vector<int> boundary;
+      igl::boundary_loop(F, boundary);
+      Eigen::VectorXi bnd = Eigen::Map<Eigen::VectorXi>(boundary.data(), boundary.size());
+      igl::map_vertices_to_circle(V, bnd, bnd_uv);
+      igl::harmonic(F, bnd, bnd_uv ,1, uv_init);
+      uv_init.conservativeResize(V.rows(), 2);
+  }
+
+  // Init empty constraint vectors
+  Eigen::VectorXi b;
+  Eigen::MatrixXd bc;
+
+  // Run one scaf iteration as reference
+  igl::SCAFData s_ref;
+  {
+      igl::scaf_precompute(V, F, uv_init, s_ref, igl::MappingEnergyType::SYMMETRIC_DIRICHLET, b, bc, 0);
+      igl::scaf_solve(s_ref, 1);
+  }
+
+  // Obtain SCAF linear system perform iteration manually
+  igl::SCAFData s_test;
+  {
+      // Set up system
+      igl::scaf_precompute(V, F, uv_init, s_test, igl::MappingEnergyType::SYMMETRIC_DIRICHLET, b, bc, 0);
+      Eigen::SparseMatrix<double> L;
+      Eigen::VectorXd rhs;
+      igl::scaf_system(s_test, L, rhs);
+
+      // Solve
+      Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
+      auto x = solver.compute(L).solve(rhs);
+
+      // Write result to uv_target
+      Eigen::MatrixXd uv_target = s_test.w_uv;
+      const int n_v_var = s_test.w_uv.rows() - s_test.frame_ids.size();
+      for (int i = 0; i < n_v_var; ++i)
+      {
+          int row = i;
+          if (i >= s_test.v_num)
+              row += s_test.frame_ids.size();
+
+          uv_target(row, 0) = x[i];
+          uv_target(row, 1) = x[i + n_v_var];
+      }
+
+      // Line search
+      auto E = [&s_test](Eigen::MatrixXd &uv) { return igl::scaf::compute_energy(s_test, uv, true); };
+      Eigen::MatrixXi w_T;
+      igl::cat(1, s_test.m_T, s_test.s_T, w_T);
+      igl::flip_avoiding_line_search(w_T, s_test.w_uv, uv_target, E, -1);
+  }
+
+  // Compare both results
+  REQUIRE(s_test.w_uv == s_ref.w_uv);
+}