Merge branch 'rayform-master' into alecjacobson

Former-commit-id: dfd5e565b3466218a00743215ff196f649382adb

include/igl/AABB.cpp

@@ -391,7 +391,7 @@ igl::AABB<DerivedV,DIM>::squared_distance(
       int i_left;
       RowVectorDIMS c_left = c;
       Scalar sqr_d_left = m_left->squared_distance(V,Ele,p,sqr_d,i_left,c_left);
-      set_min(p,sqr_d_left,i_left,c_left,sqr_d,i,c);
+      this->set_min(p,sqr_d_left,i_left,c_left,sqr_d,i,c);
       looked_left = true;
     };
     const auto & look_right = [&]()
@@ -399,7 +399,7 @@ igl::AABB<DerivedV,DIM>::squared_distance(
       int i_right;
       RowVectorDIMS c_right = c;
       Scalar sqr_d_right = m_right->squared_distance(V,Ele,p,sqr_d,i_right,c_right);
-      set_min(p,sqr_d_right,i_right,c_right,sqr_d,i,c);
+      this->set_min(p,sqr_d_right,i_right,c_right,sqr_d,i,c);
       looked_right = true;
     };
 

include/igl/copyleft/cgal/outer_hull.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "outer_hull.h"
 #include "extract_cells.h"
@@ -35,12 +35,12 @@ IGL_INLINE void igl::copyleft::cgal::outer_hull(
   // Exact types
   typedef CGAL::Epeck Kernel;
   typedef Kernel::FT ExactScalar;
-  typedef 
+  typedef
     Eigen::Matrix<
     ExactScalar,
     Eigen::Dynamic,
     Eigen::Dynamic,
-    DerivedHV::IsRowMajor> 
+    DerivedHV::IsRowMajor>
       MatrixXES;
   // Remesh self-intersections
   MatrixXES Vr;
@@ -540,4 +540,6 @@ template void igl::copyleft::cgal::outer_hull_legacy<Eigen::Matrix<CGAL::Lazy_ex
 template void igl::copyleft::cgal::outer_hull_legacy<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #include "../../barycenter.cpp"
 template void igl::barycenter<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> >&);
+#include "../../remove_unreferenced.cpp"
+template void igl::remove_unreferenced<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/flip_avoiding_line_search.cpp

@@ -6,252 +6,296 @@
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "flip_avoiding_line_search.h"
+#include "line_search.h"
 
 #include <Eigen/Dense>
 #include <vector>
-#include "line_search.h"
-#define TwoPi  2*M_PI//6.28318530717958648
-
-using namespace std;
-
-//---------------------------------------------------------------------------
-// x - array of size 3
-// In case 3 real roots: => x[0], x[1], x[2], return 3
-//         2 real roots: x[0], x[1],          return 2
-//         1 real root : x[0], x[1] ± i*x[2], return 1
-// http://math.ivanovo.ac.ru/dalgebra/Khashin/poly/index.html
-int SolveP3(std::vector<double>& x,double a,double b,double c) { // solve cubic equation x^3 + a*x^2 + b*x + c
-  double a2 = a*a;
-    double q  = (a2 - 3*b)/9;
-  double r  = (a*(2*a2-9*b) + 27*c)/54;
-    double r2 = r*r;
-  double q3 = q*q*q;
-  double A,B;
-    if(r2<q3) {
-        double t=r/sqrt(q3);
-    if( t<-1) t=-1;
-    if( t> 1) t= 1;
-        t=acos(t);
-        a/=3; q=-2*sqrt(q);
-        x[0]=q*cos(t/3)-a;
-        x[1]=q*cos((t+TwoPi)/3)-a;
-        x[2]=q*cos((t-TwoPi)/3)-a;
-        return(3);
-    } else {
-        A =-pow(fabs(r)+sqrt(r2-q3),1./3);
-    if( r<0 ) A=-A;
-    B = A==0? 0 : B=q/A;
-
-    a/=3;
-    x[0] =(A+B)-a;
-        x[1] =-0.5*(A+B)-a;
-        x[2] = 0.5*sqrt(3.)*(A-B);
-    if(fabs(x[2])<1e-14) { x[2]=x[1]; return(2); }
-        return(1);
-    }
-}
 
-double get_smallest_pos_quad_zero(double a,double b, double c) {
-  double t1,t2;
-  if (a != 0) {
-    double delta_in = pow(b,2) - 4*a*c;
-    if (delta_in < 0) {
-      return INFINITY;
+namespace igl
+{
+  namespace flip_avoiding
+  {
+    //---------------------------------------------------------------------------
+    // x - array of size 3
+    // In case 3 real roots: => x[0], x[1], x[2], return 3
+    //         2 real roots: x[0], x[1],          return 2
+    //         1 real root : x[0], x[1] ± i*x[2], return 1
+    // http://math.ivanovo.ac.ru/dalgebra/Khashin/poly/index.html
+    IGL_INLINE int SolveP3(std::vector<double>& x,double a,double b,double c)
+    { // solve cubic equation x^3 + a*x^2 + b*x + c
+      using namespace std;
+      double a2 = a*a;
+        double q  = (a2 - 3*b)/9;
+      double r  = (a*(2*a2-9*b) + 27*c)/54;
+        double r2 = r*r;
+      double q3 = q*q*q;
+      double A,B;
+        if(r2<q3)
+        {
+          double t=r/sqrt(q3);
+          if( t<-1) t=-1;
+          if( t> 1) t= 1;
+          t=acos(t);
+          a/=3; q=-2*sqrt(q);
+          x[0]=q*cos(t/3)-a;
+          x[1]=q*cos((t+(2*M_PI))/3)-a;
+          x[2]=q*cos((t-(2*M_PI))/3)-a;
+          return(3);
+        }
+        else
+        {
+          A =-pow(fabs(r)+sqrt(r2-q3),1./3);
+          if( r<0 ) A=-A;
+          B = A==0? 0 : B=q/A;
+
+          a/=3;
+          x[0] =(A+B)-a;
+          x[1] =-0.5*(A+B)-a;
+          x[2] = 0.5*sqrt(3.)*(A-B);
+          if(fabs(x[2])<1e-14)
+          {
+            x[2]=x[1]; return(2);
+          }
+          return(1);
+        }
     }
-    double delta = sqrt(delta_in);
-    t1 = (-b + delta)/ (2*a);
-    t2 = (-b - delta)/ (2*a);
-  } else {
-    t1 = t2 = -b/c;
-  }
-  assert (std::isfinite(t1));
-  assert (std::isfinite(t2));
 
-  double tmp_n = min(t1,t2);
-  t1 = max(t1,t2); t2 = tmp_n;
-  if (t1 == t2) {
-    return INFINITY; // means the orientation flips twice = doesn't flip
-  }
-  // return the smallest negative root if it exists, otherwise return infinity
-  if (t1 > 0) {
-    if (t2 > 0) {
-      return t2;
-    } else {
-      return t1;
+    IGL_INLINE double get_smallest_pos_quad_zero(double a,double b, double c)
+    {
+      using namespace std;
+      double t1,t2;
+      if (a != 0)
+      {
+        double delta_in = pow(b,2) - 4*a*c;
+        if (delta_in < 0)
+        {
+          return INFINITY;
+        }
+        double delta = sqrt(delta_in);
+        t1 = (-b + delta)/ (2*a);
+        t2 = (-b - delta)/ (2*a);
+      }
+      else
+      {
+        t1 = t2 = -b/c;
+      }
+      assert (std::isfinite(t1));
+      assert (std::isfinite(t2));
+
+      double tmp_n = min(t1,t2);
+      t1 = max(t1,t2); t2 = tmp_n;
+      if (t1 == t2)
+      {
+        return INFINITY; // means the orientation flips twice = doesn't flip
+      }
+      // return the smallest negative root if it exists, otherwise return infinity
+      if (t1 > 0)
+      {
+        if (t2 > 0)
+        {
+          return t2;
+        }
+        else
+        {
+          return t1;
+        }
+      }
+      else
+      {
+        return INFINITY;
+      }
     }
-  } else {
-    return INFINITY;
-  }
-}
-
- double get_min_pos_root_2D(const Eigen::MatrixXd& uv,const Eigen::MatrixXi& F,
-            Eigen::MatrixXd& d, int f) {
-/*
-      Finding the smallest timestep t s.t a triangle get degenerated (<=> det = 0)
-      The following code can be derived by a symbolic expression in matlab:
-
-      Symbolic matlab:
-      U11 = sym('U11');
-      U12 = sym('U12');
-      U21 = sym('U21');
-      U22 = sym('U22');
-      U31 = sym('U31');
-      U32 = sym('U32');
-
-      V11 = sym('V11');
-      V12 = sym('V12');
-      V21 = sym('V21');
-      V22 = sym('V22');
-      V31 = sym('V31');
-      V32 = sym('V32');
-
-      t = sym('t');
-
-      U1 = [U11,U12];
-      U2 = [U21,U22];
-      U3 = [U31,U32];
-
-      V1 = [V11,V12];
-      V2 = [V21,V22];
-      V3 = [V31,V32];
-
-      A = [(U2+V2*t) - (U1+ V1*t)];
-      B = [(U3+V3*t) - (U1+ V1*t)];
-      C = [A;B];
-
-      solve(det(C), t);
-      cf = coeffs(det(C),t); % Now cf(1),cf(2),cf(3) holds the coefficients for the polynom. at order c,b,a
-    */
-
-  int v1 = F(f,0); int v2 = F(f,1); int v3 = F(f,2);
-  // get quadratic coefficients (ax^2 + b^x + c)
-  #define U11 uv(v1,0)
-  #define U12 uv(v1,1)
-  #define U21 uv(v2,0)
-  #define U22 uv(v2,1)
-  #define U31 uv(v3,0)
-  #define U32 uv(v3,1)
-
-  #define V11 d(v1,0)
-  #define V12 d(v1,1)
-  #define V21 d(v2,0)
-  #define V22 d(v2,1)
-  #define V31 d(v3,0)
-  #define V32 d(v3,1)
-
-
-  double a = V11*V22 - V12*V21 - V11*V32 + V12*V31 + V21*V32 - V22*V31;
-  double b = U11*V22 - U12*V21 - U21*V12 + U22*V11 - U11*V32 + U12*V31 + U31*V12 - U32*V11 + U21*V32 - U22*V31 - U31*V22 + U32*V21;
-  double c = U11*U22 - U12*U21 - U11*U32 + U12*U31 + U21*U32 - U22*U31;
-
-  return get_smallest_pos_quad_zero(a,b,c);
-}
 
-double get_min_pos_root_3D(const Eigen::MatrixXd& uv,const Eigen::MatrixXi& F,
-            Eigen::MatrixXd& direc, int f) {
-  /*
-      Searching for the roots of:
-        +-1/6 * |ax ay az 1|
-                |bx by bz 1|
-                |cx cy cz 1|
-                |dx dy dz 1|
-      Every point ax,ay,az has a search direction a_dx,a_dy,a_dz, and so we add those to the matrix, and solve the cubic to find the step size t for a 0 volume
-      Symbolic matlab:
-        syms a_x a_y a_z a_dx a_dy a_dz % tetrahedera point and search direction
-        syms b_x b_y b_z b_dx b_dy b_dz
-        syms c_x c_y c_z c_dx c_dy c_dz
-        syms d_x d_y d_z d_dx d_dy d_dz
-        syms t % Timestep var, this is what we're looking for
-
-
-        a_plus_t = [a_x,a_y,a_z] + t*[a_dx,a_dy,a_dz];
-        b_plus_t = [b_x,b_y,b_z] + t*[b_dx,b_dy,b_dz];
-        c_plus_t = [c_x,c_y,c_z] + t*[c_dx,c_dy,c_dz];
-        d_plus_t = [d_x,d_y,d_z] + t*[d_dx,d_dy,d_dz];
-
-        vol_mat = [a_plus_t,1;b_plus_t,1;c_plus_t,1;d_plus_t,1]
-        //cf = coeffs(det(vol_det),t); % Now cf(1),cf(2),cf(3),cf(4) holds the coefficients for the polynom
-        [coefficients,terms] = coeffs(det(vol_det),t); % terms = [ t^3, t^2, t, 1], Coefficients hold the coeff we seek
-  */
-  int v1 = F(f,0); int v2 = F(f,1); int v3 = F(f,2); int v4 = F(f,3);
-  #define a_x uv(v1,0)
-  #define a_y uv(v1,1)
-  #define a_z uv(v1,2)
-  #define b_x uv(v2,0)
-  #define b_y uv(v2,1)
-  #define b_z uv(v2,2)
-  #define c_x uv(v3,0)
-  #define c_y uv(v3,1)
-  #define c_z uv(v3,2)
-  #define d_x uv(v4,0)
-  #define d_y uv(v4,1)
-  #define d_z uv(v4,2)
-
-  #define a_dx direc(v1,0)
-  #define a_dy direc(v1,1)
-  #define a_dz direc(v1,2)
-  #define b_dx direc(v2,0)
-  #define b_dy direc(v2,1)
-  #define b_dz direc(v2,2)
-  #define c_dx direc(v3,0)
-  #define c_dy direc(v3,1)
-  #define c_dz direc(v3,2)
-  #define d_dx direc(v4,0)
-  #define d_dy direc(v4,1)
-  #define d_dz direc(v4,2)
-
-  // Find solution for: a*t^3 + b*t^2 + c*d +d = 0
-  double a = a_dx*b_dy*c_dz - a_dx*b_dz*c_dy - a_dy*b_dx*c_dz + a_dy*b_dz*c_dx + a_dz*b_dx*c_dy - a_dz*b_dy*c_dx - a_dx*b_dy*d_dz + a_dx*b_dz*d_dy + a_dy*b_dx*d_dz - a_dy*b_dz*d_dx - a_dz*b_dx*d_dy + a_dz*b_dy*d_dx + a_dx*c_dy*d_dz - a_dx*c_dz*d_dy - a_dy*c_dx*d_dz + a_dy*c_dz*d_dx + a_dz*c_dx*d_dy - a_dz*c_dy*d_dx - b_dx*c_dy*d_dz + b_dx*c_dz*d_dy + b_dy*c_dx*d_dz - b_dy*c_dz*d_dx - b_dz*c_dx*d_dy + b_dz*c_dy*d_dx;
-  double b = a_dy*b_dz*c_x - a_dy*b_x*c_dz - a_dz*b_dy*c_x + a_dz*b_x*c_dy + a_x*b_dy*c_dz - a_x*b_dz*c_dy - a_dx*b_dz*c_y + a_dx*b_y*c_dz + a_dz*b_dx*c_y - a_dz*b_y*c_dx - a_y*b_dx*c_dz + a_y*b_dz*c_dx + a_dx*b_dy*c_z - a_dx*b_z*c_dy - a_dy*b_dx*c_z + a_dy*b_z*c_dx + a_z*b_dx*c_dy - a_z*b_dy*c_dx - a_dy*b_dz*d_x + a_dy*b_x*d_dz + a_dz*b_dy*d_x - a_dz*b_x*d_dy - a_x*b_dy*d_dz + a_x*b_dz*d_dy + a_dx*b_dz*d_y - a_dx*b_y*d_dz - a_dz*b_dx*d_y + a_dz*b_y*d_dx + a_y*b_dx*d_dz - a_y*b_dz*d_dx - a_dx*b_dy*d_z + a_dx*b_z*d_dy + a_dy*b_dx*d_z - a_dy*b_z*d_dx - a_z*b_dx*d_dy + a_z*b_dy*d_dx + a_dy*c_dz*d_x - a_dy*c_x*d_dz - a_dz*c_dy*d_x + a_dz*c_x*d_dy + a_x*c_dy*d_dz - a_x*c_dz*d_dy - a_dx*c_dz*d_y + a_dx*c_y*d_dz + a_dz*c_dx*d_y - a_dz*c_y*d_dx - a_y*c_dx*d_dz + a_y*c_dz*d_dx + a_dx*c_dy*d_z - a_dx*c_z*d_dy - a_dy*c_dx*d_z + a_dy*c_z*d_dx + a_z*c_dx*d_dy - a_z*c_dy*d_dx - b_dy*c_dz*d_x + b_dy*c_x*d_dz + b_dz*c_dy*d_x - b_dz*c_x*d_dy - b_x*c_dy*d_dz + b_x*c_dz*d_dy + b_dx*c_dz*d_y - b_dx*c_y*d_dz - b_dz*c_dx*d_y + b_dz*c_y*d_dx + b_y*c_dx*d_dz - b_y*c_dz*d_dx - b_dx*c_dy*d_z + b_dx*c_z*d_dy + b_dy*c_dx*d_z - b_dy*c_z*d_dx - b_z*c_dx*d_dy + b_z*c_dy*d_dx;
-  double c = a_dz*b_x*c_y - a_dz*b_y*c_x - a_x*b_dz*c_y + a_x*b_y*c_dz + a_y*b_dz*c_x - a_y*b_x*c_dz - a_dy*b_x*c_z + a_dy*b_z*c_x + a_x*b_dy*c_z - a_x*b_z*c_dy - a_z*b_dy*c_x + a_z*b_x*c_dy + a_dx*b_y*c_z - a_dx*b_z*c_y - a_y*b_dx*c_z + a_y*b_z*c_dx + a_z*b_dx*c_y - a_z*b_y*c_dx - a_dz*b_x*d_y + a_dz*b_y*d_x + a_x*b_dz*d_y - a_x*b_y*d_dz - a_y*b_dz*d_x + a_y*b_x*d_dz + a_dy*b_x*d_z - a_dy*b_z*d_x - a_x*b_dy*d_z + a_x*b_z*d_dy + a_z*b_dy*d_x - a_z*b_x*d_dy - a_dx*b_y*d_z + a_dx*b_z*d_y + a_y*b_dx*d_z - a_y*b_z*d_dx - a_z*b_dx*d_y + a_z*b_y*d_dx + a_dz*c_x*d_y - a_dz*c_y*d_x - a_x*c_dz*d_y + a_x*c_y*d_dz + a_y*c_dz*d_x - a_y*c_x*d_dz - a_dy*c_x*d_z + a_dy*c_z*d_x + a_x*c_dy*d_z - a_x*c_z*d_dy - a_z*c_dy*d_x + a_z*c_x*d_dy + a_dx*c_y*d_z - a_dx*c_z*d_y - a_y*c_dx*d_z + a_y*c_z*d_dx + a_z*c_dx*d_y - a_z*c_y*d_dx - b_dz*c_x*d_y + b_dz*c_y*d_x + b_x*c_dz*d_y - b_x*c_y*d_dz - b_y*c_dz*d_x + b_y*c_x*d_dz + b_dy*c_x*d_z - b_dy*c_z*d_x - b_x*c_dy*d_z + b_x*c_z*d_dy + b_z*c_dy*d_x - b_z*c_x*d_dy - b_dx*c_y*d_z + b_dx*c_z*d_y + b_y*c_dx*d_z - b_y*c_z*d_dx - b_z*c_dx*d_y + b_z*c_y*d_dx;
-  double d = a_x*b_y*c_z - a_x*b_z*c_y - a_y*b_x*c_z + a_y*b_z*c_x + a_z*b_x*c_y - a_z*b_y*c_x - a_x*b_y*d_z + a_x*b_z*d_y + a_y*b_x*d_z - a_y*b_z*d_x - a_z*b_x*d_y + a_z*b_y*d_x + a_x*c_y*d_z - a_x*c_z*d_y - a_y*c_x*d_z + a_y*c_z*d_x + a_z*c_x*d_y - a_z*c_y*d_x - b_x*c_y*d_z + b_x*c_z*d_y + b_y*c_x*d_z - b_y*c_z*d_x - b_z*c_x*d_y + b_z*c_y*d_x;
-
-  if (a==0) {
-    return get_smallest_pos_quad_zero(b,c,d);
-  }
-  b/=a; c/=a; d/=a; // normalize it all
-  std::vector<double> res(3);
-  int real_roots_num = SolveP3(res,b,c,d);
-  switch (real_roots_num) {
-    case 1:
-      return (res[0] >= 0) ? res[0]:INFINITY;
-    case 2: {
-      double max_root = max(res[0],res[1]); double min_root = min(res[0],res[1]);
-      if (min_root > 0) return min_root;
-      if (max_root > 0) return max_root;
-      return INFINITY;
+    IGL_INLINE double get_min_pos_root_2D(const Eigen::MatrixXd& uv,
+                                          const Eigen::MatrixXi& F,
+                                          Eigen::MatrixXd& d,
+                                          int f)
+    {
+      using namespace std;
+    /*
+          Finding the smallest timestep t s.t a triangle get degenerated (<=> det = 0)
+          The following code can be derived by a symbolic expression in matlab:
+
+          Symbolic matlab:
+          U11 = sym('U11');
+          U12 = sym('U12');
+          U21 = sym('U21');
+          U22 = sym('U22');
+          U31 = sym('U31');
+          U32 = sym('U32');
+
+          V11 = sym('V11');
+          V12 = sym('V12');
+          V21 = sym('V21');
+          V22 = sym('V22');
+          V31 = sym('V31');
+          V32 = sym('V32');
+
+          t = sym('t');
+
+          U1 = [U11,U12];
+          U2 = [U21,U22];
+          U3 = [U31,U32];
+
+          V1 = [V11,V12];
+          V2 = [V21,V22];
+          V3 = [V31,V32];
+
+          A = [(U2+V2*t) - (U1+ V1*t)];
+          B = [(U3+V3*t) - (U1+ V1*t)];
+          C = [A;B];
+
+          solve(det(C), t);
+          cf = coeffs(det(C),t); % Now cf(1),cf(2),cf(3) holds the coefficients for the polynom. at order c,b,a
+        */
+
+      int v1 = F(f,0); int v2 = F(f,1); int v3 = F(f,2);
+      // get quadratic coefficients (ax^2 + b^x + c)
+      const double& U11 = uv(v1,0);
+      const double& U12 = uv(v1,1);
+      const double& U21 = uv(v2,0);
+      const double& U22 = uv(v2,1);
+      const double& U31 = uv(v3,0);
+      const double& U32 = uv(v3,1);
+
+      const double& V11 = d(v1,0);
+      const double& V12 = d(v1,1);
+      const double& V21 = d(v2,0);
+      const double& V22 = d(v2,1);
+      const double& V31 = d(v3,0);
+      const double& V32 = d(v3,1);
+
+      double a = V11*V22 - V12*V21 - V11*V32 + V12*V31 + V21*V32 - V22*V31;
+      double b = U11*V22 - U12*V21 - U21*V12 + U22*V11 - U11*V32 + U12*V31 + U31*V12 - U32*V11 + U21*V32 - U22*V31 - U31*V22 + U32*V21;
+      double c = U11*U22 - U12*U21 - U11*U32 + U12*U31 + U21*U32 - U22*U31;
+
+      return get_smallest_pos_quad_zero(a,b,c);
     }
-    case 3:
-    default: {
-      std::sort(res.begin(),res.end());
-      if (res[0] > 0) return res[0];
-      if (res[1] > 0) return res[1];
-      if (res[2] > 0) return res[2];
-      return INFINITY;
-    }
-  }
 
-}
-
-double compute_max_step_from_singularities(const Eigen::MatrixXd& uv,
-                                            const Eigen::MatrixXi& F,
-                                            Eigen::MatrixXd& d) {
-    double max_step = INFINITY;
+    IGL_INLINE double get_min_pos_root_3D(const Eigen::MatrixXd& uv,
+                                          const Eigen::MatrixXi& F,
+                                          Eigen::MatrixXd& direc,
+                                          int f)
+    {
+      using namespace std;
+      /*
+          Searching for the roots of:
+            +-1/6 * |ax ay az 1|
+                    |bx by bz 1|
+                    |cx cy cz 1|
+                    |dx dy dz 1|
+          Every point ax,ay,az has a search direction a_dx,a_dy,a_dz, and so we add those to the matrix, and solve the cubic to find the step size t for a 0 volume
+          Symbolic matlab:
+            syms a_x a_y a_z a_dx a_dy a_dz % tetrahedera point and search direction
+            syms b_x b_y b_z b_dx b_dy b_dz
+            syms c_x c_y c_z c_dx c_dy c_dz
+            syms d_x d_y d_z d_dx d_dy d_dz
+            syms t % Timestep var, this is what we're looking for
+
+
+            a_plus_t = [a_x,a_y,a_z] + t*[a_dx,a_dy,a_dz];
+            b_plus_t = [b_x,b_y,b_z] + t*[b_dx,b_dy,b_dz];
+            c_plus_t = [c_x,c_y,c_z] + t*[c_dx,c_dy,c_dz];
+            d_plus_t = [d_x,d_y,d_z] + t*[d_dx,d_dy,d_dz];
+
+            vol_mat = [a_plus_t,1;b_plus_t,1;c_plus_t,1;d_plus_t,1]
+            //cf = coeffs(det(vol_det),t); % Now cf(1),cf(2),cf(3),cf(4) holds the coefficients for the polynom
+            [coefficients,terms] = coeffs(det(vol_det),t); % terms = [ t^3, t^2, t, 1], Coefficients hold the coeff we seek
+      */
+      int v1 = F(f,0); int v2 = F(f,1); int v3 = F(f,2); int v4 = F(f,3);
+      const double& a_x = uv(v1,0);
+      const double& a_y = uv(v1,1);
+      const double& a_z = uv(v1,2);
+      const double& b_x = uv(v2,0);
+      const double& b_y = uv(v2,1);
+      const double& b_z = uv(v2,2);
+      const double& c_x = uv(v3,0);
+      const double& c_y = uv(v3,1);
+      const double& c_z = uv(v3,2);
+      const double& d_x = uv(v4,0);
+      const double& d_y = uv(v4,1);
+      const double& d_z = uv(v4,2);
+
+      const double& a_dx = direc(v1,0);
+      const double& a_dy = direc(v1,1);
+      const double& a_dz = direc(v1,2);
+      const double& b_dx = direc(v2,0);
+      const double& b_dy = direc(v2,1);
+      const double& b_dz = direc(v2,2);
+      const double& c_dx = direc(v3,0);
+      const double& c_dy = direc(v3,1);
+      const double& c_dz = direc(v3,2);
+      const double& d_dx = direc(v4,0);
+      const double& d_dy = direc(v4,1);
+      const double& d_dz = direc(v4,2);
+
+      // Find solution for: a*t^3 + b*t^2 + c*d +d = 0
+      double a = a_dx*b_dy*c_dz - a_dx*b_dz*c_dy - a_dy*b_dx*c_dz + a_dy*b_dz*c_dx + a_dz*b_dx*c_dy - a_dz*b_dy*c_dx - a_dx*b_dy*d_dz + a_dx*b_dz*d_dy + a_dy*b_dx*d_dz - a_dy*b_dz*d_dx - a_dz*b_dx*d_dy + a_dz*b_dy*d_dx + a_dx*c_dy*d_dz - a_dx*c_dz*d_dy - a_dy*c_dx*d_dz + a_dy*c_dz*d_dx + a_dz*c_dx*d_dy - a_dz*c_dy*d_dx - b_dx*c_dy*d_dz + b_dx*c_dz*d_dy + b_dy*c_dx*d_dz - b_dy*c_dz*d_dx - b_dz*c_dx*d_dy + b_dz*c_dy*d_dx;
+
+      double b = a_dy*b_dz*c_x - a_dy*b_x*c_dz - a_dz*b_dy*c_x + a_dz*b_x*c_dy + a_x*b_dy*c_dz - a_x*b_dz*c_dy - a_dx*b_dz*c_y + a_dx*b_y*c_dz + a_dz*b_dx*c_y - a_dz*b_y*c_dx - a_y*b_dx*c_dz + a_y*b_dz*c_dx + a_dx*b_dy*c_z - a_dx*b_z*c_dy - a_dy*b_dx*c_z + a_dy*b_z*c_dx + a_z*b_dx*c_dy - a_z*b_dy*c_dx - a_dy*b_dz*d_x + a_dy*b_x*d_dz + a_dz*b_dy*d_x - a_dz*b_x*d_dy - a_x*b_dy*d_dz + a_x*b_dz*d_dy + a_dx*b_dz*d_y - a_dx*b_y*d_dz - a_dz*b_dx*d_y + a_dz*b_y*d_dx + a_y*b_dx*d_dz - a_y*b_dz*d_dx - a_dx*b_dy*d_z + a_dx*b_z*d_dy + a_dy*b_dx*d_z - a_dy*b_z*d_dx - a_z*b_dx*d_dy + a_z*b_dy*d_dx + a_dy*c_dz*d_x - a_dy*c_x*d_dz - a_dz*c_dy*d_x + a_dz*c_x*d_dy + a_x*c_dy*d_dz - a_x*c_dz*d_dy - a_dx*c_dz*d_y + a_dx*c_y*d_dz + a_dz*c_dx*d_y - a_dz*c_y*d_dx - a_y*c_dx*d_dz + a_y*c_dz*d_dx + a_dx*c_dy*d_z - a_dx*c_z*d_dy - a_dy*c_dx*d_z + a_dy*c_z*d_dx + a_z*c_dx*d_dy - a_z*c_dy*d_dx - b_dy*c_dz*d_x + b_dy*c_x*d_dz + b_dz*c_dy*d_x - b_dz*c_x*d_dy - b_x*c_dy*d_dz + b_x*c_dz*d_dy + b_dx*c_dz*d_y - b_dx*c_y*d_dz - b_dz*c_dx*d_y + b_dz*c_y*d_dx + b_y*c_dx*d_dz - b_y*c_dz*d_dx - b_dx*c_dy*d_z + b_dx*c_z*d_dy + b_dy*c_dx*d_z - b_dy*c_z*d_dx - b_z*c_dx*d_dy + b_z*c_dy*d_dx;
+
+      double c = a_dz*b_x*c_y - a_dz*b_y*c_x - a_x*b_dz*c_y + a_x*b_y*c_dz + a_y*b_dz*c_x - a_y*b_x*c_dz - a_dy*b_x*c_z + a_dy*b_z*c_x + a_x*b_dy*c_z - a_x*b_z*c_dy - a_z*b_dy*c_x + a_z*b_x*c_dy + a_dx*b_y*c_z - a_dx*b_z*c_y - a_y*b_dx*c_z + a_y*b_z*c_dx + a_z*b_dx*c_y - a_z*b_y*c_dx - a_dz*b_x*d_y + a_dz*b_y*d_x + a_x*b_dz*d_y - a_x*b_y*d_dz - a_y*b_dz*d_x + a_y*b_x*d_dz + a_dy*b_x*d_z - a_dy*b_z*d_x - a_x*b_dy*d_z + a_x*b_z*d_dy + a_z*b_dy*d_x - a_z*b_x*d_dy - a_dx*b_y*d_z + a_dx*b_z*d_y + a_y*b_dx*d_z - a_y*b_z*d_dx - a_z*b_dx*d_y + a_z*b_y*d_dx + a_dz*c_x*d_y - a_dz*c_y*d_x - a_x*c_dz*d_y + a_x*c_y*d_dz + a_y*c_dz*d_x - a_y*c_x*d_dz - a_dy*c_x*d_z + a_dy*c_z*d_x + a_x*c_dy*d_z - a_x*c_z*d_dy - a_z*c_dy*d_x + a_z*c_x*d_dy + a_dx*c_y*d_z - a_dx*c_z*d_y - a_y*c_dx*d_z + a_y*c_z*d_dx + a_z*c_dx*d_y - a_z*c_y*d_dx - b_dz*c_x*d_y + b_dz*c_y*d_x + b_x*c_dz*d_y - b_x*c_y*d_dz - b_y*c_dz*d_x + b_y*c_x*d_dz + b_dy*c_x*d_z - b_dy*c_z*d_x - b_x*c_dy*d_z + b_x*c_z*d_dy + b_z*c_dy*d_x - b_z*c_x*d_dy - b_dx*c_y*d_z + b_dx*c_z*d_y + b_y*c_dx*d_z - b_y*c_z*d_dx - b_z*c_dx*d_y + b_z*c_y*d_dx;
+
+      double d = a_x*b_y*c_z - a_x*b_z*c_y - a_y*b_x*c_z + a_y*b_z*c_x + a_z*b_x*c_y - a_z*b_y*c_x - a_x*b_y*d_z + a_x*b_z*d_y + a_y*b_x*d_z - a_y*b_z*d_x - a_z*b_x*d_y + a_z*b_y*d_x + a_x*c_y*d_z - a_x*c_z*d_y - a_y*c_x*d_z + a_y*c_z*d_x + a_z*c_x*d_y - a_z*c_y*d_x - b_x*c_y*d_z + b_x*c_z*d_y + b_y*c_x*d_z - b_y*c_z*d_x - b_z*c_x*d_y + b_z*c_y*d_x;
+
+      if (a==0)
+      {
+        return get_smallest_pos_quad_zero(b,c,d);
+      }
+      b/=a; c/=a; d/=a; // normalize it all
+      std::vector<double> res(3);
+      int real_roots_num = SolveP3(res,b,c,d);
+      switch (real_roots_num)
+      {
+        case 1:
+          return (res[0] >= 0) ? res[0]:INFINITY;
+        case 2:
+        {
+          double max_root = max(res[0],res[1]); double min_root = min(res[0],res[1]);
+          if (min_root > 0) return min_root;
+          if (max_root > 0) return max_root;
+          return INFINITY;
+        }
+        case 3:
+        default:
+        {
+          std::sort(res.begin(),res.end());
+          if (res[0] > 0) return res[0];
+          if (res[1] > 0) return res[1];
+          if (res[2] > 0) return res[2];
+          return INFINITY;
+        }
+      }
+    }
 
-    // The if statement is outside the for loops to avoid branching/ease parallelizing
-    if (uv.cols() == 2) {
-      for (int f = 0; f < F.rows(); f++) {
-        double min_positive_root = get_min_pos_root_2D(uv,F,d,f);
-        max_step = min(max_step, min_positive_root);
+    IGL_INLINE double compute_max_step_from_singularities(const Eigen::MatrixXd& uv,
+                                                          const Eigen::MatrixXi& F,
+                                                          Eigen::MatrixXd& d)
+    {
+      using namespace std;
+      double max_step = INFINITY;
+
+      // The if statement is outside the for loops to avoid branching/ease parallelizing
+      if (uv.cols() == 2)
+      {
+        for (int f = 0; f < F.rows(); f++)
+        {
+          double min_positive_root = get_min_pos_root_2D(uv,F,d,f);
+          max_step = min(max_step, min_positive_root);
+        }
       }
-    } else { // volumetric deformation
-      for (int f = 0; f < F.rows(); f++) {
-        double min_positive_root = get_min_pos_root_3D(uv,F,d,f);
-        max_step = min(max_step, min_positive_root);
+      else
+      { // volumetric deformation
+        for (int f = 0; f < F.rows(); f++)
+        {
+          double min_positive_root = get_min_pos_root_3D(uv,F,d,f);
+          max_step = min(max_step, min_positive_root);
+        }
       }
+      return max_step;
     }
-    return max_step;
- }
+  }
+}
 
 IGL_INLINE double igl::flip_avoiding_line_search(
   const Eigen::MatrixXi F,
@@ -260,12 +304,13 @@ IGL_INLINE double igl::flip_avoiding_line_search(
   std::function<double(Eigen::MatrixXd&)> energy,
   double cur_energy)
 {
-    Eigen::MatrixXd d = dst_v - cur_v;
+  using namespace std;
+  Eigen::MatrixXd d = dst_v - cur_v;
 
-    double min_step_to_singularity = compute_max_step_from_singularities(cur_v,F,d);
-    double max_step_size = min(1., min_step_to_singularity*0.8);
+  double min_step_to_singularity = igl::flip_avoiding::compute_max_step_from_singularities(cur_v,F,d);
+  double max_step_size = min(1., min_step_to_singularity*0.8);
 
-    return igl::line_search(cur_v,d,max_step_size, energy, cur_energy);
+  return igl::line_search(cur_v,d,max_step_size, energy, cur_energy);
 }
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/flip_avoiding_line_search.h

@@ -7,13 +7,13 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #ifndef IGL_FLIP_AVOIDING_LINE_SEARCH_H
 #define IGL_FLIP_AVOIDING_LINE_SEARCH_H
-#include <igl/igl_inline.h>
+#include "igl_inline.h"
 
 #include <Eigen/Dense>
 
 namespace igl
 {
-  // A bisection line search for a mesh based energy that avoids triangle flips as suggested in 
+  // A bisection line search for a mesh based energy that avoids triangle flips as suggested in
   // 		"Bijective Parameterization with Free Boundaries" (Smith J. and Schaefer S., 2015).
   //
   // The user specifies an initial vertices position (that has no flips) and target one (that my have flipped triangles).
@@ -27,7 +27,7 @@ namespace igl
   //   cur_v  						#V by dim list of variables
   //   dst_v  						#V by dim list of target vertices. This mesh may have flipped triangles
   //   energy       			    A function to compute the mesh-based energy (return an energy that is bigger than 0)
-  //   cur_energy(OPTIONAL)         The energy at the given point. Helps save redundant computations. 
+  //   cur_energy(OPTIONAL)         The energy at the given point. Helps save redundant computations.
   //							    This is optional. If not specified, the function will compute it.
   // Outputs:
   //		cur_v  						#V by dim list of variables at the new location

include/igl/line_search.cpp

@@ -15,21 +15,28 @@ IGL_INLINE double igl::line_search(
   double cur_energy)
 {
   double old_energy;
-  if (cur_energy > 0) {
+  if (cur_energy > 0)
+  {
     old_energy = cur_energy;
-  } else {
+  }
+  else
+  {
     old_energy = energy(x); // no energy was given -> need to compute the current energy
   }
   double new_energy = old_energy;
   int cur_iter = 0; int MAX_STEP_SIZE_ITER = 12;
 
-  while (new_energy >= old_energy && cur_iter < MAX_STEP_SIZE_ITER) {
+  while (new_energy >= old_energy && cur_iter < MAX_STEP_SIZE_ITER)
+  {
     Eigen::MatrixXd new_x = x + step_size * d;
 
     double cur_e = energy(new_x);
-    if ( cur_e >= old_energy) {
+    if ( cur_e >= old_energy)
+    {
       step_size /= 2;
-    } else {
+    }
+    else
+    {
       x = new_x;
       new_energy = cur_e;
     }

include/igl/line_search.h

@@ -7,7 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #ifndef IGL_LINE_SEARCH_H
 #define IGL_LINE_SEARCH_H
-#include <igl/igl_inline.h>
+#include "igl_inline.h"
 
 #include <Eigen/Dense>
 
@@ -21,7 +21,7 @@ namespace igl
   //   d  						#X by dim list of a given search direction
   //   i_step_size  			initial step size
   //   energy       			A function to compute the mesh-based energy (return an energy that is bigger than 0)
-  //   cur_energy(OPTIONAL)     The energy at the given point. Helps save redundant computations. 
+  //   cur_energy(OPTIONAL)     The energy at the given point. Helps save redundant computations.
   //							This is optional. If not specified, the function will compute it.
   // Outputs:
   //		x  						#X by dim list of variables at the new location

include/igl/slim.h

@@ -8,28 +8,21 @@
 #ifndef SLIM_H
 #define SLIM_H
 
+#include "igl_inline.h"
 #include <Eigen/Dense>
 #include <Eigen/Sparse>
 
-#include <string>
-
-#include <igl/jet.h>
-#include <igl/readOBJ.h>
-#include <igl/facet_components.h>
-#include <igl/slice.h>
-
-class WeightedGlobalLocal;
-
 namespace igl
 {
 
 // Compute a SLIM map as derived in "Scalable Locally Injective Maps" [Rabinovich et al. 2016].
-struct SLIMData {
-
+struct SLIMData
+{
   // Input
   Eigen::MatrixXd V; // #V by 3 list of mesh vertex positions
   Eigen::MatrixXi F; // #F by 3/3 list of mesh faces (triangles/tets)
-  enum SLIM_ENERGY {
+  enum SLIM_ENERGY
+  {
     ARAP,
     LOG_ARAP,
     SYMMETRIC_DIRICHLET,
@@ -40,10 +33,10 @@ struct SLIMData {
   SLIM_ENERGY slim_energy;
 
   // Optional Input
-    // soft constraints
-    Eigen::VectorXi b;
-    Eigen::MatrixXd bc;
-    double soft_const_p;
+  // soft constraints
+  Eigen::VectorXi b;
+  Eigen::MatrixXd bc;
+  double soft_const_p;
 
   double exp_factor; // used for exponential energies, ignored otherwise
   bool mesh_improvement_3d; // only supported for 3d
@@ -60,30 +53,45 @@ struct SLIMData {
   int f_num;
   double proximal_p;
 
-  WeightedGlobalLocal* wGlobalLocal;
+  Eigen::VectorXd WGL_M;
+  Eigen::VectorXd rhs;
+  Eigen::MatrixXd Ri,Ji;
+  Eigen::VectorXd W_11; Eigen::VectorXd W_12; Eigen::VectorXd W_13;
+  Eigen::VectorXd W_21; Eigen::VectorXd W_22; Eigen::VectorXd W_23;
+  Eigen::VectorXd W_31; Eigen::VectorXd W_32; Eigen::VectorXd W_33;
+  Eigen::SparseMatrix<double> Dx,Dy,Dz;
+  int f_n,v_n;
+  bool first_solve;
+  bool has_pre_calc = false;
+  int dim;
 };
 
-  // Compute necessary information to start using SLIM
-  // Inputs:
-  //		V           #V by 3 list of mesh vertex positions
-  //		F           #F by 3/3 list of mesh faces (triangles/tets)
-  //    b           list of boundary indices into V
-  //    bc          #b by dim list of boundary conditions
-  //    soft_p      Soft penalty factor (can be zero)
-  //    slim_energy Energy to minimize
-IGL_INLINE void slim_precompute(Eigen::MatrixXd& V, Eigen::MatrixXi& F, Eigen::MatrixXd& V_init, SLIMData& data,
-   SLIMData::SLIM_ENERGY slim_energy, Eigen::VectorXi& b, Eigen::MatrixXd& bc, double soft_p);
+// Compute necessary information to start using SLIM
+// Inputs:
+//		V           #V by 3 list of mesh vertex positions
+//		F           #F by 3/3 list of mesh faces (triangles/tets)
+//    b           list of boundary indices into V
+//    bc          #b by dim list of boundary conditions
+//    soft_p      Soft penalty factor (can be zero)
+//    slim_energy Energy to minimize
+IGL_INLINE void slim_precompute(Eigen::MatrixXd& V,
+                                Eigen::MatrixXi& F,
+                                Eigen::MatrixXd& V_init,
+                                SLIMData& data,
+                                SLIMData::SLIM_ENERGY slim_energy,
+                                Eigen::VectorXi& b,
+                                Eigen::MatrixXd& bc,
+                                double soft_p);
 
 // Run iter_num iterations of SLIM
 // Outputs:
 //    V_o (in SLIMData): #V by dim list of mesh vertex positions
-IGL_INLINE void slim_solve(SLIMData& data, int iter_num);
+IGL_INLINE Eigen::MatrixXd slim_solve(SLIMData& data, int iter_num);
 
-}
+} // END NAMESPACE
 
 #ifndef IGL_STATIC_LIBRARY
 #  include "slim.cpp"
 #endif
 
-
 #endif // SLIM_H

include/igl/upsample.cpp

@@ -11,31 +11,29 @@
 
 
 template <
-  typename DerivedV,
   typename DerivedF,
-  typename DerivedNV,
+  typename DerivedS,
   typename DerivedNF>
-IGL_INLINE void igl::upsample(
-  const Eigen::PlainObjectBase<DerivedV>& V,
-  const Eigen::PlainObjectBase<DerivedF>& F,
-  Eigen::PlainObjectBase<DerivedNV>& NV,
-  Eigen::PlainObjectBase<DerivedNF>& NF)
+IGL_INLINE void igl::upsample(const int n_verts,
+                              const Eigen::PlainObjectBase<DerivedF>& F,
+                              Eigen::SparseMatrix<DerivedS>& S,
+                              Eigen::PlainObjectBase<DerivedNF>& NF)
 {
-  // Use "in place" wrapper instead
-  assert(&V != &NV);
-  assert(&F != &NF);
   using namespace std;
   using namespace Eigen;
 
+  typedef Eigen::Triplet<DerivedS> Triplet_t;
+
   Eigen::Matrix<
-    typename DerivedF::Scalar,Eigen::Dynamic,Eigen::Dynamic>
-    FF,FFi;
+  typename DerivedF::Scalar,Eigen::Dynamic,Eigen::Dynamic>
+  FF,FFi;
   triangle_triangle_adjacency(F,FF,FFi);
 
   // TODO: Cache optimization missing from here, it is a mess
 
   // Compute the number and positions of the vertices to insert (on edges)
   Eigen::MatrixXi NI = Eigen::MatrixXi::Constant(FF.rows(),FF.cols(),-1);
+  Eigen::MatrixXi NIdoubles = Eigen::MatrixXi::Zero(FF.rows(), FF.cols());
   int counter = 0;
 
   for(int i=0;i<FF.rows();++i)
@@ -45,33 +43,45 @@ IGL_INLINE void igl::upsample(
       if(NI(i,j) == -1)
       {
         NI(i,j) = counter;
-        if (FF(i,j) != -1) // If it is not a border
-          NI(FF(i,j),FFi(i,j)) = counter;
+        NIdoubles(i,j) = 0;
+        if (FF(i,j) != -1) {
+          //If it is not a boundary
+          NI(FF(i,j), FFi(i,j)) = counter;
+          NIdoubles(i,j) = 1;
+        }
         ++counter;
       }
     }
   }
 
-  int n_odd = V.rows();
-  int n_even = counter;
+  const int& n_odd = n_verts;
+  const int& n_even = counter;
+  const int n_newverts = n_odd + n_even;
 
-  // Preallocate NV and NF
-  NV.resize(V.rows()+n_even,V.cols());
-  NF.resize(F.rows()*4,3);
+  //Construct vertex positions
+  std::vector<Triplet_t> tripletList;
 
   // Fill the odd vertices position
-  NV.block(0,0,V.rows(),V.cols()) = V;
+  for (int i=0; i<n_odd; ++i)
+  {
+    tripletList.emplace_back(i, i, 1.);
+  }
 
-  // Fill the even vertices position
   for(int i=0;i<FF.rows();++i)
   {
     for(int j=0;j<3;++j)
     {
-      NV.row(NI(i,j) + n_odd) = 0.5 * V.row(F(i,j)) + 0.5 * V.row(F(i,(j+1)%3));
+      if(NIdoubles(i,j)==0) {
+        tripletList.emplace_back(NI(i,j) + n_odd, F(i,j), 1./2.);
+        tripletList.emplace_back(NI(i,j) + n_odd, F(i,(j+1)%3), 1./2.);
+      }
     }
   }
+  S.resize(n_newverts, n_verts);
+  S.setFromTriplets(tripletList.begin(), tripletList.end());
 
   // Build the new topology (Every face is replaced by four)
+  NF.resize(F.rows()*4,3);
   for(int i=0; i<F.rows();++i)
   {
     VectorXi VI(6);
@@ -88,19 +98,43 @@ IGL_INLINE void igl::upsample(
     NF.row((i*4)+2) = f2;
     NF.row((i*4)+3) = f3;
   }
+}
 
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedNV,
+  typename DerivedNF>
+IGL_INLINE void igl::upsample(
+                              const Eigen::PlainObjectBase<DerivedV>& V,
+                              const Eigen::PlainObjectBase<DerivedF>& F,
+                              Eigen::PlainObjectBase<DerivedNV>& NV,
+                              Eigen::PlainObjectBase<DerivedNF>& NF,
+                              const int number_of_subdivs)
+{
+  typedef Eigen::SparseMatrix<double> SparseMat;
+  typedef Eigen::Triplet<double> Triplet_t;
+
+  NV = V;
+  NF = F;
+  for(int i=0; i<number_of_subdivs; ++i) {
+    DerivedNF tempF = NF;
+    SparseMat S;
+    upsample(NV.rows(), tempF, S, NF);
+    NV = S*NV;
+  }
 }
 
 template <
   typename MatV,
   typename MatF>
-IGL_INLINE void igl::upsample(
-  MatV& V,
-  MatF& F)
+IGL_INLINE void igl::upsample(MatV& V,
+                              MatF& F,
+                              const int number_of_subdivs)
 {
   const MatV V_copy = V;
   const MatF F_copy = F;
-  return upsample(V_copy,F_copy,V,F);
+  return upsample(V_copy,F_copy,V,F,number_of_subdivs);
 }
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/upsample.h

@@ -10,11 +10,31 @@
 #include "igl_inline.h"
 
 #include <Eigen/Core>
+#include <Eigen/Sparse>
 
 // History:
 //  changed templates from generic matrices to PlainObjectBase Alec May 7, 2011
 namespace igl
 {
+  // Subdivide without moving vertices: Given the triangle mesh [V, F],
+  // where n_verts = V.rows(), computes newV and a sparse matrix S s.t.
+  // [newV, newF] is the subdivided mesh where newV = S*V.
+  //
+  // Inputs:
+  //  n_verts an integer (number of mesh vertices)
+  //  F an m by 3 matrix of integers of triangle faces
+  // Outputs:
+  //  S a sparse matrix (will become the subdivision matrix)
+  //  newF a matrix containing the new faces
+  template <
+    typename DerivedF,
+    typename DerivedS,
+    typename DerivedNF>
+  IGL_INLINE void upsample(const int n_verts,
+                           const Eigen::PlainObjectBase<DerivedF>& F,
+                           Eigen::SparseMatrix<DerivedS>& S,
+                           Eigen::PlainObjectBase<DerivedNF>& NF);
+
   // Subdivide a mesh without moving vertices: loop subdivision but odd
   // vertices stay put and even vertices are just edge midpoints
   // 
@@ -38,18 +58,19 @@ namespace igl
     typename DerivedF,
     typename DerivedNV,
     typename DerivedNF>
-  IGL_INLINE void upsample(
-    const Eigen::PlainObjectBase<DerivedV>& V,
-    const Eigen::PlainObjectBase<DerivedF>& F,
-    Eigen::PlainObjectBase<DerivedNV>& NV,
-    Eigen::PlainObjectBase<DerivedNF>& NF);
+  IGL_INLINE void upsample(const Eigen::PlainObjectBase<DerivedV>& V,
+                           const Eigen::PlainObjectBase<DerivedF>& F,
+                           Eigen::PlainObjectBase<DerivedNV>& NV,
+                           Eigen::PlainObjectBase<DerivedNF>& NF,
+                           const int number_of_subdivs = 1);
+
   // Virtually in place wrapper
   template <
     typename MatV, 
     typename MatF>
-  IGL_INLINE void upsample(
-    MatV& V,
-    MatF& F);
+  IGL_INLINE void upsample(MatV& V,
+                           MatF& F,
+                           const int number_of_subdivs = 1);
 }
 
 #ifndef IGL_STATIC_LIBRARY