added option for uniform laplacian operator

Former-commit-id: f638f86debd52c8046b79e5cf8cf36b4b36735f4

include/igl/harmonic.cpp

@@ -9,6 +9,9 @@
 #include "cotmatrix.h"
 #include "massmatrix.h"
 #include "invert_diag.h"
+#include "adjacency_matrix.h"
+#include "sum.h"
+#include "diag.h"
 #include "min_quad_with_fixed.h"
 #include <Eigen/Sparse>
 
@@ -24,33 +27,50 @@ IGL_INLINE bool igl::harmonic(
   const Eigen::PlainObjectBase<Derivedb> & b,
   const Eigen::PlainObjectBase<Derivedbc> & bc,
   const int k,
+  const bool use_cotan_weights,
   Eigen::PlainObjectBase<DerivedW> & W)
 {
   using namespace Eigen;
   typedef typename DerivedV::Scalar Scalar;
   typedef Matrix<Scalar,Dynamic,1> VectorXS;
-  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> 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);
       break;
+    }
+    invert_diag(M,Mi);
+    Q = -L;
+    for(int p = 1;p<k;p++)
+    {
+      Q = (Q*Mi*-L).eval();
+    }
   }
-  invert_diag(M,Mi);
-  SparseMatrix<Scalar> Q = -L;
-  for(int p = 1;p<k;p++)
+  else
   {
-    Q = (Q*Mi*-L).eval();
+    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;
   }
-  const VectorXS B = VectorXS::Zero(V.rows(),1);
   min_quad_with_fixed_data<Scalar> data;
   min_quad_with_fixed_precompute(Q,b,SparseMatrix<Scalar>(),true,data);
   W.resize(V.rows(),bc.cols());
+  const VectorXS B = VectorXS::Zero(V.rows(),1);
   for(int w = 0;w<bc.cols();w++)
   {
     const VectorXS bcw = bc.col(w);
@@ -66,7 +86,7 @@ IGL_INLINE bool igl::harmonic(
 
 #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, 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, 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> >&);
 #endif

include/igl/harmonic.h

@@ -1,15 +1,15 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #ifndef IGL_HARMONIC_H
 #define IGL_HARMONIC_H
 #include "igl_inline.h"
 #include <Eigen/Core>
-namespace igl 
+namespace igl
 {
   // Compute k-harmonic weight functions "coordinates".
   //
@@ -20,6 +20,7 @@ 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
   //
@@ -35,6 +36,7 @@ namespace igl
     const Eigen::PlainObjectBase<Derivedb> & b,
     const Eigen::PlainObjectBase<Derivedbc> & bc,
     const int k,
+    const bool use_cotan_weights,
     Eigen::PlainObjectBase<DerivedW> & W);
 };
 #ifndef IGL_STATIC_LIBRARY