WA with typedef for floats -- does not crash (smooth minimum does not produce nan's), but still issues with eigen svd

Former-commit-id: af70f2fe82ef0493eab4138f554bbfc1ee2c3797

grad.h

@@ -32,12 +32,20 @@ namespace igl {
   // i, and A is the area of triangle (i,j,k). ^R90 represent a rotation of 
   // 90 degrees
   //
-  void grad(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, const Eigen::VectorXd &X, Eigen::MatrixXd &G );
+  template <typename T>
+  inline void grad(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &V,
+                   const Eigen::MatrixXi &F,
+                   const Eigen::Matrix<T, Eigen::Dynamic, 1>&X,
+                   Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &G );
 }
 
-inline void igl::grad(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, const Eigen::VectorXd &X, Eigen::MatrixXd &G)
+template <typename T>
+inline void igl::grad(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &V,
+                 const Eigen::MatrixXi &F,
+                 const Eigen::Matrix<T, Eigen::Dynamic, 1>&X,
+                 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &G )
 {
-  G = Eigen::MatrixXd::Zero(F.rows(),3);
+  G = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>::Zero(F.rows(),3);
   for (int i = 0; i<F.rows(); ++i)
   {
     // renaming indices of vertices of triangles for convenience
@@ -46,28 +54,28 @@ inline void igl::grad(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, const
     int i3 = F(i,2);
     
     // #F x 3 matrices of triangle edge vectors, named after opposite vertices
-    Eigen::RowVector3d v32 = V.row(i3) - V.row(i2);
-    Eigen::RowVector3d v13 = V.row(i1) - V.row(i3);
-    Eigen::RowVector3d v21 = V.row(i2) - V.row(i1);
+    Eigen::Matrix<T, 1, 3> v32 = V.row(i3) - V.row(i2);
+    Eigen::Matrix<T, 1, 3> v13 = V.row(i1) - V.row(i3);
+    Eigen::Matrix<T, 1, 3> v21 = V.row(i2) - V.row(i1);
     
     // area of parallelogram is twice area of triangle
     // area of parallelogram is || v1 x v2 || 
-    Eigen::RowVector3d n  = v32.cross(v13); 
+    Eigen::Matrix<T, 1, 3> n  = v32.cross(v13); 
     
     // This does correct l2 norm of rows, so that it contains #F list of twice
     // triangle areas
     double dblA = std::sqrt(n.dot(n));
     
     // now normalize normals to get unit normals
-    Eigen::RowVector3d u = n / dblA;
+    Eigen::Matrix<T, 1, 3> u = n / dblA;
     
     // rotate each vector 90 degrees around normal
     double norm21 = std::sqrt(v21.dot(v21));
     double norm13 = std::sqrt(v13.dot(v13));
-    Eigen::RowVector3d eperp21 = u.cross(v21);
+    Eigen::Matrix<T, 1, 3> eperp21 = u.cross(v21);
     eperp21 = eperp21 / std::sqrt(eperp21.dot(eperp21));
     eperp21 *= norm21;
-    Eigen::RowVector3d eperp13 = u.cross(v13);
+    Eigen::Matrix<T, 1, 3> eperp13 = u.cross(v13);
     eperp13 = eperp13 / std::sqrt(eperp13.dot(eperp13));
     eperp13 *= norm13;
     

ismanifold.h

@@ -13,38 +13,39 @@
 
 namespace igl 
 {
-    // check if the mesh is edge-manifold
-    inline bool isManifold(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F)
+  // check if the mesh is edge-manifold
+  template<typename T>
+  inline bool isManifold(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& V, const Eigen::MatrixXi& F)
+  {
+    std::vector<std::vector<int> > TTT;
+    for(int f=0;f<F.rows();++f)
+      for (int i=0;i<3;++i)
+      {
+        // v1 v2 f ei 
+        int v1 = F(f,i);
+        int v2 = F(f,(i+1)%3);
+        if (v1 > v2) std::swap(v1,v2);
+        std::vector<int> r(4);
+        r[0] = v1; r[1] = v2;
+        r[2] = f;  r[3] = i;
+        TTT.push_back(r);
+      }
+    std::sort(TTT.begin(),TTT.end());
+    
+    for(int i=2;i<TTT.size();++i)
     {
-        std::vector<std::vector<int> > TTT;
-        for(int f=0;f<F.rows();++f)
-            for (int i=0;i<3;++i)
-            {
-                // v1 v2 f ei 
-                int v1 = F(f,i);
-                int v2 = F(f,(i+1)%3);
-                if (v1 > v2) std::swap(v1,v2);
-                std::vector<int> r(4);
-                r[0] = v1; r[1] = v2;
-                r[2] = f;  r[3] = i;
-                TTT.push_back(r);
-            }
-        std::sort(TTT.begin(),TTT.end());
-        
-        for(int i=2;i<TTT.size();++i)
-        {
-            std::vector<int>& r1 = TTT[i-2];
-            std::vector<int>& r2 = TTT[i-1];
-            std::vector<int>& r3 = TTT[i];
-            if ( (r1[0] == r2[0] && r2[0] == r3[0]) 
-                && 
-                (r1[1] == r2[1] && r2[1] == r3[1]) )
-            {
-                return false;
-            }
-        }
-        return true;
+      std::vector<int>& r1 = TTT[i-2];
+      std::vector<int>& r2 = TTT[i-1];
+      std::vector<int>& r3 = TTT[i];
+      if ( (r1[0] == r2[0] && r2[0] == r3[0]) 
+          && 
+          (r1[1] == r2[1] && r2[1] == r3[1]) )
+      {
+        return false;
+      }
     }
+    return true;
+  }
 }
 
 #endif

matlabinterface.h

@@ -41,12 +41,18 @@ namespace igl
   // Send a matrix to MATLAB
   inline void mlsetmatrix(Engine** engine, std::string name, const Eigen::MatrixXd& M);
   
+  // Send a matrix to MATLAB
+  inline void mlsetmatrix(Engine** engine, std::string name, const Eigen::MatrixXf& M);
+
   // Send a matrix to MATLAB
   inline void mlsetmatrix(Engine** engine, std::string name, const Eigen::MatrixXi& M);
   
   // Receive a matrix from MATLAB
   inline void mlgetmatrix(Engine** engine, std::string name, Eigen::MatrixXd& M);
   
+  // Receive a matrix from MATLAB
+  inline void mlgetmatrix(Engine** engine, std::string name, Eigen::MatrixXf& M);
+
   // Receive a matrix from MATLAB
   inline void mlgetmatrix(Engine** engine, std::string name, Eigen::MatrixXi& M);
   
@@ -95,6 +101,24 @@ inline void igl::mlsetmatrix(Engine** mlengine, std::string name, const Eigen::M
   mxDestroyArray(A);
 }
 
+// Send a matrix to MATLAB
+inline void igl::mlsetmatrix(Engine** mlengine, std::string name, const Eigen::MatrixXf& M)
+{
+  if (*mlengine == 0)
+    mlinit(mlengine);
+  
+  mxArray *A = mxCreateDoubleMatrix(M.rows(), M.cols(), mxREAL);
+  double *pM = mxGetPr(A);
+  
+  int c = 0;
+  for(int j=0; j<M.cols();++j)
+    for(int i=0; i<M.rows();++i)
+      pM[c++] = double(M(i,j));
+  
+  engPutVariable(*mlengine, name.c_str(), A);
+  mxDestroyArray(A);
+}
+
 // Send a matrix to MATLAB
 inline void igl::mlsetmatrix(Engine** mlengine, std::string name, const Eigen::MatrixXi& M)
 {
@@ -147,6 +171,39 @@ inline void igl::mlgetmatrix(Engine** mlengine, std::string name, Eigen::MatrixX
   mxDestroyArray(ary);
 }
 
+inline void igl::mlgetmatrix(Engine** mlengine, std::string name, Eigen::MatrixXf& M)
+{
+  if (*mlengine == 0)
+    mlinit(mlengine);
+  
+  unsigned long m = 0;
+  unsigned long n = 0;
+  std::vector<double> t;
+  
+  mxArray *ary = engGetVariable(*mlengine, name.c_str());
+  if (ary == NULL)
+  {
+    m = 0;
+    n = 0;
+    M = Eigen::MatrixXf(0,0);
+  }
+  else
+  {
+    m = mxGetM(ary);
+    n = mxGetN(ary);
+    M = Eigen::MatrixXf(m,n);
+    
+    double *pM = mxGetPr(ary);
+    
+    int c = 0;
+    for(int j=0; j<M.cols();++j)
+      for(int i=0; i<M.rows();++i)
+        M(i,j) = pM[c++];
+  }
+  
+  mxDestroyArray(ary);
+}
+
 // Receive a matrix from MATLAB
 inline void igl::mlgetmatrix(Engine** mlengine, std::string name, Eigen::MatrixXi& M)
 {

moveFV.h

@@ -20,14 +20,22 @@ namespace igl
   // 
   // Output:
   // SV: scalar field defined on vertices
-  void moveFV(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, const Eigen::MatrixXd &S, Eigen::MatrixXd &SV);
+  template <typename T>
+  inline void moveFV(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &V,
+              const Eigen::MatrixXi &F,
+              const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &S,
+              Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &SV);
 }
 
-inline void igl::moveFV(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, const Eigen::MatrixXd &S, Eigen::MatrixXd &SV)
+template <typename T>
+inline void igl::moveFV(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &V,
+            const Eigen::MatrixXi &F,
+            const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &S,
+            Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &SV)
 {
   
-  SV = Eigen::MatrixXd::Zero(V.rows(),S.cols());
-  Eigen::VectorXd COUNT = Eigen::VectorXd::Zero(V.rows());
+  SV = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>::Zero(V.rows(),S.cols());
+  Eigen::Matrix<T, Eigen::Dynamic, 1> COUNT = Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(V.rows());
   for (int i = 0; i <F.rows(); ++i)
   {
     for (int j = 0; j<F.cols(); ++j)

removeDuplicates.h

@@ -22,12 +22,14 @@ namespace igl
   // Output:
   // NV, NF: new mesh without duplicate vertices
   
-  void removeDuplicates(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, Eigen::MatrixXd &NV, Eigen::MatrixXi &NF, Eigen::VectorXi &I, const double epsilon);
+  template <typename T>
+  inline void removeDuplicates(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &V, const Eigen::MatrixXi &F, Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &NV, Eigen::MatrixXi &NF, Eigen::VectorXi &I, const double epsilon = 2.2204e-15);
   
 }
 
 // Implementation
-inline void igl::removeDuplicates(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, Eigen::MatrixXd &NV, Eigen::MatrixXi &NF, Eigen::VectorXi &I, const double epsilon = 2.2204e-15)
+template <typename T>
+inline void igl::removeDuplicates(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &V, const Eigen::MatrixXi &F, Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &NV, Eigen::MatrixXi &NF, Eigen::VectorXi &I, const double epsilon)
 {
   assert (F.cols() == 3);
   
@@ -41,7 +43,7 @@ inline void igl::removeDuplicates(const Eigen::MatrixXd &V, const Eigen::MatrixX
   for (int i =0; i <n; ++i)
     VISITED[i] = false;
   
-  NV = Eigen::MatrixXd(n,V.cols());
+  NV = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>(n,V.cols());
   int count = 0;
   Eigen::VectorXd d(n);
   for (int i =0; i <n; ++i)

removeUnreferenced.h

@@ -19,12 +19,14 @@ namespace igl
   // Output:
   // NV, NF: new mesh without unreferenced vertices
   
-  void removeUnreferenced(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, Eigen::MatrixXd &NV, Eigen::MatrixXi &NF, Eigen::VectorXi &I);
+  template <typename T>
+  inline void removeUnreferenced(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &V, const Eigen::MatrixXi &F, Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &NV, Eigen::MatrixXi &NF, Eigen::VectorXi &I);
   
 }
 
 // Implementation
-inline void igl::removeUnreferenced(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, Eigen::MatrixXd &NV, Eigen::MatrixXi &NF, Eigen::VectorXi &I)
+template <typename T>
+inline void igl::removeUnreferenced(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &V, const Eigen::MatrixXi &F, Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &NV, Eigen::MatrixXi &NF, Eigen::VectorXi &I)
 {
 
   // Mark referenced vertices
@@ -42,7 +44,7 @@ inline void igl::removeUnreferenced(const Eigen::MatrixXd &V, const Eigen::Matri
   // Sum the occupied cells 
   int newsize = mark.sum();
   
-  NV = Eigen::MatrixXd(newsize,V.cols());
+  NV = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>(newsize,V.cols());
   NF = Eigen::MatrixXi(F.rows(),F.cols());
   I  = Eigen::MatrixXi(V.rows(),1);
   

tt.h

@@ -16,7 +16,8 @@ namespace igl
 {
 
   // Preprocessing
-  inline void tt_preprocess(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F, std::vector<std::vector<int> >& TTT)
+  template<typename T> 
+  inline void tt_preprocess(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& V, const Eigen::MatrixXi& F, std::vector<std::vector<int> >& TTT)
   {
     for(int f=0;f<F.rows();++f)
       for (int i=0;i<3;++i)
@@ -68,22 +69,24 @@ namespace igl
   }
 
   // Compute triangle-triangle adjacency
-  inline void tt(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F, Eigen::MatrixXi& TT)
+  template<typename T> 
+  inline void tt(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& V, const Eigen::MatrixXi& F, Eigen::MatrixXi& TT)
   {
     assert(igl::isManifold(V,F));
     std::vector<std::vector<int> > TTT;
     
-    tt_preprocess(V,F,TTT);
+    tt_preprocess<T>(V,F,TTT);
     tt_extractTT(F,TTT,TT);
   }
 
   // Compute triangle-triangle adjacency with indices
-  inline void tt(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F, Eigen::MatrixXi& TT, Eigen::MatrixXi& TTi)
+  template<typename T> 
+  inline void tt(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& V, const Eigen::MatrixXi& F, Eigen::MatrixXi& TT, Eigen::MatrixXi& TTi)
   {
     assert(igl::isManifold(V,F));
     std::vector<std::vector<int> > TTT;
     
-    tt_preprocess(V,F,TTT);
+    tt_preprocess<T>(V,F,TTT);
     tt_extractTT(F,TTT,TT);
     tt_extractTTi(F,TTT,TTi);
   }

vf.h

@@ -16,9 +16,9 @@ namespace igl
   //
   // See also: edges, cotmatrix, diag, vv
     
-  template <typename T>
+  template <typename T, typename S>
   inline void vf( 
-    const Eigen::MatrixXd & V, 
+    const Eigen::Matrix<S, Eigen::Dynamic, Eigen::Dynamic> & V, 
     const Eigen::MatrixXi & F, 
     vector<vector<T> >& VF, vector<vector<T> >& VFi);
 }
@@ -26,9 +26,9 @@ namespace igl
 // Implementation
 #include "verbose.h"
 
-template <typename T>
+template <typename T, typename S>
 inline void igl::vf(
-  const Eigen::MatrixXd & V, 
+  const Eigen::Matrix<S, Eigen::Dynamic, Eigen::Dynamic> & V, 
   const Eigen::MatrixXi & F, 
   vector<vector<T> >& VF, vector<vector<T> >& VFi)
 {