colon, full and speye functions to match matlab

Former-commit-id: 3ac845befd37976bceefbc53fd58c290080bc95f

colon.h

@@ -0,0 +1,112 @@
+#ifndef IGL_COLON_H
+#define IGL_COLON_H
+#include <Eigen/Dense>
+namespace igl
+{
+  // Colon operator like matlab's colon operator. Enumerats values between low
+  // and hi with step step.
+  // Templates:
+  //   L  should be a eigen matrix primitive type like int or double
+  //   S  should be a eigen matrix primitive type like int or double
+  //   H  should be a eigen matrix primitive type like int or double
+  //   T  should be a eigen matrix primitive type like int or double
+  // Inputs:
+  //   low  starting value if step is valid then this is *always* the first
+  //     element of I
+  //   step  step difference between sequential elements returned in I,
+  //     remember this will be cast to template T at compile time. If low<hi
+  //     then step must be positive. If low>hi then step must be negative.
+  //     Otherwise I will be set to empty.
+  //   hi  ending value, if (hi-low)%step is zero then this will be the last
+  //     element in I. If step is positive there will be no elements greater
+  //     than hi, vice versa if hi<low
+  // Output:
+  //   I  list of values from low to hi with step size step
+  template <typename L,typename S,typename H,typename T>
+  inline void colon(
+    const L low, 
+    const S step, 
+    const H hi, 
+    Eigen::Matrix<T,Eigen::Dynamic,1> & I);
+  // Same as above but step == (T)1
+  template <typename L,typename H,typename T>
+  inline void colon(
+    const L low, 
+    const H hi, 
+    Eigen::Matrix<T,Eigen::Dynamic,1> & I);
+  // Return output rather than set in reference
+  template <typename T,typename L,typename H>
+  inline Eigen::Matrix<T,Eigen::Dynamic,1> colon(
+    const L low, 
+    const H hi);
+}
+
+// Implementation
+#include <cstdio>
+
+template <typename L,typename S,typename H,typename T>
+inline void igl::colon(
+  const L low, 
+  const S step, 
+  const H hi, 
+  Eigen::Matrix<T,Eigen::Dynamic,1> & I)
+{
+  if(low < hi)
+  {
+    if(step < 0)
+    {
+      I.resize(0);
+      fprintf(stderr,"Error: colon() low(%g)<hi(%g) but step(%g)<0\n",
+        (double)low,
+        (double)hi,
+        (double)step);
+      return;
+    }
+  }
+  if(low > hi)
+  {
+    if(step > 0)
+    {
+      I.resize(0);
+      fprintf(stderr,"Error: colon() low(%g)>hi(%g) but step(%g)<0\n",
+        (double)low,
+        (double)hi,
+        (double)step);
+      return;
+    }
+  }
+  // resize output
+  int n = floor((hi-low)/step)+1;
+  I.resize(n);
+  int i = 0;
+  T v = (T)low;
+  while((low<hi && (H)v<=hi) || (low>hi && (H)v>=hi))
+  {
+    I(i) = v;
+    v = v + (T)step;
+    i++;
+  }
+  assert(i==n);
+}
+
+template <typename L,typename H,typename T>
+inline void igl::colon(
+  const L low, 
+  const H hi, 
+  Eigen::Matrix<T,Eigen::Dynamic,1> & I)
+{
+  return igl::colon(low,(T)1,hi,I);
+}
+
+template <typename T,typename L,typename H>
+inline Eigen::Matrix<T,Eigen::Dynamic,1> igl::colon(
+  const L low, 
+  const H hi)
+{
+  Eigen::Matrix<T,Eigen::Dynamic,1> I;
+  igl::colon(low,hi,I);
+  return I;
+}
+
+#endif
+

full.h

@@ -0,0 +1,39 @@
+#ifndef IGL_FULL_H
+#define IGL_FULL_H
+#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+namespace igl
+{
+  // Convert a sparsematrix into a full one
+  // Templates:
+  //   T  should be a eigen sparse matrix primitive type like int or double
+  // Input:
+  //   A  m by n sparse matrix
+  // Output:
+  //   B  m by n dense/full matrix
+  template <typename T>
+  inline void full(
+    const Eigen::SparseMatrix<T> & A,
+    Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & B);
+}
+
+// Implementation
+
+template <typename T>
+inline void igl::full(
+  const Eigen::SparseMatrix<T> & A,
+  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & B)
+{
+  B = Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>::Zero(A.rows(),A.cols());
+  // Iterate over outside
+  for(int k=0; k<A.outerSize(); ++k)
+  {
+    // Iterate over inside
+    for(typename Eigen::SparseMatrix<T>::InnerIterator it (A,k); it; ++it)
+    {
+      B(it.row(),it.col()) = it.value();
+    }
+  }
+}
+#endif

speye.h

@@ -0,0 +1,50 @@
+#ifndef IGL_SPEYE_H
+#define IGL_SPEYE_H
+
+#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
+#include <Eigen/Sparse>
+
+namespace igl
+{
+  // Builds an m by n sparse identity matrix like matlab's speye function
+  // Templates:
+  //   T  should be a eigen sparse matrix primitive type like int or double
+  // Inputs:
+  //   m  number of rows
+  //   n  number of cols
+  // Outputs:
+  //   I  m by n sparse matrix with 1's on the main diagonal
+  template <typename T>
+  inline void speye(const int n,const int m, Eigen::SparseMatrix<T> & I);
+  // Builds an n by n sparse identity matrix like matlab's speye function
+  // Templates:
+  //   T  should be a eigen sparse matrix primitive type like int or double
+  // Inputs:
+  //   n  number of rows and cols
+  // Outputs:
+  //   I  n by n sparse matrix with 1's on the main diagonal
+  template <typename T>
+  inline void speye(const int n, Eigen::SparseMatrix<T> & I);
+}
+
+template <typename T>
+inline void igl::speye(const int m, const int n, Eigen::SparseMatrix<T> & I)
+{
+  // size of diagonal
+  int d = (m<n?m:n);
+  I = Eigen::SparseMatrix<T>(m,n);
+  I.reserve(d);
+  for(int i = 0;i<d;i++)
+  {
+    I.insert(i,i) = 1.0;
+  }
+  I.finalize();
+}
+
+template <typename T>
+inline void igl::speye(const int n, Eigen::SparseMatrix<T> & I)
+{
+  return igl::speye(n,n,I);
+}
+
+#endif