bug fixes and clean up for cotmatrix,

diag, find, sum, sort, edges, adjacency_matrix sparse matrix functions to match matlab


Former-commit-id: 7c030344122c050936bed1154174da0affd961d4

adjacency_matrix.h

@@ -0,0 +1,67 @@
+#ifndef IGL_ADJACENCY_MATRIX_H
+#define IGL_ADJACENCY_MATRIX_H
+
+#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+namespace igl 
+{
+  // Constructs the graph adjacency matrix  of a given mesh (V,F)
+  // Templates:
+  //   T  should be a eigen sparse matrix primitive type like int or double
+  // Inputs:
+  //   F  #F by dim list of mesh faces (must be triangles)
+  // Outputs: 
+  //   A  max(F) by max(F) cotangent matrix, each row i corresponding to V(i,:)
+  //
+  // Example:
+  //   // Mesh in (V,F)
+  //   Eigen::SparseMatrix<double> A;
+  //   adjacency_matrix(F,A);
+  //   // sum each row 
+  //   SparseVector<double> Asum;
+  //   sum(A,1,Asum);
+  //   // Convert row sums into diagonal of sparse matrix
+  //   SparseMatrix<double> Adiag;
+  //   diag(Asum,Adiag);
+  //   // Build uniform laplacian
+  //   SparseMatrix<double> U;
+  //   U = A-Adiag;
+  //
+  // See also: edges, cotmatrix, diag
+  template <typename T>
+  void adjacency_matrix(
+    const Eigen::MatrixXi & F, 
+    Eigen::SparseMatrix<T>& A);
+}
+
+// Implementation
+#include "verbose.h"
+
+template <typename T>
+void igl::adjacency_matrix(
+  const Eigen::MatrixXi & F, 
+  Eigen::SparseMatrix<T>& A)
+{
+  Eigen::DynamicSparseMatrix<T, Eigen::RowMajor> 
+    dyn_A(F.maxCoeff()+1, F.maxCoeff()+1);
+
+  // Loop over faces
+  for(int i = 0;i<F.rows();i++)
+  {
+    // Loop over this face
+    for(int j = 0;j<F.cols();j++)
+    {
+      // Get indices of edge: s --> d
+      int s = F(i,j);
+      int d = F(i,(j+1)%F.cols());
+      dyn_A.coeffRef(s, d) = 1;
+      dyn_A.coeffRef(d, s) = 1;
+    }
+  }
+
+  A = Eigen::SparseMatrix<T>(dyn_A);
+}
+
+#endif

cotmatrix.h

@@ -1,80 +1,129 @@
-//
-//  IGL Lib - Simple C++ mesh library 
-//
-//  Copyright 2011, Daniele Panozzo. All rights reserved.
-
-#ifndef COTMATRIX_H
-#define COTMATRIX_H
+#ifndef IGL_COTMATRIX_H
+#define IGL_COTMATRIX_H
 
 #define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
 #include <Eigen/Dense>
 #include <Eigen/Sparse>
 
+// History:
+//  Used const references rather than copying the entire mesh 
+//    Alec 9 October 2011
+//  removed cotan (uniform weights) optional parameter it was building a buggy
+//    half of the uniform laplacian, please see adjacency_matrix istead 
+//    Alec 9 October 2011
 
 namespace igl 
 {
-    void computeCotWeights(Eigen::Vector3d& v1, Eigen::Vector3d& v2, Eigen::Vector3d& v3, double& cotAlpha, double& cotBeta, double& cotGamma)
-    {
-        Eigen::Vector3d v12 = v2-v1;
-        Eigen::Vector3d v13 = v3-v1;
-        Eigen::Vector3d v23 = v3-v2;
-        
-        double halfArea = (v12.cross(v13)).norm();//squaredNorm();
-        
-        //improve numerical stability
-        const double cotTolerance = 1e-10;
-        if(halfArea < cotTolerance)
-        {
-            std::cout << "Cot weights are close to singular!" << std::endl;
-            halfArea = cotTolerance;
-        }
-        
-        cotAlpha = (v12.dot(v13)) / halfArea /2;
-        cotBeta =  (v23.dot(-v12)) / halfArea /2;
-        cotGamma = (-v23.dot(-v13)) / halfArea /2;
-    }
-    
-    void cotmatrix (Eigen::MatrixXd V, Eigen::MatrixXi F, Eigen::SparseMatrix<double>& L, bool cotan = true)
+  // Constructs the cotangent stiffness matrix (discrete laplacian) for a given
+  // mesh (V,F).
+  // Inputs:
+  //   V  #V by dim list of mesh vertex positions
+  //   F  #F by 3 list of mesh faces (must be triangles)
+  // Outputs: 
+  //   L  #V by #V cotangent matrix, each row i corresponding to V(i,:)
+  //
+  // See also: adjacency_matrix
+  void cotmatrix(
+    const Eigen::MatrixXd & V, 
+    const Eigen::MatrixXi & F,
+    Eigen::SparseMatrix<double>& L);
+  // Helper function that computes the cotangent weights for each corner of a
+  // given triangle
+  // Inputs:
+  //   v1  position of corner #1 of triangle
+  //   v2  position of corner #2 of triangle
+  //   v3  position of corner #3 of triangle
+  // Outputs:
+  //   cot1  cotangent of angle at corner #1
+  //   cot2 cotangent of angle at corner #2
+  //   cot3  cotangent of angle at corner #3
+  //   
+  void computeCotWeights(
+    const Eigen::Vector3d& v1, 
+    const Eigen::Vector3d& v2, 
+    const Eigen::Vector3d& v3, 
+    double& cot1, 
+    double& cot2, 
+    double& cot3);
+}
+
+// Implementation
+
+// For error printing
+#include <cstdio>
+
+void igl::computeCotWeights(
+  const Eigen::Vector3d& v1, 
+  const Eigen::Vector3d& v2, 
+  const Eigen::Vector3d& v3, 
+  double& cot1, 
+  double& cot2, 
+  double& cot3)
+{
+  Eigen::Vector3d v12 = v2-v1;
+  Eigen::Vector3d v13 = v3-v1;
+  Eigen::Vector3d v23 = v3-v2;
+  
+  double halfArea = (v12.cross(v13)).norm();//squaredNorm();
+  
+  //improve numerical stability
+  const double cotTolerance = 1e-10;
+  if(halfArea < cotTolerance)
+  {
+    fprintf(stderr,"Cot weights are close to singular!\n");
+    halfArea = cotTolerance;
+  }
+  
+  cot1 = (v12.dot(v13)) / halfArea /2;
+  cot2 =  (v23.dot(-v12)) / halfArea /2;
+  cot3 = (-v23.dot(-v13)) / halfArea /2;
+}
+
+void igl::cotmatrix(
+  const Eigen::MatrixXd & V, 
+  const Eigen::MatrixXi & F, 
+  Eigen::SparseMatrix<double>& L)
+{
+  // Assumes vertices are 3D
+  assert(V.cols() == 3);
+  // Assumes faces are triangles
+  assert(F.cols() == 3);
+
+  Eigen::DynamicSparseMatrix<double, Eigen::RowMajor> dyn_L (V.rows(), V.rows());
+  
+  // Loop over triangles
+  for (unsigned i = 0; i < F.rows(); i++)
+  {
+    // Corner indices of this triangle
+    int vi1 = F(i,0);
+    int vi2 = F(i,1);
+    int vi3 = F(i,2);
+    // Grab corner positions of this triangle
+    Eigen::Vector3d v1(V(vi1,0), V(vi1,1), V(vi1,2));
+    Eigen::Vector3d v2(V(vi2,0), V(vi2,1), V(vi2,2));
+    Eigen::Vector3d v3(V(vi3,0), V(vi3,1), V(vi3,2));
+    // Compute cotangent of angles at each corner
+    double cot1, cot2, cot3;
+    computeCotWeights(v1, v2, v3, cot1, cot2, cot3);
+    // Throw each corner's cotangent at opposite edge (in both directions)
+    dyn_L.coeffRef(vi1, vi2) += cot3;
+    dyn_L.coeffRef(vi2, vi1) += cot3;
+    dyn_L.coeffRef(vi2, vi3) += cot1;
+    dyn_L.coeffRef(vi3, vi2) += cot1;
+    dyn_L.coeffRef(vi3, vi1) += cot2;
+    dyn_L.coeffRef(vi1, vi3) += cot2;
+  }
+
+  for (int k=0; k < dyn_L.outerSize(); ++k)
+  {
+    double tmp = 0.0f;
+    for(Eigen::DynamicSparseMatrix<double, Eigen::RowMajor>::InnerIterator it (dyn_L, k); it; ++it)
     {
-        Eigen::DynamicSparseMatrix<double, Eigen::RowMajor> dyn_L (V.rows(), V.rows());
-        
-        for (unsigned i = 0; i < F.rows(); i++)
-        {
-            int vi1 = F(i,0);
-            int vi2 = F(i,1);
-            int vi3 = F(i,2);
-            
-            if (cotan)
-            {
-                Eigen::Vector3d v1 (V(vi1,0), V(vi1,1), V(vi1,2));
-                Eigen::Vector3d v2 (V(vi2,0), V(vi2,1), V(vi2,2));
-                Eigen::Vector3d v3 (V(vi3,0), V(vi3,1), V(vi3,2));
-                double cot_a, cot_b, cot_g;
-                computeCotWeights (v1, v2, v3, cot_a, cot_b, cot_g);
-                
-                dyn_L.coeffRef (vi1, vi2) += cot_g;
-                dyn_L.coeffRef (vi2, vi1) += cot_g;
-                dyn_L.coeffRef (vi2, vi3) += cot_a;
-                dyn_L.coeffRef (vi3, vi2) += cot_a;
-                dyn_L.coeffRef (vi3, vi1) += cot_b;
-                dyn_L.coeffRef (vi1, vi3) += cot_b;
-            }
-            else
-            {
-                dyn_L.coeffRef (vi1, vi2) += 1.0;
-                dyn_L.coeffRef (vi2, vi3) += 1.0;
-                dyn_L.coeffRef (vi3, vi1) += 1.0;
-            }
-        }
-        for (int k=0; k < dyn_L.outerSize(); ++k)
-        {
-            double tmp = 0.0f;
-            for (Eigen::DynamicSparseMatrix<double, Eigen::RowMajor>::InnerIterator it (dyn_L, k); it; ++it)
-                tmp += it.value ();
-            dyn_L.coeffRef (k,k) = -tmp;
-        }
-        L = Eigen::SparseMatrix<double> (dyn_L);
+      tmp += it.value();
     }
+    dyn_L.coeffRef(k,k) = -tmp;
+  }
+  L = Eigen::SparseMatrix<double>(dyn_L);
 }
 
-#endif
+#endif

diag.h

@@ -0,0 +1,76 @@
+#ifndef IGL_DIAG_H
+#define IGL_DIAG_H
+#include <Eigen/Sparse>
+
+namespace igl
+{
+  // Either extracts the main diagonal of a matrix as a vector OR converts a
+  // vector into a matrix with vector along the main diagonal. Like matlab's
+  // diag function
+
+  // Templates:
+  //   T  should be a eigen sparse matrix primitive type like int or double
+  // Inputs:
+  //   X  an m by n sparse matrix
+  // Outputs:
+  //   V  a min(m,n) sparse vector
+  template <typename T>
+  void diag(
+    const Eigen::SparseMatrix<T>& X, 
+    Eigen::SparseVector<T>& V);
+  // Templates:
+  //   T  should be a eigen sparse matrix primitive type like int or double
+  // Inputs:
+  //  V  a m sparse vector
+  // Outputs:
+  //  X  a m by m sparse matrix
+  template <typename T>
+  void diag(
+    const Eigen::SparseVector<T>& V,
+    Eigen::SparseMatrix<T>& X);
+}
+
+// Implementation
+#include "verbose.h"
+
+template <typename T>
+void igl::diag(
+  const Eigen::SparseMatrix<T>& X, 
+  Eigen::SparseVector<T>& V)
+{
+  // Get size of input
+  int m = X.rows();
+  int n = X.cols();
+  V = Eigen::SparseVector<T>((m>n?n:m));
+
+  // 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.coeffRef(it.col()) += it.value();
+      }
+    }
+  }
+}
+
+template <typename T>
+void igl::diag(
+  const Eigen::SparseVector<T>& V,
+  Eigen::SparseMatrix<T>& X)
+{
+  // clear and resize output
+  Eigen::DynamicSparseMatrix<T, Eigen::RowMajor> dyn_X(V.size(),V.size());
+
+  // loop over non-zeros
+  for(typename SparseVector<T>::InnerIterator it(V); it; ++it)
+  {
+    dyn_X.coeffRef(it.index(),it.index()) += it.value();
+  }
+
+  X = Eigen::SparseMatrix<T>(dyn_X);
+}
+#endif

edges.h

@@ -0,0 +1,57 @@
+#ifndef IGL_EDGES_H
+#define IGL_EDGES_H
+
+#include <Eigen/Dense>
+
+namespace igl
+{
+  // Constructs a list of unique edges represented in a given mesh (V,F)
+  // Templates:
+  //   T  should be a eigen sparse matrix primitive type like int or double
+  // Inputs:
+  //   F  #F by 3 list of mesh faces (must be triangles)
+  //   or
+  //   T  #T x 4  matrix of indices of tet corners
+  // Outputs:
+  //   E #E by 2 list of edges in no particular order
+  //
+  // See also: adjacency_matrix
+  inline void edges(
+    const Eigen::MatrixXi& F, 
+    Eigen::MatrixXi& E);
+}
+
+// Implementation
+#include <map>
+#include "adjacency_matrix.h"
+
+inline void igl::edges(
+  const Eigen::MatrixXi& F, 
+  Eigen::MatrixXi& E)
+{
+  // build adjacency matrix
+  Eigen::SparseMatrix<int> A;
+  igl::adjacency_matrix(F,A);
+  // Number of non zeros should be twice number of edges
+  assert(A.nonZeros()%2 == 0);
+  // Resize to fit edges
+  E.resize(A.nonZeros()/2,2);
+  int i = 0;
+  // Iterate over outside
+  for(int k=0; k<A.outerSize(); ++k)
+  {
+    // Iterate over inside
+    for(Eigen::SparseMatrix<int>::InnerIterator it (A,k); it; ++it)
+    {
+      // only add edge in one direction
+      if(it.row()<it.col())
+      {
+        E(i,0) = it.row();
+        E(i,1) = it.col();
+        i++;
+      }
+    }
+  }
+}
+
+#endif

find.h

@@ -0,0 +1,95 @@
+#ifndef IGL_FIND_H
+#define IGL_FIND_H
+#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+namespace igl
+{
+  // Find the non-zero entries and there respective indices in a sparse matrix.
+  // Like matlab's [I,J,V] = find(X)
+  //
+  // Templates:
+  //   T  should be a eigen sparse matrix primitive type like int or double
+  // Input:
+  //   X  m by n matrix whose entries are to be found 
+  // Outputs:
+  //   I  nnz vector of row indices of non zeros entries in X
+  //   J  nnz vector of column indices of non zeros entries in X
+  //   V  nnz vector of type T non-zeros entries in X
+  //
+  template <typename T>
+  inline void find(
+    const Eigen::SparseMatrix<T>& X,
+    Eigen::Matrix<int,Eigen::Dynamic,1> & I,
+    Eigen::Matrix<int,Eigen::Dynamic,1> & J,
+    Eigen::Matrix<T,Eigen::Dynamic,1> & V);
+  // Find the non-zero entries and there respective indices in a sparse vector.
+  // Similar to matlab's [I,J,V] = find(X), but instead of [I,J] being
+  // subscripts into X, since X is a vector we just return I, a list of indices
+  // into X
+  //
+  // Templates:
+  //   T  should be a eigen sparse matrix primitive type like int or double
+  // Input:
+  //   X  vector whose entries are to be found
+  // Outputs:
+  //   I  nnz vector of indices of non zeros entries in X
+  //   V  nnz vector of type T non-zeros entries in X
+  template <typename T>
+  inline void find(
+    const Eigen::SparseVector<T>& X,
+    Eigen::Matrix<int,Eigen::Dynamic,1> & I,
+    Eigen::Matrix<T,Eigen::Dynamic,1> & V);
+}
+
+// Implementation
+#include "verbose.h"
+  
+template <typename T>
+inline void igl::find(
+  const Eigen::SparseMatrix<T>& X,
+  Eigen::Matrix<int,Eigen::Dynamic,1> & I,
+  Eigen::Matrix<int,Eigen::Dynamic,1> & J,
+  Eigen::Matrix<T,Eigen::Dynamic,1> & V)
+{
+  // Resize outputs to fit nonzeros
+  I.resize(X.nonZeros());
+  J.resize(X.nonZeros());
+  V.resize(X.nonZeros());
+
+  int i = 0;
+  // 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)
+    {
+      V(i) = it.value();
+      I(i) = it.row();
+      J(i) = it.col();
+      i++;
+    }
+  }
+}
+  
+template <typename T>
+inline void igl::find(
+  const Eigen::SparseVector<T>& X,
+  Eigen::Matrix<int,Eigen::Dynamic,1> & I,
+  Eigen::Matrix<T,Eigen::Dynamic,1> & V)
+{
+  // Resize outputs to fit nonzeros
+  I.resize(X.nonZeros());
+  V.resize(X.nonZeros());
+
+  int i = 0;
+  // loop over non-zeros
+  for(typename SparseVector<T>::InnerIterator it(X); it; ++it)
+  {
+    I(i) = it.index();
+    V(i) = it.value();
+    i++;
+  }
+}
+
+#endif

matlab-to-eigen.html

@@ -78,6 +78,41 @@ tr.d1 td
         <td></td>
       </tr>
 
+      <tr class=d0>
+        <td><pre><code>[I,J,V] = find(X)</code></pre></td>
+        <td><pre><code>igl::find(X,I,J,V)</code></pre></td>
+        <td>Matlab supports finding subscripts (I and J) as well as indices
+        (just I), but so far igl::find only supports subscripts. Also,
+        igl::find requires X to be sparse.</td>
+      </tr>
+
+      <tr class=d1>
+        <td><pre><code>X(:,j) = X(:,j) + x</code></pre></td>
+        <td><pre><code>X.col(j).array() += x</code></pre></td>
+        <td></td>
+      </tr>
+
+      <tr class=d0>
+        <td><pre><code>Adim_sum = sum(A,dim)</code></pre></td>
+        <td><pre><code>igl::sum(A,dim,Adim_sum)</code></pre></td>
+        <td>Currently the igl version only supports sparse matrix input (and
+            dim must be 1 or 2)</td>
+      </tr>
+
+      <tr class=d1>
+        <td><pre><code>D = diag(M)</code></pre></td>
+        <td><pre><code>igl::diag(M,D)</code></pre></td>
+        <td>Extract the main diagonal of a matrix. Currently igl version
+        supports sparse only.</td>
+      </tr>
+
+      <tr class=d0>
+        <td><pre><code>M = diag(D)</code></pre></td>
+        <td><pre><code>igl::diag(D,M)</code></pre></td>
+        <td>Construct new square matrix M with entries of vector D along the
+        diagonal. Currently igl version supports sparse only.</td>
+      </tr>
+
       <!-- Insert rows for each command pair -->
 
       <!-- Leave this here for copy and pasting

print_ijv.h

@@ -0,0 +1,50 @@
+#ifndef IGL_PRINT_IJV_H
+#define IGL_PRINT_IJV_H
+#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+namespace igl
+{
+  // Prints a 3 column matrix representing [I,J,V] = find(X). That is, each
+  // row is the row index, column index and value for each non zero entry. Each
+  // row is printed on a new line
+  //
+  // Templates:
+  //   T  should be a eigen sparse matrix primitive type like int or double
+  // Input:
+  //   X  m by n matrix whose entries are to be sorted
+  //   offset  optional offset for I and J indices {0}
+  template <typename T>
+  inline void print_ijv(
+    const Eigen::SparseMatrix<T>& X, 
+    const int offset=0);
+}
+
+// Implementation
+#include "find.h"
+#include <iostream>
+
+template <typename T>
+inline void igl::print_ijv(
+  const Eigen::SparseMatrix<T>& X,
+  const int offset)
+{
+  Eigen::Matrix<int,Eigen::Dynamic,1> I;
+  Eigen::Matrix<int,Eigen::Dynamic,1> J;
+  Eigen::Matrix<T,Eigen::Dynamic,1> V;
+  igl::find(X,I,J,V);
+  // Concatenate I,J,V
+  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> IJV(I.size(),3);
+  IJV.col(0) = I.cast<T>();
+  IJV.col(1) = J.cast<T>();
+  IJV.col(2) = V;
+  // Offset
+  if(offset != 0)
+  {
+    IJV.col(0).array() += offset;
+    IJV.col(1).array() += offset;
+  }
+  std::cout<<IJV;
+}
+
+#endif

sort.h

@@ -2,6 +2,7 @@
 #define IGL_SORT_H
 
 #include <vector>
+#include <Eigen/Core>
 namespace igl
 {
 

sum.h

@@ -0,0 +1,64 @@
+#ifndef IGL_SUM_H
+#define IGL_SUM_H
+#include <Eigen/Sparse>
+
+namespace igl
+{
+  // Ideally, this becomes a super overloaded function that works with sparse
+  // and dense matrices like the matlab sum function
+
+  // Sum the columns or rows of a sparse matrix
+  // Templates:
+  //   T  should be a eigen sparse matrix primitive type like int or double
+  // Inputs:
+  //   X  m by n sparse matrix
+  //   dim  dimension along which to sum (1 or 2)
+  // Output:
+  //   S  n-long sparse vector (if dim == 1) 
+  //   or
+  //   S  m-long sparse vector (if dim == 2)
+  template <typename T>
+  void sum(
+    const Eigen::SparseMatrix<T>& X, 
+    const int dim,
+    Eigen::SparseVector<T>& S);
+}
+
+template <typename T>
+void igl::sum(
+  const Eigen::SparseMatrix<T>& X, 
+  const int dim,
+  Eigen::SparseVector<T>& S)
+{
+  // dim must be 2 or 1
+  assert(dim == 1 || dim == 2);
+  // Get size of input
+  int m = X.rows();
+  int n = X.cols();
+  // resize output
+  if(dim==1)
+  {
+    S = Eigen::SparseVector<T>(n);
+  }else
+  {
+    S = Eigen::SparseVector<T>(m);
+  }
+
+  // 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(dim == 1)
+      {
+        S.coeffRef(it.col()) += it.value();
+      }else
+      {
+        S.coeffRef(it.row()) += it.value();
+      }
+    }
+  }
+
+}
+#endif