isdiagonal test, harmonic abstracted away from laplacian

Former-commit-id: ecc61370649b6d71521dd766b94af558d2134a7a

include/igl/harmonic.cpp

@@ -6,13 +6,15 @@
 // 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/.
 #include "harmonic.h"
-#include "cotmatrix.h"
-#include "massmatrix.h"
-#include "invert_diag.h"
 #include "adjacency_matrix.h"
-#include "sum.h"
+#include "cotmatrix.h"
 #include "diag.h"
+#include "invert_diag.h"
+#include "isdiag.h"
+#include "massmatrix.h"
 #include "min_quad_with_fixed.h"
+#include "speye.h"
+#include "sum.h"
 #include <Eigen/Sparse>
 
 template <
@@ -31,9 +33,7 @@ IGL_INLINE bool igl::harmonic(
 {
   using namespace Eigen;
   typedef typename DerivedV::Scalar Scalar;
-  typedef Matrix<Scalar,Dynamic,1> VectorXS;
-  SparseMatrix<Scalar> Q;
-  SparseMatrix<Scalar> L,M,Mi;
+  SparseMatrix<Scalar> L,M;
   cotmatrix(V,F,L);
   switch(F.cols())
   {
@@ -45,28 +45,7 @@ IGL_INLINE bool igl::harmonic(
       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();
-  }
-
-  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);
-    VectorXS Ww;
-    if(!min_quad_with_fixed_solve(data,B,bcw,VectorXS(),Ww))
-    {
-      return false;
-    }
-    W.col(w) = Ww;
-  }
-  return true;
+  return harmonic(L,M,b,bc,k,W);
 }
 
 template <
@@ -83,8 +62,6 @@ IGL_INLINE bool igl::harmonic(
 {
   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
@@ -93,12 +70,46 @@ IGL_INLINE bool igl::harmonic(
   // Convert row sums into diagonal of sparse matrix
   SparseMatrix<Scalar> Adiag;
   diag(Asum,Adiag);
-  // Build uniform laplacian
-  Q = -A+Adiag;
+  SparseMatrix<Scalar> L = A-Adiag;
+  SparseMatrix<Scalar> M;
+  speye(L.rows(),M);
+  return harmonic(L,M,b,bc,k,W);
+}
+
+template <
+  typename DerivedL,
+  typename DerivedM,
+  typename Derivedb,
+  typename Derivedbc,
+  typename DerivedW>
+IGL_INLINE bool igl::harmonic(
+  const Eigen::SparseMatrix<DerivedL> & L,
+  const Eigen::SparseMatrix<DerivedM> & M,
+  const Eigen::PlainObjectBase<Derivedb> & b,
+  const Eigen::PlainObjectBase<Derivedbc> & bc,
+  const int k,
+  Eigen::PlainObjectBase<DerivedW> & W)
+{
+  const int n = L.rows();
+  assert(n == L.cols() && "L must be square");
+  assert(n == M.cols() && "M must be same size as L");
+  assert(n == M.rows() && "M must be square");
+  assert(igl::isdiag(M) && "Mass matrix should be diagonal");
+
+  Eigen::SparseMatrix<DerivedM> Mi;
+  invert_diag(M,Mi);
+  Eigen::SparseMatrix<DerivedL> Q;
+  Q = -L;
+  for(int p = 1;p<k;p++)
+  {
+    Q = (Q*Mi*-L).eval();
+  }
+  typedef DerivedL Scalar;
   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);
+  min_quad_with_fixed_precompute(Q,b,Eigen::SparseMatrix<Scalar>(),true,data);
+  W.resize(n,bc.cols());
+  typedef Eigen::Matrix<Scalar,Eigen::Dynamic,1> VectorXS;
+  const VectorXS B = VectorXS::Zero(n,1);
   for(int w = 0;w<bc.cols();w++)
   {
     const VectorXS bcw = bc.col(w);

include/igl/harmonic.h

@@ -9,6 +9,7 @@
 #define IGL_HARMONIC_H
 #include "igl_inline.h"
 #include <Eigen/Core>
+#include <Eigen/Sparse>
 namespace igl
 {
   // Compute k-harmonic weight functions "coordinates".
@@ -36,30 +37,51 @@ namespace igl
     const Eigen::PlainObjectBase<Derivedbc> & bc,
     const int k,
     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,
-   const int k,
-   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,
+    const int k,
+    Eigen::PlainObjectBase<DerivedW> & W);
+  // Compute a harmonic map using a given Laplacian and mass matrix
+  //
+  // Inputs:
+  //   L  #V by #V discrete (integrated) Laplacian  
+  //   M  #V by #V mass matrix
+  //   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 #V list of weights
+  template <
+    typename DerivedL,
+    typename DerivedM,
+    typename Derivedb,
+    typename Derivedbc,
+    typename DerivedW>
+  IGL_INLINE bool harmonic(
+    const Eigen::SparseMatrix<DerivedL> & L,
+    const Eigen::SparseMatrix<DerivedM> & M,
+    const Eigen::PlainObjectBase<Derivedb> & b,
+    const Eigen::PlainObjectBase<Derivedbc> & bc,
+    const int k,
+    Eigen::PlainObjectBase<DerivedW> & W);
+};
 
 #ifndef IGL_STATIC_LIBRARY
 #include "harmonic.cpp"

include/igl/isdiag.cpp

@@ -0,0 +1,27 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 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 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "isdiag.h"
+
+template <typename Scalar>
+IGL_INLINE bool igl::isdiag(const Eigen::SparseMatrix<Scalar> & A)
+{
+  // Iterate over outside of A
+  for(int k=0; k<A.outerSize(); ++k)
+  {
+    // Iterate over inside
+    for(typename Eigen::SparseMatrix<Scalar>::InnerIterator it (A,k); it; ++it)
+    {
+      if(it.row() != it.col() && it.value()!=0)
+      {
+        return false;
+      }
+    }
+  }
+  return true;
+}
+

include/igl/isdiag.h

@@ -0,0 +1,26 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 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 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_ISDIAG_H
+#define IGL_ISDIAG_H
+#include "igl_inline.h"
+#include <Eigen/Sparse>
+namespace igl
+{
+  // Determine if a given matrix is diagonal: all non-zeros lie on the
+  // main diagonal.
+  //
+  // Inputs:
+  //   A  m by n sparse matrix
+  // Returns true iff and only if the matrix is diagonal.
+  template <typename Scalar>
+  IGL_INLINE bool isdiag(const Eigen::SparseMatrix<Scalar> & A);
+};
+#ifndef IGL_STATIC_LIBRARY
+#  include "isdiag.cpp"
+#endif
+#endif

tutorial/tutorial.md.REMOVED.git-id

@@ -1 +1 @@
-129cef4e745b13b5c567053630e25f92dfead7e1
+abaa308d298147242a29519cbbb3f30c69d0d06e