* added missing IGL_INLINE

* added namespaces in flip_avoiding_line_search


Former-commit-id: 42da81abbced3e90d385407345b7fe35b9c12a3e

include/igl/flip_avoiding_line_search.cpp

@@ -11,280 +11,290 @@
 #include <Eigen/Dense>
 #include <vector>
 
-//---------------------------------------------------------------------------
-// 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
-  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);
+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);
+        }
     }
-    else
+
+    IGL_INLINE double get_smallest_pos_quad_zero(double a,double b, double c)
     {
-      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)
+      using namespace std;
+      double t1,t2;
+      if (a != 0)
       {
-        x[2]=x[1]; return(2);
+        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);
       }
-      return(1);
-    }
-}
+      else
+      {
+        t1 = t2 = -b/c;
+      }
+      assert (std::isfinite(t1));
+      assert (std::isfinite(t2));
 
-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 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 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
+    IGL_INLINE double get_min_pos_root_2D(const Eigen::MatrixXd& uv,
+                                          const Eigen::MatrixXi& F,
+                                          Eigen::MatrixXd& d,
+                                          int f)
     {
-      return t1;
+      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);
     }
-  }
-  else
-  {
-    return INFINITY;
-  }
-}
-
- 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);
-}
 
-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:
+    IGL_INLINE double get_min_pos_root_3D(const Eigen::MatrixXd& uv,
+                                          const Eigen::MatrixXi& F,
+                                          Eigen::MatrixXd& direc,
+                                          int f)
     {
-      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;
+      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;
+        }
+      }
     }
-  }
-}
 
-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++)
+    IGL_INLINE double compute_max_step_from_singularities(const Eigen::MatrixXd& uv,
+                                                          const Eigen::MatrixXi& F,
+                                                          Eigen::MatrixXd& d)
     {
-      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);
+      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);
+        }
+      }
+      return max_step;
     }
   }
-  return max_step;
 }
 
 IGL_INLINE double igl::flip_avoiding_line_search(
@@ -297,7 +307,7 @@ IGL_INLINE double igl::flip_avoiding_line_search(
   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 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);

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/slim.cpp

@@ -39,40 +39,38 @@ namespace igl
   namespace slim
   {
     // Definitions of internal functions
-    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);
-    void buildA(igl::SLIMData& s, Eigen::SparseMatrix<double> &A);
-    void buildRhs(igl::SLIMData& s, const Eigen::SparseMatrix<double> &At);
-    void add_soft_constraints(igl::SLIMData& s, Eigen::SparseMatrix<double> &L);
-    double compute_energy(igl::SLIMData& s, Eigen::MatrixXd &V_new);
-    double compute_soft_const_energy(igl::SLIMData& s,
-                                     const Eigen::MatrixXd &V,
-                                     const Eigen::MatrixXi &F,
-                                     Eigen::MatrixXd &V_o);
-    double compute_energy_with_jacobians(igl::SLIMData& s,
-                                         const Eigen::MatrixXd &V,
-                                         const Eigen::MatrixXi &F, const Eigen::MatrixXd &Ji,
-                                         Eigen::MatrixXd &uv, Eigen::VectorXd &areas);
-
-    void solve_weighted_arap(igl::SLIMData& s,
-                             const Eigen::MatrixXd &V,
-                             const Eigen::MatrixXi &F,
-                             Eigen::MatrixXd &uv,
-                             Eigen::VectorXi &soft_b_p,
-                             Eigen::MatrixXd &soft_bc_p);
-
-    void update_weights_and_closest_rotations(igl::SLIMData& s,
-                                              const Eigen::MatrixXd &V,
-                                              const Eigen::MatrixXi &F,
-                                              Eigen::MatrixXd &uv);
-    void compute_jacobians(igl::SLIMData& s, const Eigen::MatrixXd &uv);
-    void build_linear_system(igl::SLIMData& s, Eigen::SparseMatrix<double> &L);
-    void pre_calc(igl::SLIMData& s);
+    IGL_INLINE 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);
+    IGL_INLINE void buildA(igl::SLIMData& s, Eigen::SparseMatrix<double> &A);
+    IGL_INLINE void buildRhs(igl::SLIMData& s, const Eigen::SparseMatrix<double> &At);
+    IGL_INLINE void add_soft_constraints(igl::SLIMData& s, Eigen::SparseMatrix<double> &L);
+    IGL_INLINE double compute_energy(igl::SLIMData& s, Eigen::MatrixXd &V_new);
+    IGL_INLINE double compute_soft_const_energy(igl::SLIMData& s,
+                                                const Eigen::MatrixXd &V,
+                                                const Eigen::MatrixXi &F,
+                                                Eigen::MatrixXd &V_o);
+    IGL_INLINE double compute_energy_with_jacobians(igl::SLIMData& s,
+                                                    const Eigen::MatrixXd &V,
+                                                    const Eigen::MatrixXi &F, const Eigen::MatrixXd &Ji,
+                                                    Eigen::MatrixXd &uv, Eigen::VectorXd &areas);
+    IGL_INLINE void solve_weighted_arap(igl::SLIMData& s,
+                                        const Eigen::MatrixXd &V,
+                                        const Eigen::MatrixXi &F,
+                                        Eigen::MatrixXd &uv,
+                                        Eigen::VectorXi &soft_b_p,
+                                        Eigen::MatrixXd &soft_bc_p);
+    IGL_INLINE void update_weights_and_closest_rotations( igl::SLIMData& s,
+                                                          const Eigen::MatrixXd &V,
+                                                          const Eigen::MatrixXi &F,
+                                                          Eigen::MatrixXd &uv);
+    IGL_INLINE void compute_jacobians(igl::SLIMData& s, const Eigen::MatrixXd &uv);
+    IGL_INLINE void build_linear_system(igl::SLIMData& s, Eigen::SparseMatrix<double> &L);
+    IGL_INLINE void pre_calc(igl::SLIMData& s);
 
     // Implementation
-    void compute_surface_gradient_matrix(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F,
-                                         const Eigen::MatrixXd &F1, const Eigen::MatrixXd &F2,
+    IGL_INLINE 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)
     {
 
@@ -86,7 +84,7 @@ namespace igl
       D2 = F2.col(0).asDiagonal() * Dx + F2.col(1).asDiagonal() * Dy + F2.col(2).asDiagonal() * Dz;
     }
 
-    void compute_jacobians(igl::SLIMData& s, const Eigen::MatrixXd &uv)
+    IGL_INLINE void compute_jacobians(igl::SLIMData& s, const Eigen::MatrixXd &uv)
     {
       if (s.F.cols() == 3)
       {
@@ -110,7 +108,7 @@ namespace igl
       }
     }
 
-    void update_weights_and_closest_rotations(igl::SLIMData& s,
+    IGL_INLINE void update_weights_and_closest_rotations(igl::SLIMData& s,
                                               const Eigen::MatrixXd &V,
                                               const Eigen::MatrixXi &F,
                                               Eigen::MatrixXd &uv)
@@ -385,12 +383,12 @@ namespace igl
 
     }
 
-    void solve_weighted_arap(igl::SLIMData& s,
-                             const Eigen::MatrixXd &V,
-                             const Eigen::MatrixXi &F,
-                             Eigen::MatrixXd &uv,
-                             Eigen::VectorXi &soft_b_p,
-                             Eigen::MatrixXd &soft_bc_p)
+    IGL_INLINE void solve_weighted_arap(igl::SLIMData& s,
+                                        const Eigen::MatrixXd &V,
+                                        const Eigen::MatrixXi &F,
+                                        Eigen::MatrixXd &uv,
+                                        Eigen::VectorXi &soft_b_p,
+                                        Eigen::MatrixXd &soft_bc_p)
     {
       using namespace Eigen;
 
@@ -418,7 +416,7 @@ namespace igl
     }
 
 
-    void pre_calc(igl::SLIMData& s)
+    IGL_INLINE void pre_calc(igl::SLIMData& s)
     {
       if (!s.has_pre_calc)
       {
@@ -477,7 +475,7 @@ namespace igl
       }
     }
 
-    void build_linear_system(igl::SLIMData& s, Eigen::SparseMatrix<double> &L)
+    IGL_INLINE void build_linear_system(igl::SLIMData& s, Eigen::SparseMatrix<double> &L)
     {
       // formula (35) in paper
       Eigen::SparseMatrix<double> A(s.dim * s.dim * s.f_n, s.dim * s.v_n);
@@ -499,7 +497,7 @@ namespace igl
       L.makeCompressed();
     }
 
-    void add_soft_constraints(igl::SLIMData& s, Eigen::SparseMatrix<double> &L)
+    IGL_INLINE void add_soft_constraints(igl::SLIMData& s, Eigen::SparseMatrix<double> &L)
     {
       int v_n = s.v_num;
       for (int d = 0; d < s.dim; d++)
@@ -513,17 +511,17 @@ namespace igl
       }
     }
 
-    double compute_energy(igl::SLIMData& s, Eigen::MatrixXd &V_new)
+    IGL_INLINE double compute_energy(igl::SLIMData& s, Eigen::MatrixXd &V_new)
     {
       compute_jacobians(s,V_new);
       return compute_energy_with_jacobians(s, s.V, s.F, s.Ji, V_new, s.M) +
              compute_soft_const_energy(s, s.V, s.F, V_new);
     }
 
-    double compute_soft_const_energy(igl::SLIMData& s,
-                                     const Eigen::MatrixXd &V,
-                                     const Eigen::MatrixXi &F,
-                                     Eigen::MatrixXd &V_o)
+    IGL_INLINE double compute_soft_const_energy(igl::SLIMData& s,
+                                                const Eigen::MatrixXd &V,
+                                                const Eigen::MatrixXi &F,
+                                                Eigen::MatrixXd &V_o)
     {
       double e = 0;
       for (int i = 0; i < s.b.rows(); i++)
@@ -533,10 +531,10 @@ namespace igl
       return e;
     }
 
-    double compute_energy_with_jacobians(igl::SLIMData& s,
-                                         const Eigen::MatrixXd &V,
-                                         const Eigen::MatrixXi &F, const Eigen::MatrixXd &Ji,
-                                         Eigen::MatrixXd &uv, Eigen::VectorXd &areas)
+    IGL_INLINE double compute_energy_with_jacobians(igl::SLIMData& s,
+                                                    const Eigen::MatrixXd &V,
+                                                    const Eigen::MatrixXi &F, const Eigen::MatrixXd &Ji,
+                                                    Eigen::MatrixXd &uv, Eigen::VectorXd &areas)
     {
 
       double energy = 0;
@@ -659,7 +657,7 @@ namespace igl
       return energy;
     }
 
-    void buildA(igl::SLIMData& s, Eigen::SparseMatrix<double> &A)
+    IGL_INLINE void buildA(igl::SLIMData& s, Eigen::SparseMatrix<double> &A)
     {
       // formula (35) in paper
       std::vector<Eigen::Triplet<double> > IJV;
@@ -785,7 +783,7 @@ namespace igl
       A.setFromTriplets(IJV.begin(), IJV.end());
     }
 
-    void buildRhs(igl::SLIMData& s, const Eigen::SparseMatrix<double> &At)
+    IGL_INLINE void buildRhs(igl::SLIMData& s, const Eigen::SparseMatrix<double> &At)
     {
       Eigen::VectorXd f_rhs(s.dim * s.dim * s.f_n);
       f_rhs.setZero();