clean up lscm and areamatrix

Former-commit-id: 41385b95723ed4bc969c25afdbbc72ae9d35b318

include/igl/areamatrix.cpp

@@ -12,34 +12,37 @@
 #include <iostream>
 #include <unsupported/Eigen/SparseExtra>
 
-#include <igl/boundary_vertices_sorted.h>
+//#include <igl/boundary_vertices_sorted.h>
+#include <igl/boundary_faces.h>
 
-template <typename DerivedV, typename DerivedF, typename Scalar>
+template <typename DerivedF, typename Scalar>
 IGL_INLINE void igl::areamatrix(
-  const Eigen::PlainObjectBase<DerivedV> & V,
   const Eigen::PlainObjectBase<DerivedF> & F,
   Eigen::SparseMatrix<Scalar>& A)
 {
   using namespace Eigen;
   using namespace std;
 
-	SparseMatrix<Scalar> aux (V.rows() * 2, V.rows() * 2);
-	SparseMatrix<Scalar> auxT(V.rows() * 2, V.rows() * 2);
+  // number of vertices
+  const int n = F.maxCoeff()+1;
+
+	SparseMatrix<Scalar> aux (n * 2, n * 2);
+	SparseMatrix<Scalar> auxT(n * 2, n * 2);
 
 	vector<Triplet<Scalar> > auxTripletList;
 	vector<Triplet<Scalar> > auxTTripletList;
 
-  Eigen::VectorXi bnd;
-  igl::boundary_vertices_sorted(V,F,bnd);
+  MatrixXi E;
+  boundary_faces(F,E);
 
-	for(int k = 0; k < bnd.size(); k++)
+	for(int k = 0; k < E.rows(); k++)
   {
-		int i = bnd[k];
-		int j = k + 1 == bnd.size() ? bnd[0] : bnd[k+1];
-		auxTripletList.push_back(Triplet<Scalar>(i+V.rows(), j, -0.5));
-		auxTripletList.push_back(Triplet<Scalar>(i, j+V.rows(), 0.5));
-		auxTTripletList.push_back(Triplet<Scalar>(j, i+V.rows(), -0.5));
-		auxTTripletList.push_back(Triplet<Scalar>(j+V.rows(), i, 0.5));
+		int i = E(k,0);
+		int j = E(k,1);
+		auxTripletList.push_back(Triplet<Scalar>(i+n, j, -0.5));
+		auxTripletList.push_back(Triplet<Scalar>(i, j+n, 0.5));
+		auxTTripletList.push_back(Triplet<Scalar>(j, i+n, -0.5));
+		auxTTripletList.push_back(Triplet<Scalar>(j+n, i, 0.5));
 	}
 
 	aux.setFromTriplets(auxTripletList.begin(), auxTripletList.end());
@@ -49,5 +52,5 @@ IGL_INLINE void igl::areamatrix(
 
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
-// generated by autoexplicit.sh
+template void igl::areamatrix<Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
 #endif

include/igl/areamatrix.h

@@ -14,9 +14,8 @@
 
 namespace igl
 {
-  // Constructs the symmetric area matrix A, s.t.
-  // [V.col(0)' V.col(1)'] * A * [V.col(0); V.col(1)] is the area of the planar mesh (V,F).
-  // Note: (V,F) must be a genus-0 mesh, with a single boundary
+  // Constructs the symmetric area matrix A, s.t.  [V.col(0)' V.col(1)'] * A *
+  // [V.col(0); V.col(1)] is the **vector area** of the mesh (V,F).
   //
   // Templates:
   //   DerivedV  derived type of eigen matrix for V (e.g. derived from
@@ -25,14 +24,12 @@ namespace igl
   //     MatrixXi)
   //   Scalar  scalar type for eigen sparse matrix (e.g. double)
   // Inputs:
-  //   V  #V by 2 list of mesh vertex positions
   //   F  #F by 3 list of mesh faces (must be triangles)
   // Outputs:
   //   A  #Vx2 by #Vx2 area matrix
   //
-  template <typename DerivedV, typename DerivedF, typename Scalar>
+  template <typename DerivedF, typename Scalar>
   IGL_INLINE void areamatrix(
-    const Eigen::PlainObjectBase<DerivedV> & V,
     const Eigen::PlainObjectBase<DerivedF> & F,
     Eigen::SparseMatrix<Scalar>& A);
 }

include/igl/kronecker_product.cpp

@@ -61,4 +61,5 @@ IGL_INLINE Eigen::SparseMatrix<Scalar> igl::kronecker_product(
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
 // generated by autoexplicit.sh
+template Eigen::SparseMatrix<double, 0, int> igl::kronecker_product<double>(Eigen::SparseMatrix<double, 0, int> const&, Eigen::SparseMatrix<double, 0, int> const&);
 #endif

include/igl/kronecker_product.h

@@ -19,8 +19,7 @@ namespace igl
   // Inputs:
   //   A  #M by #N sparse matrix
   //   B  #P by #Q sparse matrix
-  // Outputs:
-  //      #M*#P by #N*#Q sparse matrix
+  // Returns #M*#P by #N*#Q sparse matrix
   //
   template <typename Scalar>
   IGL_INLINE Eigen::SparseMatrix<Scalar> kronecker_product(

include/igl/lscm.cpp

@@ -7,56 +7,43 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "lscm.h"
 
-// Bug in unsupported/Eigen/SparseExtra needs iostream first
-#include <iostream>
-#include <unsupported/Eigen/SparseExtra>
-
 #include <igl/areamatrix.h>
 #include <igl/cotmatrix.h>
-#include <igl/kronecker_product.h>
-#include <igl/boundary_vertices_sorted.h>
+//#include <igl/kronecker_product.h>
+#include <igl/repdiag.h>
 #include <igl/min_quad_with_fixed.h>
+#include <igl/matlab_format.h>
+#include <iostream>
 
-IGL_INLINE Eigen::MatrixXd igl::lscm(
+IGL_INLINE void igl::lscm(
   const Eigen::MatrixXd& V,
   const Eigen::MatrixXi& F,
   const Eigen::VectorXi& b,
-  const Eigen::MatrixXd& bc)
+  const Eigen::MatrixXd& bc,
+  Eigen::MatrixXd & V_uv)
 {
   using namespace Eigen;
+  using namespace std;
   
   // Assemble the area matrix (note that A is #Vx2 by #Vx2)
   SparseMatrix<double> A;
-  igl::areamatrix(V,F,A);
+  igl::areamatrix(F,A);
 
   // Assemble the cotan laplacian matrix
   SparseMatrix<double> L;
   igl::cotmatrix(V,F,L);
 
-  SparseMatrix<double> I2(2,2);
-  I2.insert(0,0) = 1;
-  I2.insert(1,1) = 1;
+  SparseMatrix<double> L_flat;
+  repdiag(L,2,L_flat);
 
-  SparseMatrix<double> L_flat = igl::kronecker_product(I2,L);
-
-  VectorXi b_flat = b;
-  MatrixXd bc_flat = bc;
-  
-  if (b.size() < 2)
+  VectorXi b_flat(b.size()*bc.cols(),1);
+  VectorXd bc_flat(bc.size(),1);
+  for(int c = 0;c<bc.cols();c++)
   {
-    // if no boundary conditions are provided, fix two boundary points
-    Eigen::VectorXi bnd;
-    igl::boundary_vertices_sorted(V,F,bnd);
-
-    int v1 = bnd(0);
-    int v2 = bnd(round(bnd.size()/2));
-
-    b_flat.resize(4);
-    b_flat << v1, v1+1, v2, v2+1;
-    bc_flat.resize(4,1);
-    bc_flat << 0, 1, 0, 0;
+    b_flat.block(c*b.size(),0,b.rows(),1) = c*V.rows() + b.array();
+    bc_flat.block(c*bc.rows(),0,bc.rows(),1) = bc.col(c);
   }
-
+  
   // Minimize the LSCM energy
   SparseMatrix<double> Q = -0.5*L_flat + A;
   const VectorXd B_flat = VectorXd::Zero(V.rows()*2);
@@ -67,17 +54,16 @@ IGL_INLINE Eigen::MatrixXd igl::lscm(
   if(!min_quad_with_fixed_solve(data,B_flat,bc_flat,VectorXd(),W_flat))
     assert(false);
 
-  MatrixXd V_uv(V.rows(),2);
-  for (unsigned i=0;i<V.rows();++i)
+
+  assert(W_flat.rows() == V.rows()*2);
+  V_uv.resize(V.rows(),2);
+  for (unsigned i=0;i<V_uv.cols();++i)
   {
-    V_uv(i,0) = W_flat(i);
-    V_uv(i,1) = W_flat(i+V.rows());
+    V_uv.col(i) = W_flat.block(V_uv.rows()*i,0,V_uv.rows(),1);
   }
 
-  return V_uv;
 }
 
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
-// generated by autoexplicit.sh
 #endif

include/igl/lscm.h

@@ -17,24 +17,28 @@ namespace igl
   // Compute a Least-squares conformal map parametrization following the algorithm
   // presented in: Spectral Conformal Parameterization,
   //               Patrick Mullen, Yiying Tong, Pierre Alliez and Mathieu Desbrun
-  // Note: (V,F) must be a genus-0 mesh, with a single boundary
   //
   // Inputs:
   //   V  #V by 3 list of mesh vertex positions
   //   F  #F by 3 list of mesh faces (must be triangles)
   //   b  #b boundary indices into V
-  //   bc #b by #W list of boundary values
-// Outputs:
+  //   bc #b by 3 list of boundary values
+  // Outputs:
   //   UV #V by 2 list of 2D mesh vertex positions in UV space
   //
-  // Note: if b and bc are empty, lscm automatically removes the null space
-  //       by fixing two farthest points on the boundary
-
-  IGL_INLINE Eigen::MatrixXd lscm(
+  // X  Note: if b and bc are empty, lscm automatically removes the null space
+  // X        by fixing two farthest points on the boundary
+  // Boundary conditions are part of the api. It would be strange to secretly
+  // use other boundary conditions. This is also weird if there's more than one
+  // boundary loop.
+  // X  Note: (V,F) must be a genus-0 mesh, with a single boundary
+  // No longer the case. Should probably just be manifold.
+  IGL_INLINE void lscm(
     const Eigen::MatrixXd& V,
     const Eigen::MatrixXi& F,
     const Eigen::VectorXi& b,
-    const Eigen::MatrixXd& bc);
+    const Eigen::MatrixXd& bc,
+    Eigen::MatrixXd& V_uv);
 }
 
 #ifdef IGL_HEADER_ONLY

include/igl/min_quad_with_fixed.cpp

@@ -535,5 +535,6 @@ IGL_INLINE bool igl::min_quad_with_fixed_solve(
 template bool igl::min_quad_with_fixed_solve<double, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(igl::min_quad_with_fixed_data<double> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template bool igl::min_quad_with_fixed_precompute<double, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<double, 0, int> const&, bool, igl::min_quad_with_fixed_data<double>&);
 template bool igl::min_quad_with_fixed_solve<double, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(igl::min_quad_with_fixed_data<double> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template bool igl::min_quad_with_fixed_solve<double, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(igl::min_quad_with_fixed_data<double> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif
 

tutorial/502_LSCMParam/main.cpp

@@ -1,5 +1,6 @@
 #include <igl/readOFF.h>
 #include <igl/viewer/Viewer.h>
+#include <igl/boundary_vertices_sorted.h>
 
 #include <igl/lscm.h>
 
@@ -30,8 +31,16 @@ int main(int argc, char *argv[])
   // Load a mesh in OFF format
   igl::readOFF("../shared/camelhead.off", V, F);
 
+  // Fix two points on the boundary
+  VectorXi bnd,b(2,1);
+  igl::boundary_vertices_sorted(V,F,bnd);
+  b(0) = bnd(0);
+  b(1) = bnd(round(bnd.size()/2));
+  MatrixXd bc(2,2);
+  bc<<0,0,1,0;
+
   // LSCM parametrization
-  V_uv = igl::lscm(V,F,Eigen::VectorXi(),Eigen::MatrixXd());
+  igl::lscm(V,F,b,bc,V_uv);
 
   // Scale the uv
   V_uv *= 5;