* integrate code for scalable locally injective maps

* fixes for static build
* added tutorial 611


Former-commit-id: 6b7a58cabd3cadb4beb3dbb0c65031e424e14ff8

include/igl/components.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2013 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 "components.h"
 #include "adjacency_matrix.h"
@@ -90,4 +90,5 @@ IGL_INLINE void igl::components(
 template void igl::components<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::components<int, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<int, 0, int> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::components<int, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::SparseMatrix<int, 0, int> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::components<double, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/euler_characteristic.cpp

@@ -29,4 +29,5 @@ IGL_INLINE int igl::euler_characteristic(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh
+template int igl::euler_characteristic<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
 #endif

include/igl/flip_avoiding_line_search.cpp

@@ -0,0 +1,272 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Michael Rabinovich
+//
+// 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 "flip_avoiding_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;
+    }
+    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;
+  }
+}
+
+ 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;
+    }
+    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;
+
+    // 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);
+      }
+    }
+    return max_step;
+ }
+
+IGL_INLINE double igl::flip_avoiding_line_search(
+  const Eigen::MatrixXi F,
+  Eigen::MatrixXd& cur_v,
+  Eigen::MatrixXd& dst_v,
+  std::function<double(Eigen::MatrixXd&)> energy,
+  double cur_energy)
+{
+    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);
+
+    return igl::line_search(cur_v,d,max_step_size, energy, cur_energy);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+#endif

include/igl/flip_avoiding_line_search.h

@@ -0,0 +1,30 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Michael Rabinovich
+//
+// 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/.
+#ifndef IGL_FLIP_AVOIDING_LINE_SEARCH_H
+#define IGL_FLIP_AVOIDING_LINE_SEARCH_H
+#include <igl/igl_inline.h>
+
+#include <Eigen/Dense>
+
+namespace igl
+{
+  // TODO DOCUMENTATION MISSING
+  IGL_INLINE double flip_avoiding_line_search(
+    const Eigen::MatrixXi F,
+    Eigen::MatrixXd& cur_v,
+    Eigen::MatrixXd& dst_v,
+    std::function<double(Eigen::MatrixXd&)> energy,
+    double cur_energy);
+
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "flip_avoiding_line_search.cpp"
+#endif
+
+#endif

include/igl/flipped_triangles.cpp

@@ -17,12 +17,12 @@ IGL_INLINE void igl::flipped_triangles(
 {
   assert(V.cols() == 2 && "V should contain 2D positions");
   std::vector<typename DerivedX::Scalar> flip_idx;
-  for (int i = 0; i < F.rows(); i++) 
+  for (int i = 0; i < F.rows(); i++)
   {
     // https://www.cs.cmu.edu/~quake/robust.html
     typedef Eigen::Matrix<typename DerivedV::Scalar,1,2> RowVector2S;
-    RowVector2S v1_n = V.row(F(i,0)); 
-    RowVector2S v2_n = V.row(F(i,1)); 
+    RowVector2S v1_n = V.row(F(i,0));
+    RowVector2S v2_n = V.row(F(i,1));
     RowVector2S v3_n = V.row(F(i,2));
     Eigen::Matrix<typename DerivedV::Scalar,3,3> T2_Homo;
     T2_Homo.col(0) << v1_n(0),v1_n(1),1.;
@@ -30,7 +30,7 @@ IGL_INLINE void igl::flipped_triangles(
     T2_Homo.col(2) << v3_n(0),v3_n(1),1.;
     double det = T2_Homo.determinant();
     assert(det == det && "det should not be NaN");
-    if (det < 0) 
+    if (det < 0)
     {
       flip_idx.push_back(i);
     }
@@ -51,4 +51,5 @@ IGL_INLINE Eigen::VectorXi igl::flipped_triangles(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template void igl::flipped_triangles<Eigen::Matrix<double, -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<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template Eigen::Matrix<int, -1, 1, 0, -1, 1> igl::flipped_triangles<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
 #endif

include/igl/line_search.cpp

@@ -0,0 +1,43 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Michael Rabinovich
+//
+// 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 "line_search.h"
+
+IGL_INLINE double igl::line_search(
+  Eigen::MatrixXd& x,
+  const Eigen::MatrixXd& d,
+  double step_size,
+  std::function<double(Eigen::MatrixXd&)> energy,
+  double cur_energy)
+{
+  double old_energy;
+  if (cur_energy > 0) {
+    old_energy = cur_energy;
+  } 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) {
+    Eigen::MatrixXd new_x = x + step_size * d;
+
+    double cur_e = energy(new_x);
+    if ( cur_e >= old_energy) {
+      step_size /= 2;
+    } else {
+      x = new_x;
+      new_energy = cur_e;
+    }
+    cur_iter++;
+  }
+  return new_energy;
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+#endif

include/igl/line_search.h

@@ -0,0 +1,30 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Michael Rabinovich
+//
+// 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/.
+#ifndef IGL_LINE_SEARCH_H
+#define IGL_LINE_SEARCH_H
+#include <igl/igl_inline.h>
+
+#include <Eigen/Dense>
+
+namespace igl
+{
+  // TODO DOCUMENTATION MISSING
+  IGL_INLINE double line_search(
+    Eigen::MatrixXd& x,
+    const Eigen::MatrixXd& d,
+    double step_size,
+    std::function<double(Eigen::MatrixXd&)> energy,
+    double cur_energy);
+
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "line_search.cpp"
+#endif
+
+#endif

include/igl/slim.cpp

@@ -0,0 +1,852 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Michael Rabinovich
+//
+// 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 "slim.h"
+
+#include <igl/boundary_loop.h>
+#include <igl/cotmatrix.h>
+#include <igl/edge_lengths.h>
+#include <igl/grad.h>
+#include <igl/local_basis.h>
+#include <igl/readOBJ.h>
+#include <igl/repdiag.h>
+#include <igl/vector_area_matrix.h>
+#include <iostream>
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+#include <map>
+#include <set>
+#include <vector>
+
+#include "igl/arap.h"
+#include "igl/cat.h"
+#include "igl/doublearea.h"
+#include "igl/grad.h"
+#include "igl/local_basis.h"
+#include "igl/per_face_normals.h"
+#include "igl/slice_into.h"
+#include "igl/volume.h"
+#include "igl/polar_svd.h"
+
+#include <Eigen/IterativeLinearSolvers>
+#include <Eigen/Sparse>
+#include <Eigen/SparseQR>
+#include <Eigen/SparseCholesky>
+#include <Eigen/IterativeLinearSolvers>
+
+#include "flip_avoiding_line_search.h"
+
+using namespace std;
+using namespace Eigen;
+
+///////// Helper functions to compute gradient matrices
+
+void compute_surface_gradient_matrix(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F,
+                                     const Eigen::MatrixXd& F1, const Eigen::MatrixXd& F2,
+                                     Eigen::SparseMatrix<double>& D1, Eigen::SparseMatrix<double>& D2) {
+
+  Eigen::SparseMatrix<double> G;
+  igl::grad(V,F,G);
+  Eigen::SparseMatrix<double> Dx = G.block(0,0,F.rows(),V.rows());
+  Eigen::SparseMatrix<double> Dy = G.block(F.rows(),0,F.rows(),V.rows());
+  Eigen::SparseMatrix<double> Dz = G.block(2*F.rows(),0,F.rows(),V.rows());
+
+  D1 = F1.col(0).asDiagonal()*Dx + F1.col(1).asDiagonal()*Dy + F1.col(2).asDiagonal()*Dz;
+  D2 = F2.col(0).asDiagonal()*Dx + F2.col(1).asDiagonal()*Dy + F2.col(2).asDiagonal()*Dz;
+}
+
+void compute_tet_grad_matrix(const Eigen::MatrixXd& V, const Eigen::MatrixXi& T,
+                            Eigen::SparseMatrix<double>& Dx, Eigen::SparseMatrix<double>& Dy, Eigen::SparseMatrix<double>& Dz, bool uniform) {
+  using namespace Eigen;
+  assert(T.cols() == 4);
+  const int n = V.rows(); int m = T.rows();
+
+  /*
+      F = [ ...
+      T(:,1) T(:,2) T(:,3); ...
+      T(:,1) T(:,3) T(:,4); ...
+      T(:,1) T(:,4) T(:,2); ...
+      T(:,2) T(:,4) T(:,3)]; */
+  MatrixXi F(4*m,3);
+  for (int i = 0; i < m; i++) {
+    F.row(0*m + i) << T(i,0), T(i,1), T(i,2);
+    F.row(1*m + i) << T(i,0), T(i,2), T(i,3);
+    F.row(2*m + i) << T(i,0), T(i,3), T(i,1);
+    F.row(3*m + i) << T(i,1), T(i,3), T(i,2);
+  }
+  // compute volume of each tet
+  VectorXd vol; igl::volume(V,T,vol);
+
+  VectorXd A(F.rows());
+  MatrixXd N(F.rows(),3);
+  if (!uniform) {
+    // compute tetrahedron face normals
+    igl::per_face_normals(V,F,N); int norm_rows = N.rows();
+    for (int i = 0; i < norm_rows; i++)
+      N.row(i) /= N.row(i).norm();
+    igl::doublearea(V,F,A); A/=2.;
+  } else {
+    // Use a uniform tetrahedra as a reference:
+    //      V = h*[0,0,0;1,0,0;0.5,sqrt(3)/2.,0;0.5,sqrt(3)/6.,sqrt(2)/sqrt(3)]
+    //
+    // With normals
+    //         0         0    1.0000
+    //         0.8165   -0.4714   -0.3333
+    //         0          0.9428   -0.3333
+    //         -0.8165   -0.4714   -0.3333
+    for (int i = 0; i < m; i++) {
+      N.row(0*m+i) << 0,0,1;
+      double a = sqrt(2)*std::cbrt(3*vol(i));
+      A(0*m+i) = (pow(a,2)*sqrt(3))/4.;
+    }
+    for (int i = 0; i < m; i++) {
+      N.row(1*m+i) << 0.8165,-0.4714,-0.3333;
+      double a = sqrt(2)*std::cbrt(3*vol(i));
+      A(1*m+i) = (pow(a,2)*sqrt(3))/4.;
+    }
+    for (int i = 0; i < m; i++) {
+      N.row(2*m+i) << 0,0.9428,-0.3333;
+      double a = sqrt(2)*std::cbrt(3*vol(i));
+      A(2*m+i) = (pow(a,2)*sqrt(3))/4.;
+    }
+    for (int i = 0; i < m; i++) {
+      N.row(3*m+i) << -0.8165,-0.4714,-0.3333;
+      double a = sqrt(2)*std::cbrt(3*vol(i));
+      A(3*m+i) = (pow(a,2)*sqrt(3))/4.;
+    }
+
+  }
+
+  /*  G = sparse( ...
+      [0*m + repmat(1:m,1,4) ...
+       1*m + repmat(1:m,1,4) ...
+       2*m + repmat(1:m,1,4)], ...
+      repmat([T(:,4);T(:,2);T(:,3);T(:,1)],3,1), ...
+      repmat(A./(3*repmat(vol,4,1)),3,1).*N(:), ...
+      3*m,n);*/
+  std::vector<Triplet<double> > Dx_t,Dy_t,Dz_t;
+  for (int i = 0; i < 4*m; i++) {
+    int T_j; // j indexes : repmat([T(:,4);T(:,2);T(:,3);T(:,1)],3,1)
+    switch (i/m) {
+      case 0:
+        T_j = 3;
+        break;
+      case 1:
+        T_j = 1;
+        break;
+      case 2:
+        T_j = 2;
+        break;
+      case 3:
+        T_j = 0;
+        break;
+    }
+    int i_idx = i%m;
+    int j_idx = T(i_idx,T_j);
+
+    double val_before_n = A(i)/(3*vol(i_idx));
+    Dx_t.push_back(Triplet<double>(i_idx, j_idx, val_before_n * N(i,0)));
+    Dy_t.push_back(Triplet<double>(i_idx, j_idx, val_before_n * N(i,1)));
+    Dz_t.push_back(Triplet<double>(i_idx, j_idx, val_before_n * N(i,2)));
+  }
+  Dx.resize(m,n); Dy.resize(m,n); Dz.resize(m,n);
+  Dx.setFromTriplets(Dx_t.begin(), Dx_t.end());
+  Dy.setFromTriplets(Dy_t.begin(), Dy_t.end());
+  Dz.setFromTriplets(Dz_t.begin(), Dz_t.end());
+
+}
+
+// Computes the weights and solve the linear system for the quadratic proxy specified in the paper
+// The output of this is used to generate a search direction that will be fed to the Linesearch class
+class WeightedGlobalLocal {
+
+public:
+  WeightedGlobalLocal(igl::SLIMData& state);
+
+  // Compute necessary information before solving the proxy quadratic
+  void pre_calc();
+
+  // Solve the weighted proxy global step
+  // Output:
+  //    V_new #V by dim list of mesh positions (will be fed to a linesearch algorithm)
+  void solve_weighted_proxy(Eigen::MatrixXd& V_new);
+
+  // Compute the energy specified in the SLIMData structure + the soft constraint energy (in case there are soft constraints)
+  // Input:
+  //    V_new #V by dim list of mesh positions
+  virtual double compute_energy(Eigen::MatrixXd& V_new);
+
+//private:
+
+  void compute_jacobians(const Eigen::MatrixXd& V_o);
+  double compute_energy_with_jacobians(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F,
+  const Eigen::MatrixXd& Ji, Eigen::MatrixXd& V_o, Eigen::VectorXd& areas);
+  double compute_soft_const_energy(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F,
+                                             Eigen::MatrixXd& V_o);
+
+  void update_weights_and_closest_rotations(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F, Eigen::MatrixXd& uv);
+  void solve_weighted_arap(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F, Eigen::MatrixXd& uv, Eigen::VectorXi& b,
+      Eigen::MatrixXd& bc);
+
+  void build_linear_system(Eigen::SparseMatrix<double> &L);
+  void buildA(Eigen::SparseMatrix<double>& A);
+  void buildRhs(const Eigen::SparseMatrix<double>& At);
+
+  void add_soft_constraints(Eigen::SparseMatrix<double> &L);
+  void add_proximal_penalty();
+
+  igl::SLIMData& m_state;
+  Eigen::VectorXd 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;
+};
+
+//// Implementation
+
+WeightedGlobalLocal::WeightedGlobalLocal(igl::SLIMData& state) :
+                                  m_state(state) {
+}
+
+void WeightedGlobalLocal::solve_weighted_proxy(Eigen::MatrixXd& V_new) {
+
+  update_weights_and_closest_rotations(m_state.V,m_state.F,V_new);
+  solve_weighted_arap(m_state.V,m_state.F,V_new,m_state.b,m_state.bc);
+}
+
+void WeightedGlobalLocal::compute_jacobians(const Eigen::MatrixXd& uv) {
+  if (m_state.F.cols() == 3){
+    // Ji=[D1*u,D2*u,D1*v,D2*v];
+    Ji.col(0) = Dx*uv.col(0); Ji.col(1) = Dy*uv.col(0);
+    Ji.col(2) = Dx*uv.col(1); Ji.col(3) = Dy*uv.col(1);
+  } else /*tet mesh*/{
+    // Ji=[D1*u,D2*u,D3*u, D1*v,D2*v, D3*v, D1*w,D2*w,D3*w];
+    Ji.col(0) = Dx*uv.col(0); Ji.col(1) = Dy*uv.col(0); Ji.col(2) = Dz*uv.col(0);
+    Ji.col(3) = Dx*uv.col(1); Ji.col(4) = Dy*uv.col(1); Ji.col(5) = Dz*uv.col(1);
+    Ji.col(6) = Dx*uv.col(2); Ji.col(7) = Dy*uv.col(2); Ji.col(8) = Dz*uv.col(2);
+  }
+}
+
+void WeightedGlobalLocal::update_weights_and_closest_rotations(const Eigen::MatrixXd& V,
+       const Eigen::MatrixXi& F, Eigen::MatrixXd& uv) {
+  compute_jacobians(uv);
+
+  const double eps = 1e-8;
+  double exp_f = m_state.exp_factor;
+
+  if (dim==2) {
+    for(int i=0; i <Ji.rows(); ++i ) {
+    typedef Eigen::Matrix<double,2,2> Mat2;
+    typedef Eigen::Matrix<double,2,1> Vec2;
+    Mat2 ji,ri,ti,ui,vi; Vec2 sing; Vec2 closest_sing_vec;Mat2 mat_W;
+    Vec2 m_sing_new;
+    double s1,s2;
+
+    ji(0,0) = Ji(i,0); ji(0,1) = Ji(i,1);
+    ji(1,0) = Ji(i,2); ji(1,1) = Ji(i,3);
+
+    igl::polar_svd(ji,ri,ti,ui,sing,vi);
+
+    s1 = sing(0); s2 = sing(1);
+
+    // Update Weights according to energy
+    switch(m_state.slim_energy) {
+    case igl::SLIMData::ARAP: {
+      m_sing_new << 1,1;
+      break;
+    } case igl::SLIMData::SYMMETRIC_DIRICHLET: {
+        double s1_g = 2* (s1-pow(s1,-3));
+        double s2_g = 2 * (s2-pow(s2,-3));
+        m_sing_new << sqrt(s1_g/(2*(s1-1))), sqrt(s2_g/(2*(s2-1)));
+        break;
+    } case igl::SLIMData::LOG_ARAP: {
+        double s1_g = 2 * (log(s1)/s1);
+        double s2_g = 2 * (log(s2)/s2);
+        m_sing_new << sqrt(s1_g/(2*(s1-1))), sqrt(s2_g/(2*(s2-1)));
+        break;
+    } case igl::SLIMData::CONFORMAL: {
+        double s1_g = 1/(2*s2) - s2/(2*pow(s1,2));
+        double s2_g = 1/(2*s1) - s1/(2*pow(s2,2));
+
+        double geo_avg = sqrt(s1*s2);
+        double s1_min = geo_avg; double s2_min = geo_avg;
+
+        m_sing_new << sqrt(s1_g/(2*(s1-s1_min))), sqrt(s2_g/(2*(s2-s2_min)));
+
+        // change local step
+        closest_sing_vec << s1_min,s2_min;
+        ri = ui*closest_sing_vec.asDiagonal()*vi.transpose();
+        break;
+    } case igl::SLIMData::EXP_CONFORMAL: {
+        double s1_g = 2* (s1-pow(s1,-3));
+        double s2_g = 2 * (s2-pow(s2,-3));
+
+        double geo_avg = sqrt(s1*s2);
+        double s1_min = geo_avg; double s2_min = geo_avg;
+
+        double in_exp = exp_f*((pow(s1,2)+pow(s2,2))/(2*s1*s2));
+        double exp_thing = exp(in_exp);
+
+        s1_g *= exp_thing*exp_f;
+        s2_g *= exp_thing*exp_f;
+
+        m_sing_new << sqrt(s1_g/(2*(s1-1))), sqrt(s2_g/(2*(s2-1)));
+        break;
+    } case igl::SLIMData::EXP_SYMMETRIC_DIRICHLET: {
+        double s1_g = 2* (s1-pow(s1,-3));
+        double s2_g = 2 * (s2-pow(s2,-3));
+
+        double in_exp = exp_f*(pow(s1,2)+pow(s1,-2)+pow(s2,2)+pow(s2,-2));
+        double exp_thing = exp(in_exp);
+
+        s1_g *= exp_thing*exp_f;
+        s2_g *= exp_thing*exp_f;
+
+        m_sing_new << sqrt(s1_g/(2*(s1-1))), sqrt(s2_g/(2*(s2-1)));
+        break;
+      }
+    }
+
+    if (abs(s1-1) < eps) m_sing_new(0) = 1; if (abs(s2-1) < eps) m_sing_new(1) = 1;
+    mat_W = ui*m_sing_new.asDiagonal()*ui.transpose();
+
+    W_11(i) = mat_W(0,0); W_12(i) = mat_W(0,1); W_21(i) = mat_W(1,0); W_22(i) = mat_W(1,1);
+
+    // 2) Update local step (doesn't have to be a rotation, for instance in case of conformal energy)
+    Ri(i,0) = ri(0,0); Ri(i,1) = ri(1,0); Ri(i,2) = ri(0,1); Ri(i,3) = ri(1,1);
+   }
+  } else {
+    typedef Eigen::Matrix<double,3,1> Vec3; typedef Eigen::Matrix<double,3,3> Mat3;
+    Mat3 ji; Vec3 m_sing_new; Vec3 closest_sing_vec;
+    const double sqrt_2 = sqrt(2);
+    for(int i=0; i <Ji.rows(); ++i ) {
+      ji(0,0) = Ji(i,0); ji(0,1) = Ji(i,1); ji(0,2) = Ji(i,2);
+      ji(1,0) = Ji(i,3); ji(1,1) = Ji(i,4); ji(1,2) = Ji(i,5);
+      ji(2,0) = Ji(i,6); ji(2,1) = Ji(i,7); ji(2,2) = Ji(i,8);
+
+      Mat3 ri,ti,ui,vi;
+      Vec3 sing;
+      igl::polar_svd(ji,ri,ti,ui,sing,vi);
+
+      double s1 = sing(0); double s2 = sing(1); double s3 = sing(2);
+
+      // 1) Update Weights
+      switch(m_state.slim_energy) {
+        case igl::SLIMData::ARAP: {
+          m_sing_new << 1,1,1;
+          break;
+        } case igl::SLIMData::LOG_ARAP: {
+            double s1_g = 2 * (log(s1)/s1);
+            double s2_g = 2 * (log(s2)/s2);
+            double s3_g = 2 * (log(s3)/s3);
+            m_sing_new << sqrt(s1_g/(2*(s1-1))), sqrt(s2_g/(2*(s2-1))), sqrt(s3_g/(2*(s3-1)));
+            break;
+          } case igl::SLIMData::SYMMETRIC_DIRICHLET: {
+            double s1_g = 2* (s1-pow(s1,-3));
+            double s2_g = 2 * (s2-pow(s2,-3));
+            double s3_g = 2 * (s3-pow(s3,-3));
+            m_sing_new << sqrt(s1_g/(2*(s1-1))), sqrt(s2_g/(2*(s2-1))), sqrt(s3_g/(2*(s3-1)));
+            break;
+          } case igl::SLIMData::EXP_SYMMETRIC_DIRICHLET: {
+           double s1_g = 2* (s1-pow(s1,-3));
+          double s2_g = 2 * (s2-pow(s2,-3));
+          double s3_g = 2 * (s3-pow(s3,-3));
+          m_sing_new << sqrt(s1_g/(2*(s1-1))), sqrt(s2_g/(2*(s2-1))), sqrt(s3_g/(2*(s3-1)));
+
+          double in_exp = exp_f*(pow(s1,2)+pow(s1,-2)+pow(s2,2)+pow(s2,-2)+pow(s3,2)+pow(s3,-2));
+          double exp_thing = exp(in_exp);
+
+          s1_g *= exp_thing*exp_f;
+          s2_g *= exp_thing*exp_f;
+          s3_g *= exp_thing*exp_f;
+
+          m_sing_new << sqrt(s1_g/(2*(s1-1))), sqrt(s2_g/(2*(s2-1))), sqrt(s3_g/(2*(s3-1)));
+
+          break;
+        }
+        case igl::SLIMData::CONFORMAL: {
+          double common_div = 9*(pow(s1*s2*s3,5./3.));
+
+          double s1_g = (-2*s2*s3*(pow(s2,2)+pow(s3,2)-2*pow(s1,2)) ) / common_div;
+          double s2_g = (-2*s1*s3*(pow(s1,2)+pow(s3,2)-2*pow(s2,2)) ) / common_div;
+          double s3_g = (-2*s1*s2*(pow(s1,2)+pow(s2,2)-2*pow(s3,2)) ) / common_div;
+
+          double closest_s = sqrt(pow(s1,2)+pow(s3,2)) / sqrt_2;
+          double s1_min = closest_s; double s2_min = closest_s; double s3_min = closest_s;
+
+          m_sing_new << sqrt(s1_g/(2*(s1-s1_min))), sqrt(s2_g/(2*(s2-s2_min))), sqrt(s3_g/(2*(s3-s3_min)));
+
+          // change local step
+          closest_sing_vec << s1_min,s2_min,s3_min;
+          ri = ui*closest_sing_vec.asDiagonal()*vi.transpose();
+          break;
+        }
+        case igl::SLIMData::EXP_CONFORMAL: {
+          // E_conf = (s1^2 + s2^2 + s3^2)/(3*(s1*s2*s3)^(2/3) )
+          // dE_conf/ds1 = (-2*(s2*s3)*(s2^2+s3^2 -2*s1^2) ) / (9*(s1*s2*s3)^(5/3))
+          // Argmin E_conf(s1): s1 = sqrt(s1^2+s2^2)/sqrt(2)
+          double common_div = 9*(pow(s1*s2*s3,5./3.));
+
+          double s1_g = (-2*s2*s3*(pow(s2,2)+pow(s3,2)-2*pow(s1,2)) ) / common_div;
+          double s2_g = (-2*s1*s3*(pow(s1,2)+pow(s3,2)-2*pow(s2,2)) ) / common_div;
+          double s3_g = (-2*s1*s2*(pow(s1,2)+pow(s2,2)-2*pow(s3,2)) ) / common_div;
+
+          double in_exp = exp_f*( (pow(s1,2)+pow(s2,2)+pow(s3,2))/ (3*pow((s1*s2*s3),2./3)) ); ;
+          double exp_thing = exp(in_exp);
+
+          double closest_s = sqrt(pow(s1,2)+pow(s3,2)) / sqrt_2;
+          double s1_min = closest_s; double s2_min = closest_s; double s3_min = closest_s;
+
+          s1_g *= exp_thing*exp_f;
+          s2_g *= exp_thing*exp_f;
+          s3_g *= exp_thing*exp_f;
+
+          m_sing_new << sqrt(s1_g/(2*(s1-s1_min))), sqrt(s2_g/(2*(s2-s2_min))), sqrt(s3_g/(2*(s3-s3_min)));
+
+          // change local step
+          closest_sing_vec << s1_min,s2_min,s3_min;
+          ri = ui*closest_sing_vec.asDiagonal()*vi.transpose();
+        }
+      }
+      if (abs(s1-1) < eps) m_sing_new(0) = 1; if (abs(s2-1) < eps) m_sing_new(1) = 1; if (abs(s3-1) < eps) m_sing_new(2) = 1;
+      Mat3 mat_W;
+      mat_W = ui*m_sing_new.asDiagonal()*ui.transpose();
+
+      W_11(i) = mat_W(0,0);
+      W_12(i) = mat_W(0,1);
+      W_13(i) = mat_W(0,2);
+      W_21(i) = mat_W(1,0);
+      W_22(i) = mat_W(1,1);
+      W_23(i) = mat_W(1,2);
+      W_31(i) = mat_W(2,0);
+      W_32(i) = mat_W(2,1);
+      W_33(i) = mat_W(2,2);
+
+      // 2) Update closest rotations (not rotations in case of conformal energy)
+      Ri(i,0) = ri(0,0); Ri(i,1) = ri(1,0); Ri(i,2) = ri(2,0);
+      Ri(i,3) = ri(0,1); Ri(i,4) = ri(1,1); Ri(i,5) = ri(2,1);
+      Ri(i,6) = ri(0,2); Ri(i,7) = ri(1,2); Ri(i,8) = ri(2,2);
+    } // for loop end
+
+  } // if dim end
+
+}
+
+void WeightedGlobalLocal::solve_weighted_arap(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F,
+        Eigen::MatrixXd& uv, Eigen::VectorXi& soft_b_p, Eigen::MatrixXd& soft_bc_p) {
+  using namespace Eigen;
+
+  Eigen::SparseMatrix<double> L;
+  build_linear_system(L);
+
+  // solve
+  Eigen::VectorXd Uc;
+  if (dim == 2) {
+    SimplicialLDLT<SparseMatrix<double> > solver;
+    Uc = solver.compute(L).solve(rhs);
+  } else { // seems like CG performs much worse for 2D and way better for 3D
+    Eigen::VectorXd guess(uv.rows()*dim);
+    for (int i = 0; i < dim; i++) for (int j = 0; j < dim; j++) guess(uv.rows()*i + j) = uv(i,j); // flatten vector
+    ConjugateGradient<SparseMatrix<double>, Eigen::Upper> solver;
+    solver.setTolerance(1e-8);
+    Uc = solver.compute(L).solveWithGuess(rhs,guess);
+  }
+
+  for (int i = 0; i < dim; i++)
+    uv.col(i) = Uc.block(i*v_n,0,v_n,1);
+}
+
+void WeightedGlobalLocal::pre_calc() {
+  if (!has_pre_calc) {
+    v_n = m_state.v_num; f_n = m_state.f_num;
+
+    if (m_state.F.cols() == 3) {
+      dim = 2;
+      Eigen::MatrixXd F1,F2,F3;
+      igl::local_basis(m_state.V,m_state.F,F1,F2,F3);
+      compute_surface_gradient_matrix(m_state.V,m_state.F,F1,F2,Dx,Dy);
+
+      W_11.resize(f_n); W_12.resize(f_n); W_21.resize(f_n); W_22.resize(f_n);
+    } else {
+      dim = 3;
+      compute_tet_grad_matrix(m_state.V,m_state.F,Dx,Dy,Dz,
+        m_state.mesh_improvement_3d /*use normal gradient, or one from a "regular" tet*/);
+
+      W_11.resize(f_n);W_12.resize(f_n);W_13.resize(f_n);
+      W_21.resize(f_n);W_22.resize(f_n);W_23.resize(f_n);
+      W_31.resize(f_n);W_32.resize(f_n);W_33.resize(f_n);
+    }
+
+    Dx.makeCompressed();Dy.makeCompressed(); Dz.makeCompressed();
+    Ri.resize(f_n, dim*dim); Ji.resize(f_n, dim*dim);
+    rhs.resize(dim*m_state.v_num);
+
+    // flattened weight matrix
+    M.resize(dim*dim*f_n);
+    for (int i = 0; i < dim*dim; i++)
+      for (int j = 0; j < f_n; j++)
+        M(i*f_n + j) = m_state.M(j);
+
+    first_solve = true;
+    has_pre_calc = true;
+  }
+}
+
+void WeightedGlobalLocal::build_linear_system(Eigen::SparseMatrix<double> &L) {
+  // formula (35) in paper
+  Eigen::SparseMatrix<double> A(dim*dim*f_n, dim*v_n);
+  buildA(A);
+
+  Eigen::SparseMatrix<double> At = A.transpose();
+  At.makeCompressed();
+
+  Eigen::SparseMatrix<double> id_m(At.rows(),At.rows()); id_m.setIdentity();
+
+  // add proximal penalty
+  L = At*M.asDiagonal()*A + m_state.proximal_p * id_m; //add also a proximal term
+  L.makeCompressed();
+
+  buildRhs(At);
+  Eigen::SparseMatrix<double> OldL = L;
+  add_soft_constraints(L);
+  L.makeCompressed();
+}
+
+void WeightedGlobalLocal::add_soft_constraints(Eigen::SparseMatrix<double> &L) {
+  int v_n = m_state.v_num;
+  for (int d = 0; d < dim; d++) {
+    for (int i = 0; i < m_state.b.rows(); i++) {
+      int v_idx = m_state.b(i);
+      rhs(d*v_n + v_idx) += m_state.soft_const_p * m_state.bc(i,d); // rhs
+      L.coeffRef(d*v_n + v_idx, d*v_n + v_idx) += m_state.soft_const_p; // diagonal of matrix
+    }
+  }
+}
+
+double WeightedGlobalLocal::compute_energy(Eigen::MatrixXd& V_new) {
+  compute_jacobians(V_new);
+  return compute_energy_with_jacobians(m_state.V,m_state.F, Ji, V_new,m_state.M) + compute_soft_const_energy(m_state.V,m_state.F,V_new);
+}
+
+double WeightedGlobalLocal::compute_soft_const_energy(const Eigen::MatrixXd& V,
+                                                       const Eigen::MatrixXi& F,
+                                                       Eigen::MatrixXd& V_o) {
+  double e = 0;
+  for (int i = 0; i < m_state.b.rows(); i++) {
+    e += m_state.soft_const_p*(m_state.bc.row(i)-V_o.row(m_state.b(i))).squaredNorm();
+  }
+  return e;
+}
+
+double WeightedGlobalLocal::compute_energy_with_jacobians(const Eigen::MatrixXd& V,
+       const Eigen::MatrixXi& F, const Eigen::MatrixXd& Ji, Eigen::MatrixXd& uv, Eigen::VectorXd& areas) {
+
+  double energy = 0;
+  if (dim == 2) {
+    Eigen::Matrix<double,2,2> ji;
+    for (int i = 0; i < f_n; i++) {
+      ji(0,0) = Ji(i,0); ji(0,1) = Ji(i,1);
+      ji(1,0) = Ji(i,2); ji(1,1) = Ji(i,3);
+
+      typedef Eigen::Matrix<double,2,2> Mat2;
+      typedef Eigen::Matrix<double,2,1> Vec2;
+      Mat2 ri,ti,ui,vi; Vec2 sing;
+      igl::polar_svd(ji,ri,ti,ui,sing,vi);
+      double s1 = sing(0); double s2 = sing(1);
+
+      switch(m_state.slim_energy) {
+        case igl::SLIMData::ARAP: {
+          energy+= areas(i) * (pow(s1-1,2) + pow(s2-1,2));
+          break;
+        }
+        case igl::SLIMData::SYMMETRIC_DIRICHLET: {
+          energy += areas(i) * (pow(s1,2) +pow(s1,-2) + pow(s2,2) + pow(s2,-2));
+          break;
+        }
+        case igl::SLIMData::EXP_SYMMETRIC_DIRICHLET: {
+          energy += areas(i) * exp(m_state.exp_factor*(pow(s1,2) +pow(s1,-2) + pow(s2,2) + pow(s2,-2)));
+          break;
+        }
+        case igl::SLIMData::LOG_ARAP: {
+          energy += areas(i) * (pow(log(s1),2) + pow(log(s2),2));
+          break;
+        }
+        case igl::SLIMData::CONFORMAL: {
+          energy += areas(i) * ( (pow(s1,2)+pow(s2,2))/(2*s1*s2) );
+          break;
+        }
+        case igl::SLIMData::EXP_CONFORMAL: {
+          energy += areas(i) * exp(m_state.exp_factor*((pow(s1,2)+pow(s2,2))/(2*s1*s2)));
+          break;
+        }
+
+      }
+
+    }
+  } else {
+    Eigen::Matrix<double,3,3> ji;
+    for (int i = 0; i < f_n; i++) {
+      ji(0,0) = Ji(i,0); ji(0,1) = Ji(i,1); ji(0,2) = Ji(i,2);
+      ji(1,0) = Ji(i,3); ji(1,1) = Ji(i,4); ji(1,2) = Ji(i,5);
+      ji(2,0) = Ji(i,6); ji(2,1) = Ji(i,7); ji(2,2) = Ji(i,8);
+
+      typedef Eigen::Matrix<double,3,3> Mat3;
+      typedef Eigen::Matrix<double,3,1> Vec3;
+      Mat3 ri,ti,ui,vi; Vec3 sing;
+      igl::polar_svd(ji,ri,ti,ui,sing,vi);
+      double s1 = sing(0); double s2 = sing(1); double s3 = sing(2);
+
+      switch(m_state.slim_energy) {
+        case igl::SLIMData::ARAP: {
+          energy+= areas(i) * (pow(s1-1,2) + pow(s2-1,2) + pow(s3-1,2));
+          break;
+        }
+        case igl::SLIMData::SYMMETRIC_DIRICHLET: {
+          energy += areas(i) * (pow(s1,2) +pow(s1,-2) + pow(s2,2) + pow(s2,-2) + pow(s3,2) + pow(s3,-2));
+          break;
+        }
+        case igl::SLIMData::EXP_SYMMETRIC_DIRICHLET: {
+          energy += areas(i) * exp(m_state.exp_factor*(pow(s1,2) +pow(s1,-2) + pow(s2,2) + pow(s2,-2) + pow(s3,2) + pow(s3,-2)));
+          break;
+        }
+        case igl::SLIMData::LOG_ARAP: {
+          energy += areas(i) * (pow(log(s1),2) + pow(log(abs(s2)),2) + pow(log(abs(s3)),2));
+          break;
+        }
+        case igl::SLIMData::CONFORMAL: {
+          energy += areas(i) * ( ( pow(s1,2)+pow(s2,2)+pow(s3,2) ) /(3*pow(s1*s2*s3,2./3.)) );
+          break;
+        }
+        case igl::SLIMData::EXP_CONFORMAL: {
+          energy += areas(i) * exp( ( pow(s1,2)+pow(s2,2)+pow(s3,2) ) /(3*pow(s1*s2*s3,2./3.)) );
+          break;
+        }
+      }
+    }
+  }
+
+  return energy;
+}
+
+void WeightedGlobalLocal::buildA(Eigen::SparseMatrix<double>& A) {
+  // formula (35) in paper
+  std::vector<Triplet<double> > IJV;
+  if (dim == 2) {
+    IJV.reserve(4*(Dx.outerSize()+ Dy.outerSize()));
+
+    /*A = [W11*Dx, W12*Dx;
+         W11*Dy, W12*Dy;
+         W21*Dx, W22*Dx;
+         W21*Dy, W22*Dy];*/
+    for (int k=0; k<Dx.outerSize(); ++k) {
+        for (SparseMatrix<double>::InnerIterator it(Dx,k); it; ++it) {
+          int dx_r = it.row();
+          int dx_c = it.col();
+          double val = it.value();
+
+          IJV.push_back(Triplet<double>(dx_r,dx_c, val*W_11(dx_r)));
+          IJV.push_back(Triplet<double>(dx_r,v_n + dx_c, val*W_12(dx_r)));
+
+          IJV.push_back(Triplet<double>(2*f_n+dx_r,dx_c, val*W_21(dx_r)));
+          IJV.push_back(Triplet<double>(2*f_n+dx_r,v_n + dx_c, val*W_22(dx_r)));
+      }
+    }
+
+    for (int k=0; k<Dy.outerSize(); ++k) {
+      for (SparseMatrix<double>::InnerIterator it(Dy,k); it; ++it) {
+        int dy_r = it.row();
+        int dy_c = it.col();
+        double val = it.value();
+
+        IJV.push_back(Triplet<double>(f_n+dy_r,dy_c, val*W_11(dy_r)));
+        IJV.push_back(Triplet<double>(f_n+dy_r,v_n + dy_c, val*W_12(dy_r)));
+
+        IJV.push_back(Triplet<double>(3*f_n+dy_r,dy_c, val*W_21(dy_r)));
+        IJV.push_back(Triplet<double>(3*f_n+dy_r,v_n + dy_c, val*W_22(dy_r)));
+      }
+    }
+  } else {
+
+    /*A = [W11*Dx, W12*Dx, W13*Dx;
+           W11*Dy, W12*Dy, W13*Dy;
+           W11*Dz, W12*Dz, W13*Dz;
+           W21*Dx, W22*Dx, W23*Dx;
+           W21*Dy, W22*Dy, W23*Dy;
+           W21*Dz, W22*Dz, W23*Dz;
+           W31*Dx, W32*Dx, W33*Dx;
+           W31*Dy, W32*Dy, W33*Dy;
+           W31*Dz, W32*Dz, W33*Dz;];*/
+    IJV.reserve(9*(Dx.outerSize()+ Dy.outerSize() + Dz.outerSize()));
+    for (int k = 0; k < Dx.outerSize(); k++) {
+      for (SparseMatrix<double>::InnerIterator it(Dx,k); it; ++it) {
+         int dx_r = it.row();
+         int dx_c = it.col();
+         double val = it.value();
+
+         IJV.push_back(Triplet<double>(dx_r,dx_c, val*W_11(dx_r)));
+         IJV.push_back(Triplet<double>(dx_r,v_n + dx_c, val*W_12(dx_r)));
+         IJV.push_back(Triplet<double>(dx_r,2*v_n + dx_c, val*W_13(dx_r)));
+
+         IJV.push_back(Triplet<double>(3*f_n+dx_r,dx_c, val*W_21(dx_r)));
+         IJV.push_back(Triplet<double>(3*f_n+dx_r,v_n + dx_c, val*W_22(dx_r)));
+         IJV.push_back(Triplet<double>(3*f_n+dx_r,2*v_n + dx_c, val*W_23(dx_r)));
+
+         IJV.push_back(Triplet<double>(6*f_n+dx_r,dx_c, val*W_31(dx_r)));
+         IJV.push_back(Triplet<double>(6*f_n+dx_r,v_n + dx_c, val*W_32(dx_r)));
+         IJV.push_back(Triplet<double>(6*f_n+dx_r,2*v_n + dx_c, val*W_33(dx_r)));
+      }
+    }
+
+    for (int k = 0; k < Dy.outerSize(); k++) {
+      for (SparseMatrix<double>::InnerIterator it(Dy,k); it; ++it) {
+         int dy_r = it.row();
+         int dy_c = it.col();
+         double val = it.value();
+
+         IJV.push_back(Triplet<double>(f_n+dy_r,dy_c, val*W_11(dy_r)));
+         IJV.push_back(Triplet<double>(f_n+dy_r,v_n + dy_c, val*W_12(dy_r)));
+         IJV.push_back(Triplet<double>(f_n+dy_r,2*v_n + dy_c, val*W_13(dy_r)));
+
+         IJV.push_back(Triplet<double>(4*f_n+dy_r,dy_c, val*W_21(dy_r)));
+         IJV.push_back(Triplet<double>(4*f_n+dy_r,v_n + dy_c, val*W_22(dy_r)));
+         IJV.push_back(Triplet<double>(4*f_n+dy_r,2*v_n + dy_c, val*W_23(dy_r)));
+
+         IJV.push_back(Triplet<double>(7*f_n+dy_r,dy_c, val*W_31(dy_r)));
+         IJV.push_back(Triplet<double>(7*f_n+dy_r,v_n + dy_c, val*W_32(dy_r)));
+         IJV.push_back(Triplet<double>(7*f_n+dy_r,2*v_n + dy_c, val*W_33(dy_r)));
+      }
+    }
+
+    for (int k = 0; k < Dz.outerSize(); k++) {
+      for (SparseMatrix<double>::InnerIterator it(Dz,k); it; ++it) {
+         int dz_r = it.row();
+         int dz_c = it.col();
+         double val = it.value();
+
+         IJV.push_back(Triplet<double>(2*f_n + dz_r,dz_c, val*W_11(dz_r)));
+         IJV.push_back(Triplet<double>(2*f_n + dz_r,v_n + dz_c, val*W_12(dz_r)));
+         IJV.push_back(Triplet<double>(2*f_n + dz_r,2*v_n + dz_c, val*W_13(dz_r)));
+
+         IJV.push_back(Triplet<double>(5*f_n+dz_r,dz_c, val*W_21(dz_r)));
+         IJV.push_back(Triplet<double>(5*f_n+dz_r,v_n + dz_c, val*W_22(dz_r)));
+         IJV.push_back(Triplet<double>(5*f_n+dz_r,2*v_n + dz_c, val*W_23(dz_r)));
+
+         IJV.push_back(Triplet<double>(8*f_n+dz_r,dz_c, val*W_31(dz_r)));
+         IJV.push_back(Triplet<double>(8*f_n+dz_r,v_n + dz_c, val*W_32(dz_r)));
+         IJV.push_back(Triplet<double>(8*f_n+dz_r,2*v_n + dz_c, val*W_33(dz_r)));
+      }
+    }
+  }
+  A.setFromTriplets(IJV.begin(),IJV.end());
+}
+
+void WeightedGlobalLocal::buildRhs(const Eigen::SparseMatrix<double>& At) {
+  VectorXd f_rhs(dim*dim*f_n); f_rhs.setZero();
+  if (dim==2) {
+    /*b = [W11*R11 + W12*R21; (formula (36))
+         W11*R12 + W12*R22;
+         W21*R11 + W22*R21;
+         W21*R12 + W22*R22];*/
+    for (int i = 0; i < f_n; i++) {
+      f_rhs(i+0*f_n) = W_11(i) * Ri(i,0) + W_12(i)*Ri(i,1);
+      f_rhs(i+1*f_n) = W_11(i) * Ri(i,2) + W_12(i)*Ri(i,3);
+      f_rhs(i+2*f_n) = W_21(i) * Ri(i,0) + W_22(i)*Ri(i,1);
+      f_rhs(i+3*f_n) = W_21(i) * Ri(i,2) + W_22(i)*Ri(i,3);
+    }
+  } else {
+    /*b = [W11*R11 + W12*R21 + W13*R31;
+         W11*R12 + W12*R22 + W13*R32;
+         W11*R13 + W12*R23 + W13*R33;
+         W21*R11 + W22*R21 + W23*R31;
+         W21*R12 + W22*R22 + W23*R32;
+         W21*R13 + W22*R23 + W23*R33;
+         W31*R11 + W32*R21 + W33*R31;
+         W31*R12 + W32*R22 + W33*R32;
+         W31*R13 + W32*R23 + W33*R33;];*/
+    for (int i = 0; i < f_n; i++) {
+      f_rhs(i+0*f_n) = W_11(i) * Ri(i,0) + W_12(i)*Ri(i,1) + W_13(i)*Ri(i,2);
+      f_rhs(i+1*f_n) = W_11(i) * Ri(i,3) + W_12(i)*Ri(i,4) + W_13(i)*Ri(i,5);
+      f_rhs(i+2*f_n) = W_11(i) * Ri(i,6) + W_12(i)*Ri(i,7) + W_13(i)*Ri(i,8);
+      f_rhs(i+3*f_n) = W_21(i) * Ri(i,0) + W_22(i)*Ri(i,1) + W_23(i)*Ri(i,2);
+      f_rhs(i+4*f_n) = W_21(i) * Ri(i,3) + W_22(i)*Ri(i,4) + W_23(i)*Ri(i,5);
+      f_rhs(i+5*f_n) = W_21(i) * Ri(i,6) + W_22(i)*Ri(i,7) + W_23(i)*Ri(i,8);
+      f_rhs(i+6*f_n) = W_31(i) * Ri(i,0) + W_32(i)*Ri(i,1) + W_33(i)*Ri(i,2);
+      f_rhs(i+7*f_n) = W_31(i) * Ri(i,3) + W_32(i)*Ri(i,4) + W_33(i)*Ri(i,5);
+      f_rhs(i+8*f_n) = W_31(i) * Ri(i,6) + W_32(i)*Ri(i,7) + W_33(i)*Ri(i,8);
+    }
+  }
+  VectorXd uv_flat(dim*v_n);
+  for (int i = 0; i < dim; i++)
+    for (int j = 0; j < v_n; j++)
+      uv_flat(v_n*i+j) = m_state.V_o(j,i);
+
+  rhs = (At*M.asDiagonal()*f_rhs + m_state.proximal_p * uv_flat);
+}
+
+
+
+// #define TwoPi  6.28318530717958648
+// const double eps=1e-14;
+
+/// Slim Implementation
+
+IGL_INLINE void igl::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) {
+
+  data.V = V;
+  data.F = F;
+  data.V_o = V_init;
+
+  data.v_num = V.rows();
+  data.f_num = F.rows();
+
+  data.slim_energy = slim_energy;
+
+  data.b = b;
+  data.bc = bc;
+  data.soft_const_p = soft_p;
+
+  data.proximal_p = 0.0001;
+
+  igl::doublearea(V,F,data.M); data.M /= 2.;
+  data.mesh_area = data.M.sum();
+  data.mesh_improvement_3d = false; // whether to use a jacobian derived from a real mesh or an abstract regular mesh (used for mesh improvement)
+  data.exp_factor = 1.0; // param used only for exponential energies (e.g exponential symmetric dirichlet)
+
+  assert (F.cols() == 3 || F.cols() == 4);
+  data.wGlobalLocal = new WeightedGlobalLocal(data);
+
+  data.wGlobalLocal->pre_calc();
+  data.energy = data.wGlobalLocal->compute_energy(data.V_o)/data.mesh_area;
+}
+
+IGL_INLINE void igl::slim_solve(SLIMData& data, int iter_num) {
+
+  for (int i = 0; i < iter_num; i++) {
+    Eigen::MatrixXd dest_res;
+    dest_res = data.V_o;
+    data.wGlobalLocal->solve_weighted_proxy(dest_res);
+
+    double old_energy = data.energy;
+
+    std::function<double(Eigen::MatrixXd&)> compute_energy = [&](Eigen::MatrixXd& aaa) { return data.wGlobalLocal->compute_energy(aaa); };
+
+    data.energy = igl::flip_avoiding_line_search(data.F,data.V_o, dest_res, compute_energy,
+                                           data.energy*data.mesh_area)/data.mesh_area;
+  }
+}

include/igl/slim.h

@@ -0,0 +1,89 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Michael Rabinovich
+//
+// 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/.
+#ifndef SLIM_H
+#define SLIM_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 {
+
+  // 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 {
+    ARAP,
+    LOG_ARAP,
+    SYMMETRIC_DIRICHLET,
+    CONFORMAL,
+    EXP_CONFORMAL,
+    EXP_SYMMETRIC_DIRICHLET
+  };
+  SLIM_ENERGY slim_energy;
+
+  // Optional Input
+    // 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
+
+  // Output
+  Eigen::MatrixXd V_o; // #V by dim list of mesh vertex positions (dim = 2 for parametrization, 3 otherwise)
+  double energy; // objective value
+
+  // INTERNAL
+  Eigen::VectorXd M;
+  double mesh_area;
+  double avg_edge_length;
+  int v_num;
+  int f_num;
+  double proximal_p;
+
+  WeightedGlobalLocal* wGlobalLocal;
+};
+
+  // 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);
+
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "slim.cpp"
+#endif
+
+
+#endif // SLIM_H

tutorial/611_SLIM/CMakeLists.txt

@@ -0,0 +1,8 @@
+cmake_minimum_required(VERSION 2.8.12)
+project(611_SLIM)
+
+add_executable(${PROJECT_NAME}_bin
+  main.cpp)
+target_include_directories(${PROJECT_NAME}_bin PRIVATE ${LIBIGL_INCLUDE_DIRS})
+target_compile_definitions(${PROJECT_NAME}_bin PRIVATE ${LIBIGL_DEFINITIONS})
+target_link_libraries(${PROJECT_NAME}_bin ${LIBIGL_LIBRARIES} ${LIBIGL_VIEWER_EXTRA_LIBRARIES})

tutorial/611_SLIM/main.cpp

@@ -0,0 +1,394 @@
+#include <iostream>
+
+#include "igl/slim.h"
+
+#include "igl/components.h"
+#include "igl/readOBJ.h"
+#include "igl/writeOBJ.h"
+#include "igl/Timer.h"
+
+#include "igl/boundary_loop.h"
+#include "igl/map_vertices_to_circle.h"
+#include "igl/harmonic.h"
+#include <igl/serialize.h>
+#include <igl/read_triangle_mesh.h>
+#include <igl/viewer/Viewer.h>
+#include <igl/flipped_triangles.h>
+#include <igl/euler_characteristic.h>
+#include <igl/barycenter.h>
+#include <igl/adjacency_matrix.h>
+#include <igl/is_edge_manifold.h>
+#include <igl/doublearea.h>
+#include <igl/cat.h>
+
+#include <stdlib.h>
+
+#include <string>
+#include <vector>
+
+using namespace std;
+using namespace Eigen;
+
+void check_mesh_for_issues(Eigen::MatrixXd& V, Eigen::MatrixXi& F);
+void param_2d_demo_iter(igl::viewer::Viewer& viewer);
+void get_soft_constraint_for_circle(Eigen::MatrixXd& V_o, Eigen::MatrixXi& F, Eigen::VectorXi& b, Eigen::MatrixXd& bc);
+void soft_const_demo_iter(igl::viewer::Viewer& viewer);
+void deform_3d_demo_iter(igl::viewer::Viewer& viewer);
+void get_cube_corner_constraints(Eigen::MatrixXd& V_o, Eigen::MatrixXi& F, Eigen::VectorXi& b, Eigen::MatrixXd& bc);
+void display_3d_mesh(igl::viewer::Viewer& viewer);
+void int_set_to_eigen_vector(const std::set<int>& int_set, Eigen::VectorXi& vec);
+
+Eigen::MatrixXd V;
+Eigen::MatrixXi F;
+bool first_iter = true;
+igl::SLIMData sData;
+igl::Timer timer;
+
+double uv_scale_param;
+
+enum DEMO_TYPE {
+  PARAM_2D,
+  SOFT_CONST,
+  DEFORM_3D
+};
+DEMO_TYPE demo_type;
+
+bool key_down(igl::viewer::Viewer& viewer, unsigned char key, int modifier){
+  if (key == ' ') {
+    switch (demo_type) {
+      case PARAM_2D: {
+        param_2d_demo_iter(viewer);
+        break;
+      }
+      case SOFT_CONST: {
+        soft_const_demo_iter(viewer);
+        break;
+      }
+      case DEFORM_3D: {
+        deform_3d_demo_iter(viewer);
+        break;
+      }
+      default:
+        break;
+    }
+  }
+
+  return false;
+}
+
+void param_2d_demo_iter(igl::viewer::Viewer& viewer) {
+  if (first_iter) {
+    timer.start();
+    igl::read_triangle_mesh(TUTORIAL_SHARED_PATH "/face.obj", V, F);
+    check_mesh_for_issues(V,F);
+    cout << "\tMesh is valid!" << endl;
+
+    Eigen::MatrixXd uv_init;
+    Eigen::VectorXi bnd; Eigen::MatrixXd bnd_uv;
+    igl::boundary_loop(F,bnd);
+    igl::map_vertices_to_circle(V,bnd,bnd_uv);
+
+    igl::harmonic(V,F,bnd,bnd_uv,1,uv_init);
+    if (igl::flipped_triangles(uv_init,F).size() != 0) {
+      igl::harmonic(F,bnd,bnd_uv,1,uv_init); // use uniform laplacian
+    }
+
+    cout << "initialized parametrization" << endl;
+
+    sData.slim_energy = igl::SLIMData::SYMMETRIC_DIRICHLET;
+    Eigen::VectorXi b; Eigen::MatrixXd bc;
+    slim_precompute(V,F,uv_init,sData, igl::SLIMData::SYMMETRIC_DIRICHLET, b,bc,0);
+
+    uv_scale_param = 15 * (1./sqrt(sData.mesh_area));
+    viewer.data.set_mesh(V, F);
+    viewer.core.align_camera_center(V,F);
+    viewer.data.set_uv(sData.V_o*uv_scale_param);
+    viewer.data.compute_normals();
+    viewer.core.show_texture = true;
+
+    first_iter = false;
+  } else {
+    slim_solve(sData,1); // 1 iter
+    viewer.data.set_uv(sData.V_o*uv_scale_param);
+    cout << "time = " << timer.getElapsedTime() << endl;
+  }
+}
+
+void soft_const_demo_iter(igl::viewer::Viewer& viewer) {
+  if (first_iter) {
+
+    igl::read_triangle_mesh(TUTORIAL_SHARED_PATH "/circle.obj", V, F);
+
+    check_mesh_for_issues(V,F);
+    cout << "\tMesh is valid!" << endl;
+    Eigen::MatrixXd V_0 = V.block(0,0,V.rows(),2);
+
+    Eigen::VectorXi b; Eigen::MatrixXd bc;
+    get_soft_constraint_for_circle(V_0,F,b,bc);
+    double soft_const_p = 1e5;
+    slim_precompute(V,F,V_0,sData,igl::SLIMData::SYMMETRIC_DIRICHLET,b,bc,soft_const_p);
+
+    viewer.data.set_mesh(V, F);
+    viewer.core.align_camera_center(V,F);
+    viewer.data.compute_normals();
+    viewer.core.show_lines = true;
+
+    first_iter = false;
+
+  } else {
+    slim_solve(sData,1); // 1 iter
+    viewer.data.set_mesh(sData.V_o, F);
+  }
+}
+
+void deform_3d_demo_iter(igl::viewer::Viewer& viewer) {
+  if (first_iter) {
+    igl::readOBJ(TUTORIAL_SHARED_PATH "/cube_40k.obj", V, F);
+
+    Eigen::MatrixXd V_0 = V;
+    Eigen::VectorXi b; Eigen::MatrixXd bc;
+    get_cube_corner_constraints(V_0,F,b,bc);
+
+    double soft_const_p = 1e5;
+    sData.exp_factor = 5.0;
+    slim_precompute(V,F,V_0,sData,igl::SLIMData::EXP_CONFORMAL,b,bc,soft_const_p);
+    cout << "precomputed" << endl;
+
+    first_iter = false;
+    display_3d_mesh(viewer);
+
+  } else {
+    slim_solve(sData,1); // 1 iter
+    display_3d_mesh(viewer);
+  }
+}
+
+void display_3d_mesh(igl::viewer::Viewer& viewer) {
+  MatrixXd V_temp; MatrixXi F_temp;
+  Eigen::MatrixXd Barycenters;
+
+  igl::barycenter(sData.V,sData.F,Barycenters);
+  //cout << "Barycenters.rows() = " << Barycenters.rows() << endl;
+  //double t = double((key - '1')+1) / 9.0;
+  double view_depth = 10.;
+  double t = view_depth/9.;
+
+  VectorXd v = Barycenters.col(2).array() - Barycenters.col(2).minCoeff();
+  v /= v.col(0).maxCoeff();
+
+  vector<int> s;
+
+  for (unsigned i=0; i<v.size();++i)
+    if (v(i) < t)
+      s.push_back(i);
+
+  V_temp.resize(s.size()*4,3);
+  F_temp.resize(s.size()*4,3);
+
+  for (unsigned i=0; i<s.size();++i){
+    V_temp.row(i*4+0) = sData.V_o.row(sData.F(s[i],0));
+    V_temp.row(i*4+1) = sData.V_o.row(sData.F(s[i],1));
+    V_temp.row(i*4+2) = sData.V_o.row(sData.F(s[i],2));
+    V_temp.row(i*4+3) = sData.V_o.row(sData.F(s[i],3));
+    F_temp.row(i*4+0) << (i*4)+0, (i*4)+1, (i*4)+3;
+    F_temp.row(i*4+1) << (i*4)+0, (i*4)+2, (i*4)+1;
+    F_temp.row(i*4+2) << (i*4)+3, (i*4)+2, (i*4)+0;
+    F_temp.row(i*4+3) << (i*4)+1, (i*4)+2, (i*4)+3;
+  }
+  viewer.data.set_mesh(V_temp,F_temp);
+  viewer.core.align_camera_center(V_temp,F_temp);
+  viewer.data.set_face_based(true);
+  viewer.core.show_lines = true;
+}
+
+int main(int argc, char *argv[]) {
+
+  cerr << "Press space for running an iteration." << std::endl;
+  cerr << "Syntax: " << argv[0] << " demo_number (1 to 3)" << std::endl;
+  cerr << "1. 2D unconstrained parametrization" << std::endl;
+  cerr << "2. 2D deformation with positional constraints" << std::endl;
+  cerr << "3. 3D mesh deformation with positional constraints" << std::endl;
+
+  demo_type = PARAM_2D;
+
+   if (argc == 2) {
+     switch (std::atoi(argv[1])) {
+       case 1: {
+         demo_type = PARAM_2D;
+         break;
+       } case 2: {
+         demo_type = SOFT_CONST;
+         break;
+       } case 3: {
+         demo_type = DEFORM_3D;
+         break;
+       }
+       default: {
+         cerr << "Wrong demo number - Please choose one between 1-3" << std:: endl;
+         exit(1);
+       }
+     }
+   }
+
+
+  // Launch the viewer
+  igl::viewer::Viewer viewer;
+  viewer.callback_key_down = &key_down;
+
+  // Disable wireframe
+  viewer.core.show_lines = false;
+
+  // Draw checkerboard texture
+  viewer.core.show_texture = false;
+
+  // First iteration
+  key_down(viewer, ' ', 0);
+
+  viewer.launch();
+
+  return 0;
+}
+
+void check_mesh_for_issues(Eigen::MatrixXd& V, Eigen::MatrixXi& F) {
+
+  Eigen::SparseMatrix<double> A;
+  igl::adjacency_matrix(F,A);
+
+  Eigen::MatrixXi C, Ci;
+  igl::components(A, C, Ci);
+
+  int connected_components = Ci.rows();
+  if (connected_components!=1) {
+    cout << "Error! Input has multiple connected components" << endl; exit(1);
+  }
+  int euler_char = igl::euler_characteristic(V, F);
+  if (!euler_char) {
+    cout << "Error! Input does not have a disk topology, it's euler char is " << euler_char << endl; exit(1);
+  }
+  bool is_edge_manifold = igl::is_edge_manifold(V, F);
+  if (!is_edge_manifold) {
+    cout << "Error! Input is not an edge manifold" << endl; exit(1);
+  }
+
+  Eigen::VectorXd areas; igl::doublearea(V,F,areas);
+  const double eps = 1e-14;
+  for (int i = 0; i < areas.rows(); i++) {
+    if (areas(i) < eps) {
+      cout << "Error! Input has zero area faces" << endl; exit(1);
+    }
+  }
+}
+
+void get_soft_constraint_for_circle(Eigen::MatrixXd& V_o, Eigen::MatrixXi& F, Eigen::VectorXi& b, Eigen::MatrixXd& bc) {
+
+    Eigen::VectorXi bnd;
+    igl::boundary_loop(F,bnd);
+    const int B_STEPS = 22; // constraint every B_STEPS vertices of the boundary
+
+    b.resize(bnd.rows()/B_STEPS);
+    bc.resize(b.rows(),2);
+
+    int c_idx = 0;
+    for (int i = B_STEPS; i < bnd.rows(); i+=B_STEPS) {
+        b(c_idx) = bnd(i);
+        c_idx++;
+    }
+
+    bc.row(0) = V_o.row(b(0)); // keep it there for now
+    bc.row(1) = V_o.row(b(2));
+    bc.row(2) = V_o.row(b(3));
+    bc.row(3) = V_o.row(b(4));
+    bc.row(4) = V_o.row(b(5));
+
+
+    bc.row(0) << V_o(b(0),0), 0;
+    bc.row(4) << V_o(b(4),0), 0;
+    bc.row(2) << V_o(b(2),0), 0.1;
+    bc.row(3) << V_o(b(3),0), 0.05;
+    bc.row(1) << V_o(b(1),0), -0.15;
+    bc.row(5) << V_o(b(5),0), +0.15;
+}
+
+void get_cube_corner_constraints(Eigen::MatrixXd& V, Eigen::MatrixXi& F, Eigen::VectorXi& b, Eigen::MatrixXd& bc) {
+  double min_x,max_x,min_y,max_y,min_z,max_z;
+  min_x = V.col(0).minCoeff(); max_x = V.col(0).maxCoeff();
+  min_y = V.col(1).minCoeff(); max_y = V.col(1).maxCoeff();
+  min_z = V.col(2).minCoeff(); max_z = V.col(2).maxCoeff();
+
+
+  // get all cube corners
+  b.resize(8,1); bc.resize(8,3);
+  int x;
+  for (int i = 0; i < V.rows(); i++) {
+    if (V.row(i) == Eigen::RowVector3d(min_x,min_y,min_z)) b(0) = i;
+    if (V.row(i) == Eigen::RowVector3d(min_x,min_y,max_z)) b(1) = i;
+    if (V.row(i) == Eigen::RowVector3d(min_x,max_y,min_z)) b(2) = i;
+    if (V.row(i) == Eigen::RowVector3d(min_x,max_y,max_z)) b(3) = i;
+    if (V.row(i) == Eigen::RowVector3d(max_x,min_y,min_z)) b(4) = i;
+    if (V.row(i) == Eigen::RowVector3d(max_x,max_y,min_z)) b(5) = i;
+    if (V.row(i) == Eigen::RowVector3d(max_x,min_y,max_z)) b(6) = i;
+    if (V.row(i) == Eigen::RowVector3d(max_x,max_y,max_z)) b(7) = i;
+  }
+
+  // get all cube edges
+  std::set<int> cube_edge1; Eigen::VectorXi cube_edge1_vec;
+  for (int i = 0; i < V.rows(); i++) {
+    if ((V(i,0) == min_x && V(i,1) == min_y)) {
+      cube_edge1.insert(i);
+    }
+  }
+  Eigen::VectorXi edge1;
+  int_set_to_eigen_vector(cube_edge1, edge1);
+
+  std::set<int> cube_edge2; Eigen::VectorXi edge2;
+  for (int i = 0; i < V.rows(); i++) {
+    if ((V(i,0) == max_x && V(i,1) == max_y)) {
+      cube_edge2.insert(i);
+    }
+  }
+  int_set_to_eigen_vector(cube_edge2, edge2);
+  b = igl::cat(1,edge1,edge2);
+
+  std::set<int> cube_edge3; Eigen::VectorXi edge3;
+  for (int i = 0; i < V.rows(); i++) {
+    if ((V(i,0) == max_x && V(i,1) == min_y)) {
+      cube_edge3.insert(i);
+    }
+  }
+  int_set_to_eigen_vector(cube_edge3, edge3);
+  b = igl::cat(1,b,edge3);
+
+  std::set<int> cube_edge4; Eigen::VectorXi edge4;
+  for (int i = 0; i < V.rows(); i++) {
+    if ((V(i,0) == min_x && V(i,1) == max_y)) {
+      cube_edge4.insert(i);
+    }
+  }
+  int_set_to_eigen_vector(cube_edge4, edge4);
+  b = igl::cat(1,b,edge4);
+
+  bc.resize(b.rows(),3);
+  Eigen::Matrix3d m; m = Eigen::AngleAxisd(0.3 * M_PI, Eigen::Vector3d(1./sqrt(2.),1./sqrt(2.),0.)/*Eigen::Vector3d::UnitX()*/);
+  int i = 0;
+  for (; i < cube_edge1.size(); i++) {
+    Eigen::RowVector3d edge_rot_center(min_x,min_y,(min_z+max_z)/2.);
+    bc.row(i) = (V.row(b(i)) - edge_rot_center) * m + edge_rot_center;
+  }
+  for (; i < cube_edge1.size() + cube_edge2.size(); i++) {
+    Eigen::RowVector3d edge_rot_center(max_x,max_y,(min_z+max_z)/2.);
+    bc.row(i) = (V.row(b(i)) - edge_rot_center) * m.transpose() + edge_rot_center;
+  }
+  for (; i < cube_edge1.size() + cube_edge2.size() + cube_edge3.size(); i++) {
+    bc.row(i) = 0.75*V.row(b(i));
+  }
+  for (; i < b.rows(); i++) {
+    bc.row(i) = 0.75*V.row(b(i));
+  }
+}
+
+void int_set_to_eigen_vector(const std::set<int>& int_set, Eigen::VectorXi& vec) {
+  vec.resize(int_set.size()); int idx = 0;
+  for(auto f : int_set) {
+      vec(idx) = f; idx++;
+    }
+}

tutorial/CMakeLists.txt

@@ -167,6 +167,7 @@ if(TUTORIALS_CHAPTER6)
     add_subdirectory("609_Boolean")
     add_subdirectory("610_CSGTree")
   endif()
+  add_subdirectory("611_SLIM")
 endif()
 
 # Chapter 7

tutorial/images/611_SLIM.png.REMOVED.git-id

@@ -0,0 +1 @@
+9fcca7386e32522dc9425a7722ab7d7adb8c3f56

tutorial/shared/circle.obj.REMOVED.git-id

@@ -0,0 +1 @@
+8765e0b494718013bdd1305490a1d60d637f1b61

tutorial/shared/cube_40k.obj.REMOVED.git-id

@@ -0,0 +1 @@
+15badea46771f5da1ce678b5321d01b84c1ac6da

tutorial/shared/face.obj.REMOVED.git-id

@@ -0,0 +1 @@
+b07bc4c5241330e8bc15cef5d50c166e4346240c

tutorial/tutorial.html.REMOVED.git-id

@@ -1 +1 @@
-eaaebe95f03545b12c9d7bbf6689ca02aa70879f
+02cf5efb916ce8e885dcaeee5a89f91fa4e1f440

tutorial/tutorial.md.REMOVED.git-id

@@ -1 +1 @@
-8669ff1b3270a82bb3266fefdb7c381ed42e728d
+521b2da65bd08996d25fb1fab30ff66f704da312