lscm: fix sign for the area term (#1853)

Added an overload where the matrix Q can be returned (used in the test).

include/igl/lscm.cpp

@@ -18,7 +18,8 @@ IGL_INLINE bool igl::lscm(
   const Eigen::MatrixXi& F,
   const Eigen::VectorXi& b,
   const Eigen::MatrixXd& bc,
-  Eigen::MatrixXd & V_uv)
+  Eigen::MatrixXd & V_uv,
+  Eigen::SparseMatrix<double> & Q)
 {
   using namespace Eigen;
   using namespace std;
@@ -43,7 +44,7 @@ IGL_INLINE bool igl::lscm(
   }
   
   // Minimize the LSCM energy
-  SparseMatrix<double> Q = -L_flat + 2.*A;
+  Q = -L_flat - 2.*A;
   const VectorXd B_flat = VectorXd::Zero(V.rows()*2);
   igl::min_quad_with_fixed_data<double> data;
   if(!igl::min_quad_with_fixed_precompute(Q,b_flat,SparseMatrix<double>(),true,data))
@@ -67,6 +68,17 @@ IGL_INLINE bool igl::lscm(
   return true;
 }
 
+IGL_INLINE bool igl::lscm(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  const Eigen::VectorXi& b,
+  const Eigen::MatrixXd& bc,
+  Eigen::MatrixXd & V_uv)
+{
+  Eigen::SparseMatrix<double> Q;
+  return lscm(V, F, b, bc, V_uv, Q);
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 #endif

include/igl/lscm.h

@@ -32,8 +32,17 @@ namespace igl
   //   bc #b by 3 list of boundary values
   // Outputs:
   //   UV #V by 2 list of 2D mesh vertex positions in UV space
+  //   Q  #Vx2 by #Vx2 symmetric positive semi-definite matrix for computing LSCM energy
   // Returns true only on solver success.
   //
+  IGL_INLINE bool lscm( 
+      const Eigen::MatrixXd& V, 
+      const Eigen::MatrixXi& F,
+      const Eigen::VectorXi& b, 
+      const Eigen::MatrixXd& bc, 
+      Eigen::MatrixXd& V_uv,
+      Eigen::SparseMatrix<double>& Q);
+  // Wrapper where the output Q is discarded
   IGL_INLINE bool lscm( 
       const Eigen::MatrixXd& V, 
       const Eigen::MatrixXi& F,

tests/include/igl/lscm.cpp

@@ -0,0 +1,87 @@
+#include <test_common.h>
+#include <igl/grad.h>
+#include <igl/per_face_normals.h>
+#include <igl/doublearea.h>
+#include <igl/boundary_loop.h>
+#include <igl/lscm.h>
+#include <igl/EPS.h>
+
+double compute_lscm_energy(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  const Eigen::VectorXd& W_flat)
+{
+  const int nV = V.rows();
+  const int nF = F.rows();
+
+  assert(V.cols() == 3);
+  assert(F.cols() == 3);
+  assert(V_uv.rows() == nV);
+  assert(V_uv.cols() == 2);
+
+  // Compute gradient
+  Eigen::SparseMatrix<double> G;
+  igl::grad(V, F, G);
+
+  Eigen::VectorXd nablaU_flat = G * W_flat.head(nV);
+  Eigen::VectorXd nablaV_flat = G * W_flat.tail(nV);
+
+  Eigen::MatrixXd nablaU(nF, 3);
+  Eigen::MatrixXd nablaV(nF, 3);
+  nablaU << nablaU_flat.segment(0, nF), nablaU_flat.segment(nF, nF), nablaU_flat.segment(2 * nF, nF);
+  nablaV << nablaV_flat.segment(0, nF), nablaV_flat.segment(nF, nF), nablaV_flat.segment(2 * nF, nF);
+
+  // Compute face normal
+  Eigen::MatrixXd N;
+  igl::per_face_normals(V, F, N);
+
+  // Compute area
+  Eigen::VectorXd dblA;
+  igl::doublearea(V, F, dblA);
+
+  // Sum up
+  double sum = 0;
+
+  for (int i = 0; i < nF; ++i) {
+    Eigen::Vector3d N_i = N.row(i);
+    Eigen::Vector3d nablaU_i = nablaU.row(i);
+    Eigen::Vector3d nablaV_i = nablaV.row(i);
+
+    // Rotate nablaU_i about N_i by 90 degrees ccw
+    nablaU_i = N_i.cross(nablaU_i);
+
+    sum += 0.5 * dblA[i] * (nablaU_i - nablaV_i).squaredNorm();
+  }
+
+  return 0.5 * sum;
+}
+
+TEST_CASE("lscm: lscm_energy_check", "[igl]")
+{
+  // Load a mesh in OFF format
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  igl::read_triangle_mesh(test_common::data_path("camelhead.off"), V, F);
+
+  // Fix two points on the boundary
+  Eigen::VectorXi bnd, b(2, 1);
+  igl::boundary_loop(F, bnd);
+  b(0) = bnd(0);
+  b(1) = bnd(bnd.size() / 2);
+  Eigen::MatrixXd bc(2, 2);
+  bc << 0, 0, 1, 0;
+
+  // LSCM parametrization
+  Eigen::MatrixXd V_uv;
+  Eigen::SparseMatrix<double> Q;
+  igl::lscm(V, F, b, bc, V_uv, Q);
+
+  Eigen::VectorXd W_flat(2 * V.rows());
+  W_flat << V_uv.col(0), V_uv.col(1);
+
+  // Consistency check
+  double e1 = compute_lscm_energy(V, F, W_flat);
+  double e2 = 0.5 * W_flat.transpose() * Q * W_flat;
+
+  REQUIRE(std::abs(e1 - e2) < igl::EPS<double>());
+}