rewrite the uniform laplacian as overloading function and restore the previous tutorial

Former-commit-id: 65b5d854cad7f839089539a3e424602dd4b3d66e

include/igl/harmonic.cpp

@@ -26,47 +26,32 @@ IGL_INLINE bool igl::harmonic(
   const Eigen::PlainObjectBase<DerivedF> & F,
   const Eigen::PlainObjectBase<Derivedb> & b,
   const Eigen::PlainObjectBase<Derivedbc> & bc,
-  const int k,
-  const bool use_cotan_weights,
+  int k,
   Eigen::PlainObjectBase<DerivedW> & W)
 {
   using namespace Eigen;
   typedef typename DerivedV::Scalar Scalar;
   typedef Matrix<Scalar,Dynamic,1> VectorXS;
   SparseMatrix<Scalar> Q;
-  if(use_cotan_weights == true){
-    SparseMatrix<Scalar> L,M,Mi;
-    cotmatrix(V,F,L);
-    switch(F.cols())
-    {
-      case 3:
-        massmatrix(V,F,MASSMATRIX_TYPE_VORONOI,M);
-        break;
-      case 4:
-      default:
-        massmatrix(V,F,MASSMATRIX_TYPE_BARYCENTRIC,M);
+  SparseMatrix<Scalar> L,M,Mi;
+  cotmatrix(V,F,L);
+  switch(F.cols())
+  {
+    case 3:
+      massmatrix(V,F,MASSMATRIX_TYPE_VORONOI,M);
       break;
-    }
-    invert_diag(M,Mi);
-    Q = -L;
-    for(int p = 1;p<k;p++)
-    {
-      Q = (Q*Mi*-L).eval();
-    }
+    case 4:
+    default:
+      massmatrix(V,F,MASSMATRIX_TYPE_BARYCENTRIC,M);
+    break;
   }
-  else
+  invert_diag(M,Mi);
+  Q = -L;
+  for(int p = 1;p<k;p++)
   {
-    SparseMatrix<Scalar> A;
-    adjacency_matrix(F,A);
-    // sum each row
-    SparseVector<Scalar> Asum;
-    sum(A,1,Asum);
-    // Convert row sums into diagonal of sparse matrix
-    SparseMatrix<Scalar> Adiag;
-    diag(Asum,Adiag);
-    // Build uniform laplacian
-    Q = -A+Adiag;
+    Q = (Q*Mi*-L).eval();
   }
+
   min_quad_with_fixed_data<Scalar> data;
   min_quad_with_fixed_precompute(Q,b,SparseMatrix<Scalar>(),true,data);
   W.resize(V.rows(),bc.cols());
@@ -84,9 +69,57 @@ IGL_INLINE bool igl::harmonic(
   return true;
 }
 
+template <
+  typename DerivedF,
+  typename Derivedb,
+  typename Derivedbc,
+  typename DerivedW>
+IGL_INLINE bool igl::harmonic(
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  const Eigen::PlainObjectBase<Derivedb> & b,
+  const Eigen::PlainObjectBase<Derivedbc> & bc,
+  int k,
+  Eigen::PlainObjectBase<DerivedW> & W)
+{
+  using namespace Eigen;
+  typedef typename Derivedbc::Scalar Scalar;
+  typedef Matrix<Scalar,Dynamic,1> VectorXS;
+  SparseMatrix<Scalar> Q;
+  SparseMatrix<Scalar> A;
+  adjacency_matrix(F,A);
+  // sum each row
+  SparseVector<Scalar> Asum;
+  sum(A,1,Asum);
+  // Convert row sums into diagonal of sparse matrix
+  SparseMatrix<Scalar> Adiag;
+  diag(Asum,Adiag);
+  // Build uniform laplacian
+  Q = -A+Adiag;
+  min_quad_with_fixed_data<Scalar> data;
+  min_quad_with_fixed_precompute(Q,b,SparseMatrix<Scalar>(),true,data);
+  W.resize(A.rows(),bc.cols());
+  const VectorXS B = VectorXS::Zero(A.rows(),1);
+  for(int w = 0;w<bc.cols();w++)
+  {
+    const VectorXS bcw = bc.col(w);
+    VectorXS Ww;
+    if(!min_quad_with_fixed_solve(data,B,bcw,VectorXS(),Ww))
+    {
+      return false;
+    }
+    W.col(w) = Ww;
+  }
+  return true;
+}
+
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template bool igl::harmonic<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
-template bool igl::harmonic<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-template bool igl::harmonic<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template bool igl::harmonic<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template bool igl::harmonic<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template bool igl::harmonic<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+
+template bool igl::harmonic<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template bool igl::harmonic<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template bool igl::harmonic<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/harmonic.h

@@ -20,7 +20,6 @@ namespace igl
   //   b  #b boundary indices into V
   //   bc #b by #W list of boundary values
   //   k  power of harmonic operation (1: harmonic, 2: biharmonic, etc)
-  //   use_cotan_weights  flag to indicate whether the laplacian operator is uniform
   // Outputs:
   //   W  #V by #W list of weights
   //
@@ -36,9 +35,32 @@ namespace igl
     const Eigen::PlainObjectBase<Derivedb> & b,
     const Eigen::PlainObjectBase<Derivedbc> & bc,
     const int k,
-    const bool use_cotan_weights,
     Eigen::PlainObjectBase<DerivedW> & W);
-};
+
+ // Compute harmonic map using uniform laplacian operator
+ //
+ //
+ // Inputs:
+ //   F  #F by simplex-size list of element indices
+ //   b  #b boundary indices into V
+ //   bc #b by #W list of boundary values
+ //   k  power of harmonic operation (1: harmonic, 2: biharmonic, etc)
+ // Outputs:
+ //   W  #V by #W list of weights
+ //
+ template <
+   typename DerivedF,
+   typename Derivedb,
+   typename Derivedbc,
+   typename DerivedW>
+ IGL_INLINE bool harmonic(
+   const Eigen::PlainObjectBase<DerivedF> & F,
+   const Eigen::PlainObjectBase<Derivedb> & b,
+   const Eigen::PlainObjectBase<Derivedbc> & bc,
+   int k,
+   Eigen::PlainObjectBase<DerivedW> & W);
+ };
+
 #ifndef IGL_STATIC_LIBRARY
 #include "harmonic.cpp"
 #endif

tutorial/401_BiharmonicDeformation/main.cpp

@@ -30,11 +30,11 @@ bool pre_draw(igl::viewer::Viewer & viewer)
   {
     MatrixXd D;
     MatrixXd D_bc = U_bc_anim - V_bc;
-    igl::harmonic(V,F,b,D_bc,2,true,D);
+    igl::harmonic(V,F,b,D_bc,2,D);
     U = V+D;
   }else
   {
-    igl::harmonic(V,F,b,U_bc_anim,2,true,U);
+    igl::harmonic(V,F,b,U_bc_anim,2.,U);
   }
   viewer.data.set_vertices(U);
   viewer.data.compute_normals();

tutorial/402_PolyharmonicDeformation/main.cpp

@@ -21,7 +21,7 @@ bool pre_draw(igl::viewer::Viewer & viewer)
   using namespace Eigen;
   if(resolve)
   {
-    igl::harmonic(V,F,b,bc,k,true,Z);
+    igl::harmonic(V,F,b,bc,k,Z);
     resolve = false;
   }
   U.col(2) = z_max*Z;

tutorial/403_BoundedBiharmonicWeights/main.cpp

@@ -32,7 +32,7 @@
 
 #include "tutorial_shared_path.h"
 
-typedef 
+typedef
   std::vector<Eigen::Quaterniond,Eigen::aligned_allocator<Eigen::Quaterniond> >
   RotationList;
 
@@ -78,7 +78,7 @@ bool pre_draw(igl::viewer::Viewer & viewer)
     MatrixXd CT;
     MatrixXi BET;
     igl::deform_skeleton(C,BE,T,CT,BET);
-    
+
     viewer.data.set_vertices(U);
     viewer.data.set_edges(CT,BET,sea_green);
     viewer.data.compute_normals();

tutorial/501_HarmonicParam/main.cpp

@@ -44,7 +44,7 @@ int main(int argc, char *argv[])
   igl::map_vertices_to_circle(V,bnd,bnd_uv);
 
   // Harmonic parametrization for the internal vertices
-  igl::harmonic(V,F,bnd,bnd_uv,1,true,V_uv);
+  igl::harmonic(V,F,bnd,bnd_uv,1,V_uv);
 
   // Scale UV to make the texture more clear
   V_uv *= 5;

tutorial/503_ARAPParam/main.cpp

@@ -52,7 +52,7 @@ int main(int argc, char *argv[])
   Eigen::MatrixXd bnd_uv;
   igl::map_vertices_to_circle(V,bnd,bnd_uv);
 
-  igl::harmonic(V,F,bnd,bnd_uv,1,true,initial_guess);
+  igl::harmonic(V,F,bnd,bnd_uv,1,initial_guess);
 
   // Add dynamic regularization to avoid to specify boundary conditions
   igl::ARAPData arap_data;