merge

cmake/recipes/external/catch2.cmake

@@ -8,7 +8,7 @@ include(FetchContent)
 FetchContent_Declare(
     catch2
     GIT_REPOSITORY https://github.com/catchorg/Catch2.git
-    GIT_TAG v2.13.3
+    GIT_TAG v2.13.8
     GIT_SHALLOW TRUE
 )
 FetchContent_MakeAvailable(catch2)

cmake/recipes/external/embree.cmake

@@ -8,7 +8,7 @@ include(FetchContent)
 FetchContent_Declare(
     embree
     GIT_REPOSITORY https://github.com/embree/embree.git
-    GIT_TAG        v3.12.1
+    GIT_TAG        v3.13.3
     GIT_SHALLOW    TRUE
 )
 
@@ -45,7 +45,6 @@ add_library(embree::embree ALIAS embree)
 
 # Some order for IDEs
 set_target_properties(embree PROPERTIES FOLDER "ThirdParty//Embree")
-set_target_properties(algorithms PROPERTIES FOLDER "ThirdParty//Embree")
 set_target_properties(lexers PROPERTIES FOLDER "ThirdParty//Embree")
 set_target_properties(math PROPERTIES FOLDER "ThirdParty//Embree")
 set_target_properties(simd PROPERTIES FOLDER "ThirdParty//Embree")

include/igl/FastWindingNumberForSoups.h

@@ -1,4 +1,4 @@
-// This header created by issuing: `echo "// This header created by issuing: \`$BASH_COMMAND\` $(echo "" | cat - LICENSE README.md | sed -e "s#^..*#\/\/ &#") $(echo "" | cat - SYS_Types.h SYS_Math.h VM_SSEFunc.h VM_SIMD.h UT_Array.h UT_ArrayImpl.h UT_SmallArray.h UT_FixedVector.h UT_ParallelUtil.h UT_BVH.h UT_BVHImpl.h UT_SolidAngle.h UT_Array.cpp UT_SolidAngle.cpp | sed -e "s/^#.*include  *\".*$//g")" > ../FastWindingNumberForSoups.h` 
+// This header created by issuing: `echo "// This header created by issuing: \`$BASH_COMMAND\` $(echo "" | cat - LICENSE README.md | sed -e "s#^..*#\/\/ &#") $(echo "" | cat - SYS_Types.h SYS_Math.h VM_SSEFunc.h VM_SIMDFunc.h VM_SIMD.h UT_Array.h UT_ArrayImpl.h UT_SmallArray.h UT_FixedVector.h UT_ParallelUtil.h UT_BVH.h UT_BVHImpl.h UT_SolidAngle.h UT_Array.cpp UT_SolidAngle.cpp | sed -e "s/^#.*include  *\".*$//g")" > ~/Repos/libigl/include/igl/FastWindingNumberForSoups.h` 
 // MIT License
 
 // Copyright (c) 2018 Side Effects Software Inc.
@@ -48,7 +48,7 @@
 
 // Create a single header using:
 
-//     echo "// This header created by issuing: \`$BASH_COMMAND\` $(echo "" | cat - LICENSE README.md | sed -e "s#^..*#\/\/&#") $(echo "" | cat - SYS_Types.h SYS_Math.h VM_SSEFunc.h VM_SIMD.h UT_Array.h UT_ArrayImpl.h UT_SmallArray.h UT_FixedVector.h UT_ParallelUtil.h UT_BVH.h UT_BVHImpl.h UT_SolidAngle.h UT_Array.cpp UT_SolidAngle.cpp | sed -e "s/^#.*include  *\".*$//g")" 
+//     echo "// This header created by issuing: \`$BASH_COMMAND\` $(echo "" | cat - LICENSE README.md | sed -e "s#^..*#\/\/ &#") $(echo "" | cat - SYS_Types.h SYS_Math.h VM_SSEFunc.h VM_SIMD.h UT_Array.h UT_ArrayImpl.h UT_SmallArray.h UT_FixedVector.h UT_ParallelUtil.h UT_BVH.h UT_BVHImpl.h UT_SolidAngle.h UT_Array.cpp UT_SolidAngle.cpp | sed -e "s/^#.*include  *\".*$//g")" 
 /*
  * Copyright (c) 2018 Side Effects Software Inc.
  *
@@ -357,6 +357,7 @@ static inline fpreal64 SYSabs(fpreal64 a) { return ::fabs(a); }
  */
 
 #pragma once
+#ifdef __SSE__
 
 #ifndef __VM_SSEFunc__
 #define __VM_SSEFunc__
@@ -734,6 +735,407 @@ vm_allbits(const v4si &a)
 }}
 
 #endif
+#endif
+#pragma once
+#ifndef __SSE__
+#ifndef __VM_SIMDFunc__
+#define __VM_SIMDFunc__
+
+
+
+#include <cmath>
+
+namespace igl { namespace FastWindingNumber {
+
+struct v4si {
+	int32 v[4];
+};
+
+struct v4sf {
+	float v[4];
+};
+
+static SYS_FORCE_INLINE v4sf V4SF(const v4si &v) {
+	static_assert(sizeof(v4si) == sizeof(v4sf) && alignof(v4si) == alignof(v4sf), "v4si and v4sf must be compatible");
+	return *(const v4sf*)&v;
+}
+
+static SYS_FORCE_INLINE v4si V4SI(const v4sf &v) {
+	static_assert(sizeof(v4si) == sizeof(v4sf) && alignof(v4si) == alignof(v4sf), "v4si and v4sf must be compatible");
+	return *(const v4si*)&v;
+}
+
+static SYS_FORCE_INLINE int32 conditionMask(bool c) {
+	return c ? int32(0xFFFFFFFF) : 0;
+}
+
+static SYS_FORCE_INLINE v4sf
+VM_SPLATS(float f) {
+	return v4sf{{f, f, f, f}};
+}
+
+static SYS_FORCE_INLINE v4si
+VM_SPLATS(uint32 i) {
+	return v4si{{int32(i), int32(i), int32(i), int32(i)}};
+}
+
+static SYS_FORCE_INLINE v4si
+VM_SPLATS(int32 i) {
+	return v4si{{i, i, i, i}};
+}
+
+static SYS_FORCE_INLINE v4sf
+VM_SPLATS(float a, float b, float c, float d) {
+	return v4sf{{a, b, c, d}};
+}
+
+static SYS_FORCE_INLINE v4si
+VM_SPLATS(uint32 a, uint32 b, uint32 c, uint32 d) {
+	return v4si{{int32(a), int32(b), int32(c), int32(d)}};
+}
+
+static SYS_FORCE_INLINE v4si
+VM_SPLATS(int32 a, int32 b, int32 c, int32 d) {
+	return v4si{{a, b, c, d}};
+}
+
+static SYS_FORCE_INLINE v4si
+VM_LOAD(const int32 v[4]) {
+	return v4si{{v[0], v[1], v[2], v[3]}};
+}
+
+static SYS_FORCE_INLINE v4sf
+VM_LOAD(const float v[4]) {
+	return v4sf{{v[0], v[1], v[2], v[3]}};
+}
+
+
+static inline v4si VM_ICMPEQ(v4si a, v4si b) {
+	return v4si{{
+		conditionMask(a.v[0] == b.v[0]),
+		conditionMask(a.v[1] == b.v[1]),
+		conditionMask(a.v[2] == b.v[2]),
+		conditionMask(a.v[3] == b.v[3])
+	}};
+}
+
+static inline v4si VM_ICMPGT(v4si a, v4si b) {
+	return v4si{{
+		conditionMask(a.v[0] > b.v[0]),
+		conditionMask(a.v[1] > b.v[1]),
+		conditionMask(a.v[2] > b.v[2]),
+		conditionMask(a.v[3] > b.v[3])
+	}};
+}
+
+static inline v4si VM_ICMPLT(v4si a, v4si b) {
+	return v4si{{
+		conditionMask(a.v[0] < b.v[0]),
+		conditionMask(a.v[1] < b.v[1]),
+		conditionMask(a.v[2] < b.v[2]),
+		conditionMask(a.v[3] < b.v[3])
+	}};
+}
+
+static inline v4si VM_IADD(v4si a, v4si b) {
+	return v4si{{
+		(a.v[0] + b.v[0]),
+		(a.v[1] + b.v[1]),
+		(a.v[2] + b.v[2]),
+		(a.v[3] + b.v[3])
+	}};
+}
+
+static inline v4si VM_ISUB(v4si a, v4si b) {
+	return v4si{{
+		(a.v[0] - b.v[0]),
+		(a.v[1] - b.v[1]),
+		(a.v[2] - b.v[2]),
+		(a.v[3] - b.v[3])
+	}};
+}
+
+static inline v4si VM_OR(v4si a, v4si b) {
+	return v4si{{
+		(a.v[0] | b.v[0]),
+		(a.v[1] | b.v[1]),
+		(a.v[2] | b.v[2]),
+		(a.v[3] | b.v[3])
+	}};
+}
+
+static inline v4si VM_AND(v4si a, v4si b) {
+	return v4si{{
+		(a.v[0] & b.v[0]),
+		(a.v[1] & b.v[1]),
+		(a.v[2] & b.v[2]),
+		(a.v[3] & b.v[3])
+	}};
+}
+
+static inline v4si VM_ANDNOT(v4si a, v4si b) {
+	return v4si{{
+		((~a.v[0]) & b.v[0]),
+		((~a.v[1]) & b.v[1]),
+		((~a.v[2]) & b.v[2]),
+		((~a.v[3]) & b.v[3])
+	}};
+}
+
+static inline v4si VM_XOR(v4si a, v4si b) {
+	return v4si{{
+		(a.v[0] ^ b.v[0]),
+		(a.v[1] ^ b.v[1]),
+		(a.v[2] ^ b.v[2]),
+		(a.v[3] ^ b.v[3])
+	}};
+}
+
+static SYS_FORCE_INLINE int
+VM_EXTRACT(const v4si v, int index) {
+	return v.v[index];
+}
+
+static SYS_FORCE_INLINE float
+VM_EXTRACT(const v4sf v, int index) {
+	return v.v[index];
+}
+
+static SYS_FORCE_INLINE v4si
+VM_INSERT(v4si v, int32 value, int index) {
+	v.v[index] = value;
+	return v;
+}
+
+static SYS_FORCE_INLINE v4sf
+VM_INSERT(v4sf v, float value, int index) {
+	v.v[index] = value;
+	return v;
+}
+
+static inline v4si VM_CMPEQ(v4sf a, v4sf b) {
+	return v4si{{
+		conditionMask(a.v[0] == b.v[0]),
+		conditionMask(a.v[1] == b.v[1]),
+		conditionMask(a.v[2] == b.v[2]),
+		conditionMask(a.v[3] == b.v[3])
+	}};
+}
+
+static inline v4si VM_CMPNE(v4sf a, v4sf b) {
+	return v4si{{
+		conditionMask(a.v[0] != b.v[0]),
+		conditionMask(a.v[1] != b.v[1]),
+		conditionMask(a.v[2] != b.v[2]),
+		conditionMask(a.v[3] != b.v[3])
+	}};
+}
+
+static inline v4si VM_CMPGT(v4sf a, v4sf b) {
+	return v4si{{
+		conditionMask(a.v[0] > b.v[0]),
+		conditionMask(a.v[1] > b.v[1]),
+		conditionMask(a.v[2] > b.v[2]),
+		conditionMask(a.v[3] > b.v[3])
+	}};
+}
+
+static inline v4si VM_CMPLT(v4sf a, v4sf b) {
+	return v4si{{
+		conditionMask(a.v[0] < b.v[0]),
+		conditionMask(a.v[1] < b.v[1]),
+		conditionMask(a.v[2] < b.v[2]),
+		conditionMask(a.v[3] < b.v[3])
+	}};
+}
+
+static inline v4si VM_CMPGE(v4sf a, v4sf b) {
+	return v4si{{
+		conditionMask(a.v[0] >= b.v[0]),
+		conditionMask(a.v[1] >= b.v[1]),
+		conditionMask(a.v[2] >= b.v[2]),
+		conditionMask(a.v[3] >= b.v[3])
+	}};
+}
+
+static inline v4si VM_CMPLE(v4sf a, v4sf b) {
+	return v4si{{
+		conditionMask(a.v[0] <= b.v[0]),
+		conditionMask(a.v[1] <= b.v[1]),
+		conditionMask(a.v[2] <= b.v[2]),
+		conditionMask(a.v[3] <= b.v[3])
+	}};
+}
+
+static inline v4sf VM_ADD(v4sf a, v4sf b) {
+	return v4sf{{
+		(a.v[0] + b.v[0]),
+		(a.v[1] + b.v[1]),
+		(a.v[2] + b.v[2]),
+		(a.v[3] + b.v[3])
+	}};
+}
+
+static inline v4sf VM_SUB(v4sf a, v4sf b) {
+	return v4sf{{
+		(a.v[0] - b.v[0]),
+		(a.v[1] - b.v[1]),
+		(a.v[2] - b.v[2]),
+		(a.v[3] - b.v[3])
+	}};
+}
+
+static inline v4sf VM_NEG(v4sf a) {
+	return v4sf{{
+		(-a.v[0]),
+		(-a.v[1]),
+		(-a.v[2]),
+		(-a.v[3])
+	}};
+}
+
+static inline v4sf VM_MUL(v4sf a, v4sf b) {
+	return v4sf{{
+		(a.v[0] * b.v[0]),
+		(a.v[1] * b.v[1]),
+		(a.v[2] * b.v[2]),
+		(a.v[3] * b.v[3])
+	}};
+}
+
+static inline v4sf VM_DIV(v4sf a, v4sf b) {
+	return v4sf{{
+		(a.v[0] / b.v[0]),
+		(a.v[1] / b.v[1]),
+		(a.v[2] / b.v[2]),
+		(a.v[3] / b.v[3])
+	}};
+}
+
+static inline v4sf VM_MADD(v4sf a, v4sf b, v4sf c) {
+	return v4sf{{
+		(a.v[0] * b.v[0]) + c.v[0],
+		(a.v[1] * b.v[1]) + c.v[1],
+		(a.v[2] * b.v[2]) + c.v[2],
+		(a.v[3] * b.v[3]) + c.v[3]
+	}};
+}
+
+static inline v4sf VM_ABS(v4sf a) {
+	return v4sf{{
+		(a.v[0] < 0) ? -a.v[0] : a.v[0],
+		(a.v[1] < 0) ? -a.v[1] : a.v[1],
+		(a.v[2] < 0) ? -a.v[2] : a.v[2],
+		(a.v[3] < 0) ? -a.v[3] : a.v[3]
+	}};
+}
+
+static inline v4sf VM_MAX(v4sf a, v4sf b) {
+	return v4sf{{
+		(a.v[0] < b.v[0]) ? b.v[0] : a.v[0],
+		(a.v[1] < b.v[1]) ? b.v[1] : a.v[1],
+		(a.v[2] < b.v[2]) ? b.v[2] : a.v[2],
+		(a.v[3] < b.v[3]) ? b.v[3] : a.v[3]
+	}};
+}
+
+static inline v4sf VM_MIN(v4sf a, v4sf b) {
+	return v4sf{{
+		(a.v[0] > b.v[0]) ? b.v[0] : a.v[0],
+		(a.v[1] > b.v[1]) ? b.v[1] : a.v[1],
+		(a.v[2] > b.v[2]) ? b.v[2] : a.v[2],
+		(a.v[3] > b.v[3]) ? b.v[3] : a.v[3]
+	}};
+}
+
+static inline v4sf VM_INVERT(v4sf a) {
+	return v4sf{{
+		(1.0f/a.v[0]),
+		(1.0f/a.v[1]),
+		(1.0f/a.v[2]),
+		(1.0f/a.v[3])
+	}};
+}
+
+static inline v4sf VM_SQRT(v4sf a) {
+	return v4sf{{
+		std::sqrt(a.v[0]),
+		std::sqrt(a.v[1]),
+		std::sqrt(a.v[2]),
+		std::sqrt(a.v[3])
+	}};
+}
+
+static inline v4si VM_INT(v4sf a) {
+	return v4si{{
+		int32(a.v[0]),
+		int32(a.v[1]),
+		int32(a.v[2]),
+		int32(a.v[3])
+	}};
+}
+
+static inline v4sf VM_IFLOAT(v4si a) {
+	return v4sf{{
+		float(a.v[0]),
+		float(a.v[1]),
+		float(a.v[2]),
+		float(a.v[3])
+	}};
+}
+
+static SYS_FORCE_INLINE void VM_P_FLOOR() {}
+
+static SYS_FORCE_INLINE int32 singleIntFloor(float f) {
+	// Casting to int32 usually truncates toward zero, instead of rounding down,
+	// so subtract one if the result is above f.
+	int32 i = int32(f);
+	i -= (float(i) > f);
+	return i;
+}
+static inline v4si VM_FLOOR(v4sf a) {
+	return v4si{{
+		singleIntFloor(a.v[0]),
+		singleIntFloor(a.v[1]),
+		singleIntFloor(a.v[2]),
+		singleIntFloor(a.v[3])
+	}};
+}
+
+static SYS_FORCE_INLINE void VM_E_FLOOR() {}
+
+static SYS_FORCE_INLINE bool vm_allbits(v4si a) {
+	return (
+		(a.v[0] == -1) && 
+		(a.v[1] == -1) && 
+		(a.v[2] == -1) && 
+		(a.v[3] == -1)
+	);
+}
+
+int SYS_FORCE_INLINE _mm_movemask_ps(const v4si& v) {
+	return (
+		int(v.v[0] < 0) |
+		(int(v.v[1] < 0)<<1) |
+		(int(v.v[2] < 0)<<2) |
+		(int(v.v[3] < 0)<<3)
+	);
+}
+
+int SYS_FORCE_INLINE _mm_movemask_ps(const v4sf& v) {
+	// Use std::signbit just in case it needs to distinguish between +0 and -0
+	// or between positive and negative NaN values (e.g. these could really
+	// be integers instead of floats).
+	return (
+		int(std::signbit(v.v[0])) |
+		(int(std::signbit(v.v[1]))<<1) |
+		(int(std::signbit(v.v[2]))<<2) |
+		(int(std::signbit(v.v[3]))<<3)
+	);
+}
+}}
+#endif
+#endif
 /*
  * Copyright (c) 2018 Side Effects Software Inc.
  *
@@ -770,6 +1172,8 @@ vm_allbits(const v4si &a)
 //#define FORCE_NON_SIMD
 
 
+
+
 namespace igl { namespace FastWindingNumber {
 
 class v4uf;
@@ -986,11 +1390,13 @@ public:
         return base;
     }
 
+#ifdef __SSE__
     template <int A, int B, int C, int D>
     SYS_FORCE_INLINE v4uf swizzle() const
     { 
         return VM_SHUFFLE<A,B,C,D>(vector);
     }
+#endif
 
     SYS_FORCE_INLINE v4uu isFinite() const
     {

include/igl/IndexComparison.h

@@ -7,7 +7,6 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #ifndef IGL_INDEXCOMPARISON_H
 #define IGL_INDEXCOMPARISON_H
-#include <iostream>
 namespace igl{
   // Comparison struct used by sort
   // http://bytes.com/topic/c/answers/132045-sort-get-index

include/igl/bbw.cpp

@@ -118,12 +118,8 @@ IGL_INLINE bool igl::bbw(
     }
     W.col(i) = Wi;
   };
-#ifdef WIN32
-  for (int i = 0; i < m; ++i)
-    optimize_weight(i);
-#else
+
   parallel_for(m,optimize_weight,2);
-#endif
   if(error)
   {
     return false;

include/igl/doublearea.cpp

@@ -219,7 +219,7 @@ IGL_INLINE void igl::doublearea_quad(
   }
 
   // Compute areas
-  Eigen::VectorXd doublearea_tri;
+  Eigen::Matrix<typename DeriveddblA::Scalar, Eigen::Dynamic, 1> doublearea_tri;
   igl::doublearea(V,Ft,doublearea_tri);
 
   dblA.resize(F.rows(),1);

include/igl/fast_winding_number.cpp

@@ -1,6 +1,5 @@
 #include "fast_winding_number.h"
 #include "octree.h"
-#include "knn.h"
 #include "parallel_for.h"
 #include "PI.h"
 #include <vector>

include/igl/linprog.cpp

@@ -10,7 +10,6 @@
 #include "slice_into.h"
 #include "find.h"
 #include "colon.h"
-#include <iostream>
 
 //#define IGL_LINPROG_VERBOSE
 IGL_INLINE bool igl::linprog(
@@ -242,7 +241,7 @@ IGL_INLINE bool igl::linprog(
     s*=1e6*c.array().abs().maxCoeff();
     s.head(n) = c;
   }
-  x.resize(std::max(B.maxCoeff()+1,n));
+  x.setZero(std::max(B.maxCoeff()+1,n));
   igl::slice_into(xb,B,x);
   x = x.head(n).eval();
   return success;
@@ -281,9 +280,12 @@ IGL_INLINE bool igl::linprog(
   MatrixXd BSP(0,2*n+m);
   if(p>0)
   {
-    MatrixXd BS(p,2*n);
-    BS<<B,MatrixXd::Zero(p,n);
+    // B ∈ ℝ^(p × n)
+    MatrixXd BS(p,n+m);
+    BS<<B,MatrixXd::Zero(p,m);
+    // BS ∈ ℝ^(p × n+m)
     BSP = BS*P;
+    // BSP ∈ ℝ^(p × 2n+m)
   }
 
   VectorXd fSP = VectorXd::Ones(2*n+m);
@@ -296,7 +298,24 @@ IGL_INLINE bool igl::linprog(
   bb<<bS,c;
 
   VectorXd xxs;
+  // min   ccᵀxxs
+  // s.t.  AA xxs ≤ bb
+  //          xxs ≥ 0
+  //        
+  // x = x⁺ - x⁻
+  //
+  //    P
+  // .--^---.
+  // [I -I 0  [x⁺   = [x
+  //  0  0 I]  x⁻      s]
+  //           s]
+  // Pᵀ [xᵀ sᵀ] = xxsᵀ
+  //
+  // min  [fᵀ -fᵀ 𝟙ᵀ] [x⁺;x⁻;s]
+  // s.t.  AA [x⁺;x⁻;s] ≤ b
+  // s.t.  [x⁺;x⁻;s] ≥ 0
   bool ret = linprog(cc,AA,bb,0,xxs);
+  // x = P(1:n,:) xxs
   x = P.block(0,0,n,2*n+m)*xxs;
   return ret;
 }

include/igl/linprog.h

@@ -13,7 +13,7 @@ namespace igl
 {
   // Solve a linear program given in "standard form"
   //
-  // min  f'x
+  // min  c'x
   // s.t. A(    1:k,:) x <= b(1:k)
   //      A(k+1:end,:) x = b(k+1:end)
   //   ** x >= 0 **
@@ -33,7 +33,7 @@ namespace igl
     const Eigen::MatrixXd & A,
     const Eigen::VectorXd & b,
     const int k,
-    Eigen::VectorXd & f);
+    Eigen::VectorXd & x);
   
   // Wrapper in friendlier general form (no implicit bounds on x)
   //

include/igl/mapping_energy_with_jacobians.cpp

@@ -125,7 +125,7 @@ IGL_INLINE double igl::mapping_energy_with_jacobians(
         }
         case igl::MappingEnergyType::EXP_CONFORMAL:
         {
-          energy += areas(i) * exp((pow(s1, 2) + pow(s2, 2) + pow(s3, 2)) / (3 * pow(s1 * s2 * s3, 2. / 3.)));
+          energy += areas(i) * exp(exp_factor * (pow(s1, 2) + pow(s2, 2) + pow(s3, 2)) / (3 * pow(s1 * s2 * s3, 2. / 3.)));
           break;
         }
       }

include/igl/marching_cubes.cpp

@@ -52,8 +52,8 @@ IGL_INLINE void igl::marching_cubes(
     const unsigned i = xyz2i(x,y,z);
 
     //Make a local copy of the values at the cube's corners
-    static Eigen::Matrix<Scalar,8,1> cS;
-    static Eigen::Matrix<Index,8,1> cI;
+    Eigen::Matrix<Scalar,8,1> cS;
+    Eigen::Matrix<Index,8,1> cI;
     //Find which vertices are inside of the surface and which are outside
     for(int c = 0; c < 8; c++)
     {
@@ -62,7 +62,7 @@ IGL_INLINE void igl::marching_cubes(
       cS(c) = S(ic);
     }
 
-  march_cube(GV,cS,cI,isovalue,V,n,F,m,E2V);
+    march_cube(GV,cS,cI,isovalue,V,n,F,m,E2V);
 
   };
 
@@ -110,11 +110,12 @@ IGL_INLINE void igl::marching_cubes(
   Index m = 0;
 
   // march over cubes
+
+  //Make a local copy of the values at the cube's corners
+  Eigen::Matrix<Scalar, 8, 1> cS;
+  Eigen::Matrix<Index, 8, 1> cI;
   for(Index c = 0;c<GI.rows();c++)
   {
-    //Make a local copy of the values at the cube's corners
-    static Eigen::Matrix<Scalar,8,1> cS;
-    static Eigen::Matrix<Index,8,1> cI;
     for(int v = 0; v < 8; v++)
     {
       cI(v) = GI(c,v);

include/igl/octree.cpp

@@ -1,6 +1,5 @@
 #include "octree.h"
 #include <vector>
-#include <queue>
 
 namespace igl {
   template <typename DerivedP, typename IndexType, typename DerivedCH,

include/igl/opengl/ViewerCore.cpp

@@ -363,9 +363,9 @@ IGL_INLINE void igl::opengl::ViewerCore::draw_labels(
   float width  = viewport(2);
   float height = viewport(3);
   float text_shift_scale_factor = orthographic ? 0.01 : 0.03;
-  float render_scale = orthographic ? 0.6 : 1.7;
+  float render_scale = (orthographic ? 0.6 : 1.7) * data.label_size;
   glUniform1f(glGetUniformLocation(data.meshgl.shader_text, "TextShiftFactor"), text_shift_scale_factor);
-  glUniform3f(glGetUniformLocation(data.meshgl.shader_text, "TextColor"), 0, 0, 0);
+  glUniform3f(glGetUniformLocation(data.meshgl.shader_text, "TextColor"), data.label_color(0), data.label_color(1), data.label_color(2));
   glUniform2f(glGetUniformLocation(data.meshgl.shader_text, "CellSize"), 1.0f / 16, (300.0f / 384) / 6);
   glUniform2f(glGetUniformLocation(data.meshgl.shader_text, "CellOffset"), 0.5 / 256.0, 0.5 / 256.0);
   glUniform2f(glGetUniformLocation(data.meshgl.shader_text, "RenderSize"), 

include/igl/opengl/ViewerData.cpp

@@ -35,6 +35,7 @@ IGL_INLINE igl::opengl::ViewerData::ViewerData()
   use_matcap        (false),
   point_size(30),
   line_width(0.5f),
+  label_size(1),
   line_color(0,0,0,1),
   label_color(0,0,0.04,1),
   shininess(35.0f),

include/igl/opengl/ViewerData.h

@@ -276,6 +276,7 @@ public:
   float point_size;
   // line_width is NOT SUPPORTED on Mac OS and Windows
   float line_width;
+  float label_size;
   Eigen::Matrix<float, 4, 1, Eigen::DontAlign> line_color;
   Eigen::Matrix<float, 4, 1, Eigen::DontAlign> label_color;
 

include/igl/readPLY.cpp

@@ -171,6 +171,7 @@ IGL_INLINE bool readPLY(
   try
   {
     std::vector<uint8_t> fileBufferBytes;
+    // read_file_binary will call fclose
     read_file_binary(fp,fileBufferBytes);
     FileMemoryStream stream((char*)fileBufferBytes.data(), fileBufferBytes.size());
     return readPLY(stream,V,F,E,N,UV,VD,Vheader,FD,Fheader,ED,Eheader,comments);
@@ -179,6 +180,7 @@ IGL_INLINE bool readPLY(
   {
     std::cerr << "ReadPLY error: " << e.what() << std::endl;
   }
+  fclose(fp);
   return false;
 }
 
@@ -448,7 +450,7 @@ IGL_INLINE bool readPLY(
   else
   {
     FD.resize(faces->count, _face_header.size());
-    tinyply_buffer_to_matrix(*_face_data, FD, faces->count, 1);
+    tinyply_buffer_to_matrix(*_face_data, FD, faces->count, _face_header.size());
   }
 
   /// convert edge data:

include/igl/slice.cpp

@@ -188,6 +188,8 @@ IGL_INLINE DerivedX igl::slice(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+template Eigen::Matrix<double, -1, 1, 0, -1, 1> igl::slice<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::DenseBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&);
+template Eigen::Matrix<int, -1, 1, 0, -1, 1> igl::slice<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&);
 // generated by autoexplicit.sh
 template void igl::slice<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Block<Eigen::Matrix<int, -1, -1, 0, -1, -1>, -1, 1, true>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::DenseBase<Eigen::Block<Eigen::Matrix<int, -1, -1, 0, -1, -1>, -1, 1, true> > const&, int, Eigen::Matrix<int, -1, 1, 0, -1, 1>&);
 // generated by autoexplicit.sh

include/igl/triangle/scaf.h

@@ -46,7 +46,7 @@ namespace igl
       Eigen::VectorXd m_M; // mesh area or volume
       Eigen::VectorXd s_M; // scaffold area or volume
       Eigen::VectorXd w_M; // area/volume weights for whole
-      double mesh_measure; // area or volume
+      double mesh_measure = 0; // area or volume
       double proximal_p = 0;
 
       Eigen::VectorXi frame_ids;

include/igl/writePLY.cpp

@@ -147,7 +147,7 @@ bool writePLY(
         _ed.resize(ED.size());
         Eigen::Map<Eigen::Matrix<EDScalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > >( &_ed[0], ED.rows(), ED.cols() ) = ED;
 
-        file.add_properties_to_element("face", FDheader,
+        file.add_properties_to_element("edge", EDheader,
             tynyply_type<EDScalar>(), ED.rows(), reinterpret_cast<uint8_t*>( &_ed[0] ), tinyply::Type::INVALID, 0);
     }
 

tests/include/igl/linprog.cpp

@@ -0,0 +1,41 @@
+#include <test_common.h>
+#include <igl/linprog.h>
+
+TEST_CASE("linprog: 2D-inequality", "[igl]" )
+{
+  Eigen::VectorXd f(2); f << -2, -1;
+  Eigen::VectorXd b(3); b << 8, 14, 4;
+  Eigen::MatrixXd A(3, 2); A << 2, -1, 1, 2, -1, 1;
+  Eigen::MatrixXd B = Eigen::MatrixXd(0,2);
+  Eigen::VectorXd c = Eigen::VectorXd(0);
+  Eigen::VectorXd x;
+  igl::linprog(f, A, b, B, c, x);
+  Eigen::VectorXd x_correct(2); x_correct<<6,4;
+  test_common::assert_near( x, x_correct, 1e-10);
+}
+
+TEST_CASE("linprog: 2D-inequality+2-equality", "[igl]" )
+{
+  Eigen::VectorXd f(2); f << -2, -1;
+  Eigen::VectorXd b(3); b << 8, 14, 4;
+  Eigen::MatrixXd A(3, 2); A << 2, -1, 1, 2, -1, 1;
+  Eigen::MatrixXd B = Eigen::MatrixXd(2,2);B<<1,0,0,1;
+  Eigen::VectorXd c = Eigen::VectorXd(2);c<<4,4;
+  Eigen::VectorXd x;
+  igl::linprog(f, A, b, B, c, x);
+  Eigen::VectorXd x_correct(2); x_correct<<4,4;
+  test_common::assert_near( x, x_correct, 1e-10);
+}
+
+TEST_CASE("linprog: 2D-inequality+1-equality", "[igl]" )
+{
+  Eigen::VectorXd f(2); f << -2, -1;
+  Eigen::VectorXd b(3); b << 8, 14, 4;
+  Eigen::MatrixXd A(3, 2); A << 2, -1, 1, 2, -1, 1;
+  Eigen::MatrixXd B = Eigen::MatrixXd(1,2);B<<1,0;
+  Eigen::VectorXd c = Eigen::VectorXd(1);c<<4;
+  Eigen::VectorXd x;
+  igl::linprog(f, A, b, B, c, x);
+  Eigen::VectorXd x_correct(2); x_correct<<4,5;
+  test_common::assert_near( x, x_correct, 1e-10);
+}

tests/include/igl/readPLY.cpp

@@ -94,6 +94,10 @@ TEST_CASE("readPLY: quad_cube.ply", "[igl]")
 "property float t\n"
 "element face 6\n"
 "property list uchar uint vertex_indices\n"
+"property uchar red\n"
+"property uchar green\n"
+"property uchar blue\n"
+"property uchar alpha\n"
 "end_header\n"
 "1.000000 1.000000 1.000000 0.000000 -0.000000 1.000000 0.625000 0.500000\n"
 "-1.000000 1.000000 1.000000 0.000000 -0.000000 1.000000 0.875000 0.500000\n"
@@ -119,12 +123,13 @@ TEST_CASE("readPLY: quad_cube.ply", "[igl]")
 "-1.000000 1.000000 1.000000 0.000000 1.000000 0.000000 0.625000 0.250000\n"
 "1.000000 1.000000 1.000000 0.000000 1.000000 0.000000 0.625000 0.500000\n"
 "1.000000 1.000000 -1.000000 0.000000 1.000000 0.000000 0.375000 0.500000\n"
-"4 0 1 2 3\n"
-"4 4 5 6 7\n"
-"4 8 9 10 11\n"
-"4 12 13 14 15\n"
-"4 16 17 18 19\n"
-"4 20 21 22 23\n";
+"4 0 1 2 3 0 0 255 255\n"
+"4 4 5 6 7 0 0 255 255\n"
+"4 8 9 10 11 0 0 255 255\n"
+"4 12 13 14 15 0 0 255 255\n"
+"4 16 17 18 19 0 0 255 255\n"
+"4 20 21 22 23 0 0 255 255\n";
+
 
     std::ofstream f("quad_cube.ply");
     f.write(ply_quad_cube,strlen(ply_quad_cube));
@@ -137,4 +142,47 @@ TEST_CASE("readPLY: quad_cube.ply", "[igl]")
     REQUIRE (V.cols() == 3);
     REQUIRE (F.rows() == 6);
     REQUIRE (F.cols() == 4);
+
+    Eigen::MatrixXd N,UV,VD,FD,ED;
+    Eigen::MatrixXi E;
+    std::vector<std::string> headerV;
+    std::vector<std::string> headerF, headerE, comments;
+ 
+    REQUIRE (igl::readPLY("quad_cube.ply", V, F, E, N, UV, 
+              VD, headerV,
+              FD, headerF, 
+              ED, headerE, 
+              comments));
+    
+    REQUIRE (V.rows() == 24);
+    REQUIRE (V.cols() == 3);
+    REQUIRE (F.rows() == 6);
+    REQUIRE (F.cols() == 4);
+
+    REQUIRE (E.rows() == 0);
+    REQUIRE (E.cols() == 0);
+
+    REQUIRE (N.rows() == 24);
+    REQUIRE (N.cols() == 3);
+
+    REQUIRE (UV.rows() == 24);
+    REQUIRE (UV.cols() == 2);
+
+    REQUIRE (VD.rows() == 0);
+    REQUIRE (VD.cols() == 0);
+    REQUIRE (headerV.empty());
+
+    REQUIRE (FD.rows() == 6);
+    REQUIRE (FD.cols() == 4);
+    REQUIRE (headerF.size() == 4);
+
+    REQUIRE (headerF[0] == "red");
+    REQUIRE (headerF[1] == "green");
+    REQUIRE (headerF[2] == "blue");
+    REQUIRE (headerF[3] == "alpha");
+
+    REQUIRE (ED.rows() == 0);
+    REQUIRE (ED.cols() == 0);
+
+    REQUIRE (headerE.size() == 0);
 }

tests/include/igl/writePLY.cpp

@@ -1,6 +1,7 @@
 #include <test_common.h>
 #include <igl/readPLY.h>
 #include <igl/writePLY.h>
+#include <igl/edges.h>
 #include <fstream>
 #include <string>
 #include <vector>
@@ -31,12 +32,37 @@ TEST_CASE("writePLY: bunny.ply", "[igl]")
     Eigen::VectorXd face_data(F1.rows());
     for(size_t i=0;i<F1.rows();++i)
         face_data(i)=(double)i;
-    Eigen::MatrixXd FD2(F1.rows(),FD1.cols()+1);
+
+
+
     // there is no face data in the input file
+    Eigen::MatrixXd FD2(F1.rows(),FD1.cols()+1);
     FD2<<face_data;
 
+    //input file have no edge data
+    REQUIRE (E1.rows() == 0);
+    REQUIRE (E1.cols() == 0);
+    REQUIRE (ED1.rows() == 0);
+    REQUIRE (Eheader1.empty());    
+
+    Eigen::MatrixXi E2;
+    std::vector<std::string> Eheader2;
+
+    //generate edges
+    igl::edges(F1,E2);
+
+    //generate edge data
+    Eheader2.push_back("edge_data");
+    Eigen::VectorXd edge_data(E2.rows());
+    for(size_t i=0;i<E2.rows();++i)
+        edge_data(i)=(double)i;
+    
+    // there is no edge data in the input file
+    Eigen::MatrixXd ED2(E2.rows(),1);
+    ED2<<edge_data;
+
     // test that saving preserves all the data, including new data column
-    REQUIRE (igl::writePLY("writePLY_test_bunny.ply", V1, F1, E1, N1, UV1, VD2, Vheader1, FD2,Fheader1, ED1, Eheader1, comments1, igl::FileEncoding::Binary));
+    REQUIRE (igl::writePLY("writePLY_test_bunny.ply", V1, F1, E2, N1, UV1, VD2, Vheader1, FD2,Fheader1, ED2, Eheader2, comments1, igl::FileEncoding::Binary));
 
     Eigen::MatrixXd V,N,UV,VD,FD,ED;
     Eigen::MatrixXi F,E;
@@ -49,9 +75,10 @@ TEST_CASE("writePLY: bunny.ply", "[igl]")
     REQUIRE (V.cols() == 3);
     REQUIRE (F.rows() == 69451);
     REQUIRE (F.cols() == 3);
-    // no edge data
-    REQUIRE (E.rows() == 0);
-    REQUIRE (E.cols() == 0);
+
+    // generated edges
+    REQUIRE (E.rows() == E2.rows());
+    REQUIRE (E.cols() == E2.cols());
 
     // no normals or texture coordinates
     REQUIRE (N.rows() == 0);
@@ -72,7 +99,7 @@ TEST_CASE("writePLY: bunny.ply", "[igl]")
     REQUIRE (Vheader[1] == "intensity" );
     REQUIRE (Vheader[2] == "dummy_data" );
 
-    // Face should have only one column
+    // Face datashould have only one column
     REQUIRE (Fheader.size() == 1);
     REQUIRE (Fheader[0] == "face_data" );
     REQUIRE (FD.rows() == F.rows());
@@ -82,9 +109,16 @@ TEST_CASE("writePLY: bunny.ply", "[igl]")
     for(size_t i=0;i<F.rows();++i)
         REQUIRE (FD(i,0) == (double)i);
 
-    REQUIRE (ED.rows() == 0);
-    REQUIRE (ED.cols() == 0);
-    REQUIRE (Eheader.size() == 0);
+    // Edge data should have only one column
+    REQUIRE (Eheader.size() == 1);
+    REQUIRE (Eheader[0] == "edge_data" );
+    REQUIRE (ED.rows() == E.rows());
+    REQUIRE (ED.cols() == 1);
+
+    // the dummy column contents check
+    for(size_t i=0;i<E.rows();++i)
+        REQUIRE (ED(i,0) == (double)i);
+
 
     // there are comments
     REQUIRE (comments.size() == 2);