Fixed and generalized cotmatrix, write .mesh file, cat like matlab's cat

Former-commit-id: 7e5b490a5864fda38e071b106db729b48f607b46

cat.h

@@ -0,0 +1,193 @@
+#ifndef IGL_CAT_H
+#define IGL_CAT_H
+
+#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
+#include <Eigen/Sparse>
+#include <Eigen/Dense>
+
+namespace igl
+{
+  // If you're using Dense matrices you might be better off using the << operator
+
+  // This is an attempt to act like matlab's cat function.
+
+  // Perform concatenation of a two matrices along a single dimension
+  // If dim == 1, then C = [A;B]. If dim == 2 then C = [A B]
+  // 
+  // Template:
+  //   Scalar  scalar data type for sparse matrices like double or int
+  //   Mat  matrix type for all matrices (e.g. MatrixXd, SparseMatrix)
+  //   MatC  matrix type for ouput matrix (e.g. MatrixXd) needs to support
+  //     resize
+  // Inputs:
+  //   A  first input matrix
+  //   B  second input matrix
+  //   dim  dimension along which to concatenate, 0 or 1
+  // Outputs:
+  //   C  output matrix
+  //   
+  template <typename Scalar>
+  void cat(
+      const int dim, 
+      const Eigen::SparseMatrix<Scalar> & A, 
+      const Eigen::SparseMatrix<Scalar> & B, 
+      Eigen::SparseMatrix<Scalar> & C);
+  template <typename Derived, class MatC>
+  void cat(
+    const int dim,
+    const Eigen::MatrixBase<Derived> & A, 
+    const Eigen::MatrixBase<Derived> & B,
+    MatC & C);
+  // Wrapper that returns C
+  template <class Mat>
+  Mat cat(const int dim, const Mat & A, const Mat & B);
+
+  // Note: Maybe we can autogenerate a bunch of overloads D = cat(int,A,B,C),
+  // E = cat(int,A,B,C,D), etc. 
+
+  // Concatenate a "matrix" of blocks
+  // C = [A0;A1;A2;...;An] where Ai = [A[i][0] A[i][1] ... A[i][m]];
+  //
+  // Inputs:
+  //   A  a matrix (vector of row vectors)
+  // Output:
+  //   C
+  template <class Mat>
+    void cat(const std::vector<std::vector< Mat > > & A, Mat & C);
+}
+
+// Implementation
+
+// Sparse matrices need to be handled carefully. Because C++ does not 
+// Template:
+//   Scalar  sparse matrix scalar type, e.g. double
+template <typename Scalar>
+void igl::cat(
+    const int dim, 
+    const Eigen::SparseMatrix<Scalar> & A, 
+    const Eigen::SparseMatrix<Scalar> & B, 
+    Eigen::SparseMatrix<Scalar> & C)
+{
+  assert(dim == 1 || dim == 2);
+  using namespace Eigen;
+  // Special case if B or A is empty
+  if(A.size() == 0)
+  {
+    C = B;
+    return;
+  }
+  if(B.size() == 0)
+  {
+    C = A;
+    return;
+  }
+
+  DynamicSparseMatrix<Scalar, RowMajor> dyn_C;
+  if(dim == 1)
+  {
+    assert(A.cols() == B.cols());
+    dyn_C.resize(A.rows()+B.rows(),A.cols());
+  }else if(dim == 2)
+  {
+    assert(A.rows() == B.rows());
+    dyn_C.resize(A.rows(),A.cols()+B.cols());
+  }else
+  {
+    fprintf(stderr,"cat.h: Error: Unsupported dimension %d\n",dim);
+  }
+
+  dyn_C.reserve(A.nonZeros()+B.nonZeros());
+
+  // Iterate over outside of A
+  for(int k=0; k<A.outerSize(); ++k)
+  {
+    // Iterate over inside
+    for(typename SparseMatrix<Scalar>::InnerIterator it (A,k); it; ++it)
+    {
+      dyn_C.coeffRef(it.row(),it.col()) += it.value();
+    }
+  }
+
+  // Iterate over outside of B
+  for(int k=0; k<B.outerSize(); ++k)
+  {
+    // Iterate over inside
+    for(typename SparseMatrix<Scalar>::InnerIterator it (B,k); it; ++it)
+    {
+      int r = (dim == 1 ? A.rows()+it.row() : it.row());
+      int c = (dim == 2 ? A.cols()+it.col() : it.col());
+      dyn_C.coeffRef(r,c) += it.value();
+    }
+  }
+
+  C = SparseMatrix<Scalar>(dyn_C);
+}
+
+template <typename Derived, class MatC>
+void igl::cat(
+  const int dim,
+  const Eigen::MatrixBase<Derived> & A, 
+  const Eigen::MatrixBase<Derived> & B,
+  MatC & C)
+{
+  assert(dim == 1 || dim == 2);
+  // Special case if B or A is empty
+  if(A.size() == 0)
+  {
+    C = B;
+    return;
+  }
+  if(B.size() == 0)
+  {
+    C = A;
+    return;
+  }
+
+  if(dim == 1)
+  {
+    assert(A.cols() == B.cols());
+    C.resize(A.rows()+B.rows(),A.cols());
+    C << A,B;
+  }else if(dim == 2)
+  {
+    assert(A.rows() == B.rows());
+    C.resize(A.rows(),A.cols()+B.cols());
+    C << A,B;
+  }else
+  {
+    fprintf(stderr,"cat.h: Error: Unsupported dimension %d\n",dim);
+  }
+}
+
+template <class Mat>
+Mat igl::cat(const int dim, const Mat & A, const Mat & B)
+{
+  assert(dim == 1 || dim == 2);
+  Mat C;
+  igl::cat(dim,A,B,C);
+  return C;
+}
+
+template <class Mat>
+void cat(const std::vector<std::vector< Mat > > & A, Mat & C)
+{
+  assert(dim == 1 || dim == 2);
+  using namespace igl;
+  using namespace std;
+  // Start with empty matrix
+  C.resize(0,0);
+  for(typename vector<vector< Mat > >::const_iterator rit = A.begin(); rit != A.end(); rit++)
+  {
+    // Concatenate each row horizontally
+    // Start with empty matrix
+    Mat row(0,0);
+    for(typename vector<vector< Mat > >::iterator cit = A.begin(); rit != A.end(); rit++)
+    {
+      row = cat(2,row,*cit);
+    }
+    // Concatenate rows vertically
+    C = cat(1,C,row);
+  }
+}
+
+#endif

cotmatrix.h

@@ -16,119 +16,95 @@ namespace igl
 {
   // Constructs the cotangent stiffness matrix (discrete laplacian) 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 3 list of mesh faces (must be triangles)
+  //   F  #F by simplex_size 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
+  //
+  // Known bugs: off by 1e-16 on regular grid. I think its a problem of
+  // arithmetic order in cotangent.h: C(i,e) = (arithmetic)/dblA/4
+  template <typename DerivedV, typename DerivedF, typename Scalar>
   inline 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
-  //   
-  inline void computeCotWeights(
-    const Eigen::Vector3d& v1, 
-    const Eigen::Vector3d& v2, 
-    const Eigen::Vector3d& v3, 
-    double& cot1, 
-    double& cot2, 
-    double& cot3);
+    const Eigen::MatrixBase<DerivedV> & V, 
+    const Eigen::MatrixBase<DerivedF> & F, 
+    Eigen::SparseMatrix<Scalar>& L);
 }
 
 // Implementation
 
 // For error printing
 #include <cstdio>
+#include "cotangent.h"
 
-inline 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;
-}
-
+template <typename DerivedV, typename DerivedF, typename Scalar>
 inline void igl::cotmatrix(
-  const Eigen::MatrixXd & V, 
-  const Eigen::MatrixXi & F, 
-  Eigen::SparseMatrix<double>& L)
+  const Eigen::MatrixBase<DerivedV> & V, 
+  const Eigen::MatrixBase<DerivedF> & F, 
+  Eigen::SparseMatrix<Scalar>& L)
 {
-  // Assumes vertices are 3D or 2D
-  assert((V.cols() == 3) || (V.cols() == 2));
-  int dim = V.cols();
-  // Assumes faces are triangles
-  assert(F.cols() == 3);
+  using namespace igl;
+  using namespace Eigen;
 
-  Eigen::DynamicSparseMatrix<double, Eigen::RowMajor> dyn_L (V.rows(), V.rows());
-  // This is important! it could decrease the comptuation time by a factor of 2
-  // Laplacian for a closed 2d manifold mesh will have on average 7 entries per
-  // row
-  dyn_L.reserve(7*V.rows());
-  
-  // Loop over triangles
-  for (unsigned i = 0; i < F.rows(); i++)
+  DynamicSparseMatrix<Scalar, RowMajor> dyn_L (V.rows(), V.rows());
+  Matrix<int,Dynamic,2> edges;
+  int simplex_size = F.cols();
+  // 3 for triangles, 4 for tets
+  assert(simplex_size == 3 || simplex_size == 4);
+  if(simplex_size == 3)
   {
-    // 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), (dim==2?0:V(vi1,2)));
-    Eigen::Vector3d v2(V(vi2,0), V(vi2,1), (dim==2?0:V(vi2,2)));
-    Eigen::Vector3d v3(V(vi3,0), V(vi3,1), (dim==2?0: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;
+    // This is important! it could decrease the comptuation time by a factor of 2
+    // Laplacian for a closed 2d manifold mesh will have on average 7 entries per
+    // row
+    dyn_L.reserve(7*V.rows());
+    edges.resize(3,2);
+    edges << 
+      1,2,
+      2,0,
+      0,1;
+  }else if(simplex_size == 4)
+  {
+    dyn_L.reserve(17*V.rows());
+    edges.resize(6,2);
+    edges << 
+      1,2,
+      2,0,
+      0,1,
+      3,0,
+      3,1,
+      3,2;
+  }else
+  {
+    return;
   }
-
-  for (int k=0; k < dyn_L.outerSize(); ++k)
+  // Gather cotangents
+  Matrix<Scalar,Dynamic,Dynamic> C;
+  cotangent(V,F,C);
+  
+  // Loop over triangles
+  for(int i = 0; i < F.rows(); i++)
   {
-    double tmp = 0.0f;
-    for(Eigen::DynamicSparseMatrix<double, Eigen::RowMajor>::InnerIterator it (dyn_L, k); it; ++it)
+    // loop over edges of element
+    for(int e = 0;e<edges.rows();e++)
     {
-      tmp += it.value();
+      int source = F(i,edges(e,0));
+      int dest = F(i,edges(e,1));
+      dyn_L.coeffRef(source,dest) += C(i,e);
+      dyn_L.coeffRef(dest,source) += C(i,e);
+      dyn_L.coeffRef(source,source) += -C(i,e);
+      dyn_L.coeffRef(dest,dest) += -C(i,e);
     }
-    dyn_L.coeffRef(k,k) = -tmp;
   }
-  L = Eigen::SparseMatrix<double>(dyn_L);
+    // Corner indices of this triangle
+  L = SparseMatrix<Scalar>(dyn_L);
 }
-
 #endif

matlab-to-eigen.html

@@ -227,6 +227,12 @@ tr.gotcha2 td
         <td></td>
       </tr>
 
+      <tr class=d1>
+        <td><pre><code>M = mode(X,dim)</code></pre></td>
+        <td><pre><code>igl::mode(X,dim,M)</code></pre></td>
+        <td></td>
+      </tr>
+
       <!-- Insert rows for each command pair -->
 
       <!-- Leave this here for copy and pasting

mode.h

@@ -0,0 +1,69 @@
+#ifndef IGL_MODE_H
+#define IGL_MODE_H
+#include <Eigen/Dense>
+namespace igl
+{
+  // Takes mode of coefficients in a matrix along a given dension
+  //
+  // Templates:
+  //   T  should be a eigen matrix primitive type like int or double
+  // Inputs:
+  //   X  m by n original matrix
+  //   d  dension along which to take mode, m or n
+  // Outputs:
+  //   M  vector containing mode along dension d, if d==1 then this will be a
+  //     n-long vector if d==2 then this will be a m-long vector
+  template <typename T>
+  void mode(
+    const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & X,
+    const int d, 
+    Eigen::Matrix<T,Eigen::Dynamic,1> & M);
+}
+
+// Implementation 
+#include <vector>
+
+template <typename T>
+void igl::mode(
+  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & X,
+  const int d, 
+  Eigen::Matrix<T,Eigen::Dynamic,1> & M)
+{
+  assert(d==1 || d==2);
+  using namespace std;
+  int m = X.rows();
+  int n = X.cols();
+  M.resize((d==1)?n:m,1);
+  for(int i = 0;i<((d==2)?m:n);i++)
+  {
+    vector<int> counts(((d==2)?n:m),0);
+    for(int j = 0;j<((d==2)?n:m);j++)
+    {
+      T v = (d==2)?X(i,j):X(j,i);
+      for(int k = 0;k<((d==2)?n:m);k++)
+      {
+        T u = (d==2)?X(i,k):X(k,i);
+        if(v == u)
+        {
+          counts[k]++;
+        }
+      }
+    }
+    assert(counts.size() > 0);
+    int max_count = -1;
+    int max_count_j = -1;
+    int j =0;
+    for(vector<int>::iterator it = counts.begin();it<counts.end();it++)
+    {
+      if(max_count < *it)
+      {
+        max_count = *it;
+        max_count_j = j;
+      }
+      j++;
+    }
+    M(i,0) = (d==2)?X(i,max_count_j):X(max_count_j,i);
+  }
+}
+
+#endif

repdiag.h

@@ -34,6 +34,9 @@ namespace igl
     const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & A,
     const int d,
     Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & B);
+  // Wrapper with B as output
+  template <class Mat>
+  inline Mat repdiag(const Mat & A, const int d);
 }
 
 // Implementation
@@ -88,4 +91,12 @@ inline void igl::repdiag(
   }
 }
 
+// Wrapper with B as output
+template <class Mat>
+inline Mat igl::repdiag(const Mat & A, const int d)
+{
+  Mat B;
+  repdiag(A,d,B);
+  return B;
+}
 #endif

writeMESH.h

@@ -0,0 +1,122 @@
+#ifndef IGL_WRITEMESH_H
+#define IGL_WRITEMESH_H
+
+#include <string>
+#include <vector>
+
+namespace igl
+{
+  // save a tetrahedral volume mesh to a .mesh file
+  //
+  // Templates:
+  //   Scalar  type for positions and vectors (will be cast as double)
+  //   Index  type for indices (will be cast to int)
+  // Input:
+  //   mesh_file_name  path of .mesh file
+  // Outputs:
+  //   V  double matrix of vertex positions  #V by 3
+  //   T  #T list of tet indices into vertex positions
+  //   F  #F list of face indices into vertex positions
+  template <typename Scalar, typename Index>
+  inline bool writeMESH(
+    const std::string mesh_file_name,
+    std::vector<std::vector<Scalar > > & V,
+    std::vector<std::vector<Index > > & T,
+    std::vector<std::vector<Index > > & F);
+
+  // Input:
+  //   mesh_file_name  path of .mesh file
+  // Outputs:
+  //   V  eigen double matrix #V by 3
+  //   T  eigen int matrix #T by 4
+  //   F  eigen int matrix #F by 3
+  inline bool writeMESH(
+    const std::string str,
+    Eigen::MatrixXd& V,
+    Eigen::MatrixXi& T,
+    Eigen::MatrixXi& F);
+}
+
+// Implementation
+#include <cstdio>
+#include "verbose.h"
+
+template <typename Scalar, typename Index>
+inline bool igl::writeMESH(
+  const std::string mesh_file_name,
+  std::vector<std::vector<Scalar > > & V,
+  std::vector<std::vector<Index > > & T,
+  std::vector<std::vector<Index > > & F)
+{
+  // not implemented but should be
+  assert(false);
+}
+
+#include <Eigen/Core>
+
+inline bool igl::writeMESH(
+  const std::string str,
+  Eigen::MatrixXd& V,
+  Eigen::MatrixXi& T,
+  Eigen::MatrixXi& F)
+{
+  using namespace std;
+  using namespace igl;
+  using namespace Eigen;
+
+  FILE * mesh_file = fopen(str.c_str(),"w");
+  if(NULL==mesh_file)
+  {
+    fprintf(stderr,"IOError: %s could not be opened...",str.c_str());
+    return false;
+  }
+  // print header
+  fprintf(mesh_file,"MeshVersionFormatted 1\n");
+  fprintf(mesh_file,"Dimension 3\n");
+  // print tet vertices
+  fprintf(mesh_file,"Vertices\n");
+  // print number of tet vertices
+  size_t number_of_tet_vertices = V.rows();
+  fprintf(mesh_file,"%ld\n",number_of_tet_vertices);
+  // loop over tet vertices
+  for(size_t i = 0;i<number_of_tet_vertices;i++)
+  {
+    // print position of ith tet vertex
+    fprintf(mesh_file,"%lg %lg %lg 1\n",
+      (double)V(i,0),
+      (double)V(i,1),
+      (double)V(i,2));
+  }
+  verbose("WARNING: save_mesh() assumes that vertices have"
+      " same indices in surface as volume...\n");
+  // print faces
+  fprintf(mesh_file,"Triangles\n");
+  // print number of triangles
+  size_t number_of_triangles = F.rows();
+  fprintf(mesh_file,"%ld\n",number_of_triangles);
+  // loop over faces
+  for(int i = 0;i<number_of_triangles;i++)
+  {
+    // loop over vertices in face
+    fprintf(mesh_file,"%ld %ld %ld 1\n", F(i,0)+1, F(i,1)+1, F(i,2)+1);
+  }
+  // print tetrahedra
+  fprintf(mesh_file,"Tetrahedra\n");
+  size_t number_of_tetrahedra = T.rows();
+  // print number of tetrahedra
+  fprintf(mesh_file,"%ld\n",number_of_tetrahedra);
+  // loop over tetrahedra
+  for(size_t i = 0; i < number_of_tetrahedra;i++)
+  {
+    // mesh standard uses 1-based indexing
+    fprintf(mesh_file, "%ld %ld %ld %ld 1\n",
+      T(i,0)+1,
+      T(i,1)+1,
+      T(i,2)+1,
+      T(i,3)+1);
+  }
+  fclose(mesh_file);
+  return true;
+}
+
+#endif