massmatrix

Former-commit-id: fb3534f2bddb9b57d54429bc07c2a50a723fffaa

include/igl/diag.cpp

@@ -11,6 +11,7 @@ IGL_INLINE void igl::diag(
   int m = X.rows();
   int n = X.cols();
   V = Eigen::SparseVector<T>((m>n?n:m));
+  V.reserve(V.size());
 
   // Iterate over outside
   for(int k=0; k<X.outerSize(); ++k)
@@ -26,6 +27,30 @@ IGL_INLINE void igl::diag(
   }
 }
 
+template <typename T,typename DerivedV>
+IGL_INLINE void igl::diag(
+  const Eigen::SparseMatrix<T>& X, 
+  Eigen::MatrixBase<DerivedV> & V)
+{
+  // Get size of input
+  int m = X.rows();
+  int n = X.cols();
+  V.derived().resize((m>n?n:m),1);
+
+  // Iterate over outside
+  for(int k=0; k<X.outerSize(); ++k)
+  {
+    // Iterate over inside
+    for(typename Eigen::SparseMatrix<T>::InnerIterator it (X,k); it; ++it)
+    {
+      if(it.col() == it.row())
+      {
+        V(it.col()) = it.value();
+      }
+    }
+  }
+}
+
 template <typename T>
 IGL_INLINE void igl::diag(
   const Eigen::SparseVector<T>& V,
@@ -33,16 +58,34 @@ IGL_INLINE void igl::diag(
 {
   // clear and resize output
   Eigen::DynamicSparseMatrix<T, Eigen::RowMajor> dyn_X(V.size(),V.size());
-
+  dyn_X.reserve(V.size());
   // loop over non-zeros
   for(typename Eigen::SparseVector<T>::InnerIterator it(V); it; ++it)
   {
     dyn_X.coeffRef(it.index(),it.index()) += it.value();
   }
+  X = Eigen::SparseMatrix<T>(dyn_X);
+}
 
+template <typename T, typename DerivedV>
+IGL_INLINE void igl::diag(
+  const Eigen::MatrixBase<DerivedV> & V,
+  Eigen::SparseMatrix<T>& X)
+{
+  assert(V.rows() == 1 || V.cols() == 1);
+  // clear and resize output
+  Eigen::DynamicSparseMatrix<T, Eigen::RowMajor> dyn_X(V.size(),V.size());
+  dyn_X.reserve(V.size());
+  // loop over non-zeros
+  for(int i = 0;i<V.size();i++)
+  {
+    dyn_X.coeffRef(i,i) += V[i];
+    i++;
+  }
   X = Eigen::SparseMatrix<T>(dyn_X);
 }
 
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
+template void igl::diag<double, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/diag.h

@@ -20,6 +20,10 @@ namespace igl
   IGL_INLINE void diag(
     const Eigen::SparseMatrix<T>& X, 
     Eigen::SparseVector<T>& V);
+  template <typename T,typename DerivedV>
+  IGL_INLINE void diag(
+    const Eigen::SparseMatrix<T>& X, 
+    Eigen::MatrixBase<DerivedV>& V);
   // Templates:
   //   T  should be a eigen sparse matrix primitive type like int or double
   // Inputs:
@@ -30,6 +34,10 @@ namespace igl
   IGL_INLINE void diag(
     const Eigen::SparseVector<T>& V,
     Eigen::SparseMatrix<T>& X);
+  template <typename T, typename DerivedV>
+  IGL_INLINE void diag(
+    const Eigen::MatrixBase<DerivedV>& V,
+    Eigen::SparseMatrix<T>& X);
 }
 
 #ifdef IGL_HEADER_ONLY

include/igl/massmatrix.cpp

@@ -0,0 +1,153 @@
+#include "massmatrix.h"
+#include "normalize_rows.h"
+#include "sparse.h"
+#include "repmat.h"
+#include <Eigen/Geometry>
+#include <iostream>
+
+template <typename DerivedV, typename DerivedF, typename Scalar>
+IGL_INLINE void igl::massmatrix(
+  const Eigen::MatrixBase<DerivedV> & V, 
+  const Eigen::MatrixBase<DerivedF> & F, 
+  const MassMatrixType type,
+  Eigen::SparseMatrix<Scalar>& M)
+{
+  using namespace Eigen;
+  using namespace std;
+  assert(type!=MASSMATRIX_FULL);
+
+  const int n = V.rows();
+  const int m = F.rows();
+  const int simplex_size = F.cols();
+
+  Matrix<int,Dynamic,1> MI;
+  Matrix<int,Dynamic,1> MJ;
+  Matrix<Scalar,Dynamic,1> MV;
+  if(simplex_size == 3)
+  {
+    // Triangles
+    // edge lengths numbered same as opposite vertices
+    Matrix<Scalar,Dynamic,3> l(m,3);
+    // loop over faces
+    for(int i = 0;i<m;i++)
+    {
+      l(i,0) = sqrt((V.row(F(i,1))-V.row(F(i,2))).array().pow(2).sum());
+      l(i,1) = sqrt((V.row(F(i,2))-V.row(F(i,0))).array().pow(2).sum());
+      l(i,2) = sqrt((V.row(F(i,0))-V.row(F(i,1))).array().pow(2).sum());
+    }
+    // semiperimeters
+    Matrix<Scalar,Dynamic,1> s = l.rowwise().sum()*0.5;
+    assert(s.rows() == m);
+    // Heron's forumal for area
+    Matrix<Scalar,Dynamic,1> dblA(m);
+    for(int i = 0;i<m;i++)
+    {
+      dblA(i) = 2.0*sqrt(s(i)*(s(i)-l(i,0))*(s(i)-l(i,1))*(s(i)-l(i,2)));
+    }
+
+    switch(type)
+    {
+      case MASSMATRIX_BARYCENTRIC:
+        // diagonal entries for each face corner
+        MI.resize(m*3,1); MJ.resize(m*3,1); MV.resize(m*3,1);
+        MI.block(0*m,0,m,1) = F.col(0);
+        MI.block(1*m,0,m,1) = F.col(1);
+        MI.block(2*m,0,m,1) = F.col(2);
+        MJ = MI;
+        repmat(dblA,3,1,MV);
+        MV.array() /= 6.0;
+        break;
+      case MASSMATRIX_VORONOI:
+        {
+          // diagonal entries for each face corner
+          // http://www.alecjacobson.com/weblog/?p=874
+          MI.resize(m*3,1); MJ.resize(m*3,1); MV.resize(m*3,1);
+          MI.block(0*m,0,m,1) = F.col(0);
+          MI.block(1*m,0,m,1) = F.col(1);
+          MI.block(2*m,0,m,1) = F.col(2);
+          MJ = MI;
+
+          // Holy shit this needs to be cleaned up and optimized
+          Matrix<Scalar,Dynamic,3> cosines(m,3);
+          cosines.col(0) = 
+            (l.col(2).array().pow(2)+l.col(1).array().pow(2)-l.col(0).array().pow(2))/(l.col(1).array()*l.col(2).array()*2.0);
+          cosines.col(1) = 
+            (l.col(0).array().pow(2)+l.col(2).array().pow(2)-l.col(1).array().pow(2))/(l.col(2).array()*l.col(0).array()*2.0);
+          cosines.col(2) = 
+            (l.col(1).array().pow(2)+l.col(0).array().pow(2)-l.col(2).array().pow(2))/(l.col(0).array()*l.col(1).array()*2.0);
+          Matrix<Scalar,Dynamic,3> barycentric = cosines.array() * l.array();
+
+          WRONG NORMALIZE THIS IS L2 NORMALIZE
+
+          //normalize_rows(barycentric,barycentric);
+          cout<<"barycentric=["<<barycentric<<"];"<<endl;
+          Matrix<Scalar,Dynamic,3> partial = barycentric;
+          partial.col(0).array() *= dblA.array() * 0.5;
+          partial.col(1).array() *= dblA.array() * 0.5;
+          partial.col(2).array() *= dblA.array() * 0.5;
+          cout<<"partial=["<<partial<<"];"<<endl;
+          Matrix<Scalar,Dynamic,3> quads(partial.rows(),partial.cols());
+          quads.col(0) = (partial.col(1)+partial.col(2))*0.5;
+          quads.col(1) = (partial.col(2)+partial.col(0))*0.5;
+          quads.col(2) = (partial.col(0)+partial.col(1))*0.5;
+          cout<<"quads=["<<quads<<"];"<<endl;
+
+          Matrix<Scalar,Dynamic,3> quads_clamped(quads.rows(),quads.cols());
+          quads_clamped.col(0) = (cosines.col(0).array()<0).select( 0.25*dblA,quads.col(0));
+          quads_clamped.col(1) = (cosines.col(0).array()<0).select(0.125*dblA,quads.col(1));
+          quads_clamped.col(2) = (cosines.col(0).array()<0).select(0.125*dblA,quads.col(2));
+          quads_clamped.col(0) = (cosines.col(1).array()<0).select(0.125*dblA,quads.col(0));
+          quads_clamped.col(1) = (cosines.col(1).array()<0).select(0.125*dblA,quads.col(1));
+          quads_clamped.col(2) = (cosines.col(1).array()<0).select( 0.25*dblA,quads.col(2));
+          quads_clamped.col(0) = (cosines.col(2).array()<0).select(0.125*dblA,quads.col(0));
+          quads_clamped.col(1) = (cosines.col(2).array()<0).select(0.125*dblA,quads.col(1));
+          quads_clamped.col(2) = (cosines.col(2).array()<0).select( 0.25*dblA,quads.col(2));
+
+          MV.block(0*m,0,m,1) = quads_clamped.col(0);
+          MV.block(1*m,0,m,1) = quads_clamped.col(1);
+          MV.block(2*m,0,m,1) = quads_clamped.col(2);
+          
+          break;
+        }
+      case MASSMATRIX_FULL:
+        assert(false);
+        break;
+      default:
+        assert(false);
+    }
+
+  }else if(simplex_size == 4)
+  {
+    assert(V.cols() == 3);
+    assert(type == MASSMATRIX_BARYCENTRIC);
+    MI.resize(m*4,1); MJ.resize(m*4,1); MV.resize(m*4,1);
+    MI.block(0*m,0,m,1) = F.col(0);
+    MI.block(1*m,0,m,1) = F.col(1);
+    MI.block(2*m,0,m,1) = F.col(2);
+    MI.block(3*m,0,m,1) = F.col(3);
+    MJ = MI;
+    // loop over tets
+    for(int i = 0;i<m;i++)
+    {
+      // http://en.wikipedia.org/wiki/Tetrahedron#Volume
+      Matrix<Scalar,3,1> v0m3 = V.row(F(i,0)) - V.row(F(i,3));
+      Matrix<Scalar,3,1> v1m3 = V.row(F(i,1)) - V.row(F(i,3));
+      Matrix<Scalar,3,1> v2m3 = V.row(F(i,2)) - V.row(F(i,3));
+      Scalar v = fabs(v0m3.dot(v1m3.cross(v2m3)))/6.0;
+      MV(i+0*m) = v/4.0;
+      MV(i+1*m) = v/4.0;
+      MV(i+2*m) = v/4.0;
+      MV(i+3*m) = v/4.0;
+    }
+  }else
+  {
+    // Unsupported simplex size
+    assert(false);
+  }
+  sparse(MI,MJ,MV,n,n,M);
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+template void igl::massmatrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::MassMatrixType, Eigen::SparseMatrix<double, 0, int>&);
+#endif

include/igl/massmatrix.h

@@ -0,0 +1,53 @@
+#ifndef IGL_MASSMATRIX_H
+#define IGL_MASSMATRIX_H
+#include "igl_inline.h"
+
+#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+namespace igl 
+{
+
+  enum MassMatrixType
+  {
+    MASSMATRIX_BARYCENTRIC,
+    MASSMATRIX_VORONOI,
+    MASSMATRIX_FULL
+  };
+#define NUM_MASSMATRIXTYPE 3
+
+  // Constructs the mass (area) matrix for a given mesh (V,F).
+  //
+  // Templates:
+  //   DerivedV  derived type of eigen matrix for V (e.g. derived from
+  //     MatirxXd)
+  //   DerivedF  derived type of eigen matrix for F (e.g. derived from
+  //     MatrixXi)
+  //   Scalar  scalar type for eigen sparse matrix (e.g. double)
+  // Inputs:
+  //   V  #V by dim list of mesh vertex positions
+  //   F  #F by simplex_size list of mesh faces (must be triangles)
+  //   type  one of the following ints:
+  //     IGL_MASSMATRIX_BARYCENTRIC  barycentric
+  //     IGL_MASSMATRIX_VORONOI voronoi-hybrid {default}
+  //     IGL_MASSMATRIX_FULL full {not implemented}
+  // Outputs: 
+  //   M  #V by #V mass matrix
+  //
+  // See also: adjacency_matrix
+  //
+  template <typename DerivedV, typename DerivedF, typename Scalar>
+  IGL_INLINE void massmatrix(
+    const Eigen::MatrixBase<DerivedV> & V, 
+    const Eigen::MatrixBase<DerivedF> & F, 
+    const MassMatrixType type,
+    Eigen::SparseMatrix<Scalar>& M);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "massmatrix.cpp"
+#endif
+
+#endif
+

include/igl/normalize_rows.cpp

@@ -2,8 +2,8 @@
 
 template <typename DerivedV>
 IGL_INLINE void igl::normalize_rows(
-                               const Eigen::PlainObjectBase<DerivedV>& A,
-                               Eigen::PlainObjectBase<DerivedV> & B)
+  const Eigen::PlainObjectBase<DerivedV>& A,
+  Eigen::PlainObjectBase<DerivedV> & B)
 {
   // Resize output
   B.resize(A.rows(),A.cols());
@@ -19,4 +19,5 @@ IGL_INLINE void igl::normalize_rows(
 // generated by autoexplicit.sh
 template void igl::normalize_rows<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::normalize_rows<Eigen::Matrix<double, -1, 3, 1, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&);
+template void igl::normalize_rows<Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 #endif

include/igl/normalize_rows.h

@@ -13,8 +13,8 @@ namespace igl
   //   B  #rows by k input matrix, can be the same as A
   template <typename DerivedV>
   IGL_INLINE void normalize_rows(
-                                 const Eigen::PlainObjectBase<DerivedV>& A,
-                                 Eigen::PlainObjectBase<DerivedV> & B);
+   const Eigen::PlainObjectBase<DerivedV>& A,
+   Eigen::PlainObjectBase<DerivedV> & B);
 }
 
 #ifdef IGL_HEADER_ONLY

include/igl/per_vertex_normals.cpp

@@ -5,9 +5,9 @@
 
 template <typename DerivedV, typename DerivedF>
 IGL_INLINE void igl::per_vertex_normals(
-                                   const Eigen::PlainObjectBase<DerivedV>& V,
-                                   const Eigen::PlainObjectBase<DerivedF>& F,
-                                   Eigen::PlainObjectBase<DerivedV> & N)
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedV> & N)
 {
   Eigen::PlainObjectBase<DerivedV> PFN;
   igl::per_face_normals(V,F,PFN);
@@ -30,10 +30,10 @@ IGL_INLINE void igl::per_vertex_normals(
 
 template <typename DerivedV, typename DerivedF>
 IGL_INLINE void igl::per_vertex_normals(
-                                        const Eigen::PlainObjectBase<DerivedV>& V,
-                                        const Eigen::PlainObjectBase<DerivedF>& F,
-                                        const Eigen::PlainObjectBase<DerivedV>& FN,
-                                        Eigen::PlainObjectBase<DerivedV> & N)
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  const Eigen::PlainObjectBase<DerivedV>& FN,
+  Eigen::PlainObjectBase<DerivedV> & N)
 {
   // Resize for output
   N = Eigen::PlainObjectBase<DerivedV>::Zero(V.rows(),3);

include/igl/repmat.cpp

@@ -55,4 +55,5 @@ IGL_INLINE void igl::repmat(
 // Explicit template specialization
 // generated by autoexplicit.sh
 template void igl::repmat<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::repmat<Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, int, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/sparse.cpp

@@ -49,4 +49,5 @@ IGL_INLINE void igl::sparse(
 
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
+template void igl::sparse<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, double>(Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1> const&, unsigned long, unsigned long, Eigen::SparseMatrix<double, 0, int>&);
 #endif