Merge branch 'master-upstream' into python_bindings

Former-commit-id: 249916650cf815145cb2b998dc78621fce2ef7e1

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>
 
@@ -29,6 +32,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;
   cotmatrix(V,F,L);
   switch(F.cols())
@@ -39,18 +43,19 @@ IGL_INLINE bool igl::harmonic(
     case 4:
     default:
       massmatrix(V,F,MASSMATRIX_TYPE_BARYCENTRIC,M);
-      break;
+    break;
   }
   invert_diag(M,Mi);
-  SparseMatrix<Scalar> Q = -L;
+  Q = -L;
   for(int p = 1;p<k;p++)
   {
     Q = (Q*Mi*-L).eval();
   }
-  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);
@@ -64,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,
+  const 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, 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

@@ -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".
   //
@@ -36,7 +36,31 @@ 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);
+ };
+
 #ifndef IGL_STATIC_LIBRARY
 #include "harmonic.cpp"
 #endif

tutorial/401_BiharmonicDeformation/main.cpp

@@ -34,7 +34,7 @@ bool pre_draw(igl::viewer::Viewer & viewer)
     U = V+D;
   }else
   {
-    igl::harmonic(V,F,b,U_bc_anim,2,U);
+    igl::harmonic(V,F,b,U_bc_anim,2.,U);
   }
   viewer.data.set_vertices(U);
   viewer.data.compute_normals();
@@ -66,7 +66,7 @@ int main(int argc, char *argv[])
   VectorXi S;
   igl::readDMAT(TUTORIAL_SHARED_PATH "/decimated-max-selection.dmat",S);
   igl::colon<int>(0,V.rows()-1,b);
-  b.conservativeResize(stable_partition( b.data(), b.data()+b.size(), 
+  b.conservativeResize(stable_partition( b.data(), b.data()+b.size(),
    [&S](int i)->bool{return S(i)>=0;})-b.data());
 
   // Boundary conditions directly on deformed positions

tutorial/402_PolyharmonicDeformation/main.cpp

@@ -68,7 +68,7 @@ int main(int argc, char *argv[])
   VectorXb is_inner = (V.rowwise().norm().array()-0.15)<1e-15;
   VectorXb in_b = is_outer.array() || is_inner.array();
   igl::colon<int>(0,V.rows()-1,b);
-  b.conservativeResize(stable_partition( b.data(), b.data()+b.size(), 
+  b.conservativeResize(stable_partition( b.data(), b.data()+b.size(),
    [&in_b](int i)->bool{return in_b(i);})-b.data());
   bc.resize(b.size(),1);
   for(int bi = 0;bi<b.size();bi++)

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

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