bad types in lex (#2544)

* swept volume overload

* swept volume: take transform list, templatize, expose SignedDistanceType

Replace the transform(t)+steps interface of swept_volume,
swept_volume_signed_distance and swept_volume_bounding_box with a list of
rigid transformations passed directly. The transform(t)+steps overloads are
removed rather than kept.

Templatize on Eigen::MatrixBase inputs / Eigen::PlainObjectBase outputs per
libigl style. The transform list is templated on both scalar and allocator so
std::vector<Eigen::Affine3d> and the aligned_allocator spelling both bind.

Expose SignedDistanceType on swept_volume and swept_volume_signed_distance,
dispatching like signed_distance_3 (pseudonormal / winding number / fast
winding number / unsigned) with the precomputation hoisted out of the
per-time-step loop.

swept_volume's isolevel is now a distance (typename DerivedV::Scalar) rather
than a count of grid cells. The padding cancelled out of the grid spacing
already (h == diag/(grid_res-1) regardless of pad), so pad is now derived from
the requested distance instead. Passing isolevel = k*h for integer k
reproduces the old isolevel_grid = k grid exactly. The doc had claimed
"distance level to be contoured" all along while the size_t type and
isolevel_grid implementation meant cells.

707_SweptVolume preloads the motion as a list of transforms sampled uniformly
over t in [0,1] and passes that to swept_volume. Its isolevel of 0.1 had been
truncating to 0 by the size_t conversion, so the tutorial was dilating by
nothing; it is now 10% of the bunny's largest side.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

---------

Co-authored-by: Claude Opus 5 <noreply@anthropic.com>

.github/workflows/continuous.yml

@@ -25,14 +25,14 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        os: [ubuntu-24.04, macos-latest]
+        os: [ubuntu-24.04, macos-15]
         config: [Release]
         build-params: [ {static:  ON, tutorials:  ON, tests:  ON }, {static: OFF, tutorials: OFF, tests:  ON }, {static: OFF, tutorials:  ON, tests: OFF }]
     env:
       IGL_NUM_THREADS: 1  # See https://github.com/libigl/libigl/pull/996
     steps:
       - name: Checkout repository
-        uses: actions/checkout@v1
+        uses: actions/checkout@v4
         with:
           fetch-depth: 10
 
@@ -56,7 +56,7 @@ jobs:
         if: runner.os == 'macOS'
         uses: maxim-lobanov/setup-xcode@v1
         with:
-          xcode-version: '16.3.0'
+          xcode-version: latest-stable
 
       - name: Cache Build
         id: cache-build
@@ -118,31 +118,15 @@ jobs:
       ]
     steps:
       - name: Checkout repository
-        uses: actions/checkout@v1
+        uses: actions/checkout@v4
         with:
           fetch-depth: 10
 
       - name: Install Ninja
         uses: seanmiddleditch/gha-setup-ninja@master
 
-      - name: Set env variable for sccache
-        run: |
-          echo "appdata=$env:LOCALAPPDATA" >> ${env:GITHUB_ENV}
-
-      - name: Cache build
-        id: cache-build
-        uses: actions/cache@v4
-        with:
-          path: ${{ env.appdata }}\Mozilla\sccache
-          key: ${{ runner.os }}-${{ matrix.config }}-${{ matrix.build-params.static }}-cache
-
-      - name: Prepare sccache
-        run: |
-          iwr -useb 'https://raw.githubusercontent.com/scoopinstaller/install/master/install.ps1' -outfile 'install.ps1'
-          .\install.ps1 -RunAsAdmin
-          scoop install sccache --global
-          # Scoop modifies the PATH so we make it available for the next steps of the job
-          echo "${env:PATH}" >> ${env:GITHUB_PATH}
+      - name: Install sccache
+        uses: mozilla-actions/sccache-action@v0.0.9
 
       - name: Configure and build
         shell: cmd
@@ -174,4 +158,3 @@ jobs:
 
       - name: Tests
         run: cd build; ctest --show-only; ctest --verbose -j2
-

.gitignore

@@ -51,3 +51,4 @@ LibiglOptions.cmake
 dox/
 latex/
 scripts/
+CLAUDE.md

CMakeLists.txt

@@ -99,6 +99,7 @@ option(LIBIGL_USE_STATIC_LIBRARY "Use libigl as static library" ${LIBIGL_TOPLEVE
 
 # Permissive modules. These modules are available under MPL2 license, and their dependencies are available
 # under a permissive or public domain license.
+option(LIBIGL_CYCODEBASE       "Build target igl::cycodebase"       ${LIBIGL_TOPLEVEL_PROJECT})
 option(LIBIGL_EMBREE           "Build target igl::embree"           ${LIBIGL_TOPLEVEL_PROJECT})
 option(LIBIGL_GLFW             "Build target igl::glfw"             ${LIBIGL_TOPLEVEL_PROJECT})
 option(LIBIGL_IMGUI            "Build target igl::imgui"            ${LIBIGL_TOPLEVEL_PROJECT})

cmake/igl/modules/cycodebase.cmake

@@ -0,0 +1,26 @@
+# 1. Define module
+igl_add_library(igl_cycodebase)
+
+# 2. Include headers
+include(GNUInstallDirs)
+target_include_directories(igl_cycodebase ${IGL_SCOPE}
+    $<BUILD_INTERFACE:${libigl_SOURCE_DIR}/include>
+    $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
+)
+
+# 3. Target sources
+file(GLOB INC_FILES "${libigl_SOURCE_DIR}/include/igl/cycodebase/*.h")
+file(GLOB SRC_FILES "${libigl_SOURCE_DIR}/include/igl/cycodebase/*.cpp")
+igl_target_sources(igl_cycodebase ${INC_FILES} ${SRC_FILES})
+
+# 4. Dependencies
+include(cycodebase)
+target_link_libraries(igl_cycodebase ${IGL_SCOPE}
+    igl::core
+    cyCodeBase::cyCodeBase
+)
+
+# 5. Unit tests
+file(GLOB SRC_FILES "${libigl_SOURCE_DIR}/tests/include/igl/cycodebase/*.cpp")
+igl_add_test(igl_cycodebase ${SRC_FILES})
+

cmake/libigl.cmake

@@ -13,6 +13,7 @@ include(igl_windows)
 
 # Libigl permissive modules
 igl_include(core)
+igl_include_optional(cycodebase)
 igl_include_optional(embree)
 igl_include_optional(opengl)
 igl_include_optional(glfw)

cmake/recipes/external/cycodebase.cmake

@@ -0,0 +1,18 @@
+if(TARGET cycodebase::cycodebase)
+    return()
+endif()
+
+FetchContent_Declare(
+  cyCodeBase
+  GIT_REPOSITORY https://github.com/cemyuksel/cyCodeBase/
+  GIT_TAG e36f3cffca65eb12a8a071f0443128b7de6ed75d
+)
+FetchContent_Populate(cyCodeBase)
+add_library(cyCodeBase_interface INTERFACE)
+target_include_directories(cyCodeBase_interface INTERFACE ${cycodebase_SOURCE_DIR})
+
+if(NOT (CMAKE_SYSTEM_PROCESSOR MATCHES "x86_64|AMD64|i[3-6]86"))
+  target_compile_definitions(cyCodeBase_interface INTERFACE CY_NO_INTRIN_H)
+endif()
+
+add_library(cyCodeBase::cyCodeBase ALIAS cyCodeBase_interface)

cmake/recipes/external/eigen.cmake

@@ -8,7 +8,7 @@ include(FetchContent)
 FetchContent_Declare(
     eigen
     GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git
-    GIT_TAG tags/3.4.0
+    GIT_TAG tags/5.0.1
     GIT_SHALLOW TRUE
 )
 FetchContent_GetProperties(eigen)

cmake/recipes/external/libigl_tutorial_data.cmake

@@ -8,7 +8,7 @@ include(FetchContent)
 FetchContent_Declare(
     libigl_tutorial_data
     GIT_REPOSITORY https://github.com/libigl/libigl-tutorial-data
-    GIT_TAG        6700bf49000cd64199835fb40323e5ca9c7796ab
+    GIT_TAG        644dd4104843b6d736745d9dafbd70bf8d175648
 )
 FetchContent_MakeAvailable(libigl_tutorial_data)
 

include/igl/FastWindingNumberForSoups.h

@@ -281,7 +281,7 @@ static constexpr inline fpreal32 SYSmax(fpreal32 a, fpreal32 b)	{ return h_max(a
 static constexpr inline fpreal64 SYSmin(fpreal64 a, fpreal64 b)	{ return h_min(a,b); }
 static constexpr inline fpreal64 SYSmax(fpreal64 a, fpreal64 b)	{ return h_max(a,b); }
 
-// Some systems have size_t as a seperate type from uint.  Some don't.
+// Some systems have size_t as a separate type from uint.  Some don't.
 #if (defined(LINUX) && defined(IA64)) || defined(MBSD)
 static constexpr inline size_t SYSmin(size_t a, size_t b)		{ return h_min(a,b); }
 static constexpr inline size_t SYSmax(size_t a, size_t b)		{ return h_max(a,b); }
@@ -742,408 +742,408 @@ vm_allbits(const v4si &a)
 
 #endif
 #endif
-#pragma once
-#ifndef __SSE__
-#ifndef __VM_SIMDFunc__
-#define __VM_SIMDFunc__
-
+#pragma once
+#ifndef __SSE__
+#ifndef __VM_SIMDFunc__
+#define __VM_SIMDFunc__
+
+
+
+#include <cmath>
 
-
-#include <cmath>
-
 namespace igl {
   /// @private
-  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
+  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.
  *
@@ -3758,7 +3758,7 @@ namespace UT_Thread { inline int getNumProcessors() {
 /////
 ///// Requirements for the Range functor are:
 /////   - the requirements of the tbb Range Concept
-/////   - UT_estimatorNumItems<Range> must return the the estimated number of work items
+/////   - UT_estimatorNumItems<Range> must return the estimated number of work items
 /////     for the range. When Range::size() is not the correct estimate, then a 
 /////     (partial) specialization of UT_estimatorNumItemsimatorRange must be provided
 /////     for the type Range.

include/igl/HalfEdgeIterator.cpp

@@ -8,6 +8,8 @@
 
 #include "HalfEdgeIterator.h"
 
+#include <cassert>
+
 template <typename DerivedF, typename DerivedFF, typename DerivedFFi>
 IGL_INLINE igl::HalfEdgeIterator<DerivedF,DerivedFF,DerivedFFi>::HalfEdgeIterator(
     const Eigen::MatrixBase<DerivedF>& _F,

include/igl/MshLoader.cpp

@@ -114,10 +114,31 @@ IGL_INLINE igl::MshLoader::MshLoader(const std::string &filename) {
     fin.close();
 }
 
+IGL_INLINE int igl::MshLoader::node_dense_index(int node_tag) const {
+    const auto it = m_node_tag_to_dense.find(node_tag);
+    if (it == m_node_tag_to_dense.end()) {
+        std::stringstream err_msg;
+        err_msg << "Unknown node tag: " << node_tag;
+        throw std::runtime_error(err_msg.str());
+    }
+	return it->second;
+}
+
+IGL_INLINE int igl::MshLoader::element_dense_index(int elem_tag) const {
+	const auto it = m_element_tag_to_dense.find(elem_tag);
+    if (it == m_element_tag_to_dense.end()) {
+		std::stringstream err_msg;
+        err_msg << "Unknown element tag: " << elem_tag;
+		throw std::runtime_error(err_msg.str());
+    }
+	return it->second;
+}
+
 IGL_INLINE void igl::MshLoader::parse_nodes(std::ifstream& fin) {
     size_t num_nodes;
     fin >> num_nodes;
     m_nodes.resize(num_nodes*3);
+    m_node_tag_to_dense.clear();
 
     if (m_binary) {
 		size_t stride = (4+3*m_data_size);
@@ -127,23 +148,37 @@ IGL_INLINE void igl::MshLoader::parse_nodes(std::ifstream& fin) {
         fin.read(data, num_bytes);
 
         for (size_t i=0; i<num_nodes; i++) {
-            int node_idx;
-			memcpy(&node_idx, data+i*stride, sizeof(int));
-			node_idx-=1;
-			// directly move into vector storage
-			// this works only when m_data_size==sizeof(Float)==sizeof(double)
-			memcpy(&m_nodes[node_idx*3], data+i*stride + 4, m_data_size*3);
+            int node_tag;
+			memcpy(&node_tag, data+i*stride, sizeof(int));
+
+            if (node_tag <= 0) {
+                throw std::runtime_error("Invalid node tag");
+            }
+            if (m_node_tag_to_dense.find(node_tag) != m_node_tag_to_dense.end()) {
+                throw std::runtime_error("Duplicate node tag");
+            }
+            m_node_tag_to_dense[node_tag] = static_cast<int>(i);
+            // directly move into vector storage
+            // this works only when m_data_size==sizeof(Float)==sizeof(double)
+			memcpy(&m_nodes[i*3], data+i*stride + 4, m_data_size*3);
         }
         delete [] data;
     } else {
-        int node_idx;
+        int node_tag;
         for (size_t i=0; i<num_nodes; i++) {
-            fin >> node_idx;
-            node_idx -= 1;
+            fin >> node_tag;
+
+            if (node_tag <= 0) {
+				throw std::runtime_error("Invalid node tag");
+            }
+            if (m_node_tag_to_dense.find(node_tag) != m_node_tag_to_dense.end()) {
+                throw std::runtime_error("Duplicate node tag");
+            }
+            m_node_tag_to_dense[node_tag] = static_cast<int>(i);
             // here it's 3D node explicitly
-            fin >> m_nodes[node_idx*3]
-                >> m_nodes[node_idx*3+1]
-                >> m_nodes[node_idx*3+2];
+            fin >> m_nodes[i*3]
+                >> m_nodes[i*3+1]
+                >> m_nodes[i*3+2];
         }
     }
 }
@@ -152,6 +187,7 @@ IGL_INLINE void igl::MshLoader::parse_elements(std::ifstream& fin) {
     m_elements_tags.resize(2); //hardcoded to have 2 tags
     size_t num_elements;
     fin >> num_elements;
+    m_element_tag_to_dense.clear();
 
     size_t nodes_per_element;
 
@@ -168,15 +204,24 @@ IGL_INLINE void igl::MshLoader::parse_elements(std::ifstream& fin) {
 
             // store node info
             for (size_t i=0; i<num_elems; i++) {
-                int elem_idx;
+                int elem_tag;
 
                 // all elements in the segment share the same elem_type and number of nodes per element
                 m_elements_types.push_back(elem_type);
                 m_elements_lengths.push_back(nodes_per_element);
 
-                fin.read((char*)&elem_idx, sizeof(int));
-                elem_idx -= 1;
-                m_elements_ids.push_back(elem_idx);
+                fin.read((char*)&elem_tag, sizeof(int));
+
+                if (elem_tag <= 0) {
+                    throw std::runtime_error("Invalid element tag");
+                }
+                if (m_element_tag_to_dense.find(elem_tag) != m_element_tag_to_dense.end()) {
+                    throw std::runtime_error("Duplicate element tag");
+                }
+                m_element_tag_to_dense[elem_tag] = static_cast<int>(m_elements_ids.size());
+
+                elem_tag -= 1;
+                m_elements_ids.push_back(elem_tag);
 
                 // read first two tags
                 for (size_t j=0; j<num_tags; j++) {
@@ -191,10 +236,10 @@ IGL_INLINE void igl::MshLoader::parse_elements(std::ifstream& fin) {
                 m_elements_nodes_idx.push_back(m_elements.size());
                 // Element values.
                 for (size_t j=0; j<nodes_per_element; j++) {
-                    int idx;
-                    fin.read((char*)&idx, sizeof(int));
+                    int node_tag;
+                    fin.read((char*)&node_tag, sizeof(int));
                     
-                    m_elements.push_back(idx-1);
+                    m_elements.push_back(node_dense_index(node_tag));
                 }
             }
             elem_read += num_elems;
@@ -202,8 +247,16 @@ IGL_INLINE void igl::MshLoader::parse_elements(std::ifstream& fin) {
     } else {
         for (size_t i=0; i<num_elements; i++) {
             // Parse per element header
-            int elem_num, elem_type, num_tags;
-            fin >> elem_num >> elem_type >> num_tags;
+            int elem_tag, elem_type, num_tags;
+            fin >> elem_tag >> elem_type >> num_tags;
+
+            if (elem_tag <= 0) {
+                throw std::runtime_error("Invalid element tag");
+            }
+            if (m_element_tag_to_dense.find(elem_tag) != m_element_tag_to_dense.end()) {
+                throw std::runtime_error("Duplicate element tag");
+            }
+            m_element_tag_to_dense[elem_tag] = static_cast<int>(m_elements_ids.size());
 
             // read tags.
             for (size_t j=0; j<num_tags; j++) {
@@ -218,14 +271,14 @@ IGL_INLINE void igl::MshLoader::parse_elements(std::ifstream& fin) {
             m_elements_types.push_back(elem_type);
             m_elements_lengths.push_back(nodes_per_element);
 
-            elem_num -= 1;
-            m_elements_ids.push_back(elem_num);
+            elem_tag -= 1;
+            m_elements_ids.push_back(elem_tag);
             m_elements_nodes_idx.push_back(m_elements.size());
             // Parse node idx.
             for (size_t j=0; j<nodes_per_element; j++) {
-                int idx;
-                fin >> idx;
-                m_elements.push_back(idx-1); // msh index starts from 1.
+                int node_tag;
+                fin >> node_tag;
+                m_elements.push_back(node_dense_index(node_tag)); // msh index starts from 1.
             }
         }
     }
@@ -274,7 +327,7 @@ IGL_INLINE void igl::MshLoader::parse_node_field( std::ifstream& fin ) {
     int num_components    = int_tags[1];
     int num_entries       = int_tags[2];
 
-    std::vector<Float> field( num_entries*num_components );
+    std::vector<Float> field((m_nodes.size()/3)*num_components);
 
     if (m_binary) {
         size_t num_bytes = (num_components * m_data_size + 4) * num_entries;
@@ -282,23 +335,20 @@ IGL_INLINE void igl::MshLoader::parse_node_field( std::ifstream& fin ) {
         igl::_msh_eat_white_space(fin);
         fin.read(data, num_bytes);
         for (size_t i=0; i<num_entries; i++) {
-			int node_idx;
-			memcpy(&node_idx,&data[i*(4+num_components*m_data_size)],4);
-			
-            if(node_idx<1) throw std::runtime_error("Negative or zero index");
-            node_idx -= 1;
-			
-            if(node_idx>=num_entries) throw std::runtime_error("Index too big");
+			int node_tag;
+			memcpy(&node_tag,&data[i*(4+num_components*m_data_size)],4);
+
+            const int node_idx = node_dense_index(node_tag);
             size_t base_idx = i*(4+num_components*m_data_size) + 4;
             // TODO: make this work when m_data_size != sizeof(double) ?
 			memcpy(&field[node_idx*num_components], &data[base_idx], num_components*m_data_size);
         }
         delete [] data;
     } else {
-        int node_idx;
+        int node_tag;
         for (size_t i=0; i<num_entries; i++) {
-            fin >> node_idx;
-            node_idx -= 1;
+            fin >> node_tag;
+            const int node_idx = node_dense_index(node_tag);
             for (size_t j=0; j<num_components; j++) {
                 fin >> field[node_idx*num_components+j];
             }
@@ -346,7 +396,7 @@ IGL_INLINE void igl::MshLoader::parse_element_field(std::ifstream& fin) {
     std::string fieldname = str_tags[0];
     int num_components = int_tags[1];
     int num_entries = int_tags[2];
-    std::vector<Float> field(num_entries*num_components);
+    std::vector<Float> field(m_elements_ids.size()*num_components);
 
     if (m_binary) {
         size_t num_bytes = (num_components * m_data_size + 4) * num_entries;
@@ -354,20 +404,20 @@ IGL_INLINE void igl::MshLoader::parse_element_field(std::ifstream& fin) {
         igl::_msh_eat_white_space(fin);
         fin.read(data, num_bytes);
         for (int i=0; i<num_entries; i++) {
-			int elem_idx;
+			int elem_tag;
 			// works with sizeof(int)==4
-			memcpy(&elem_idx, &data[i*(4+num_components*m_data_size)],4);
-            elem_idx -= 1;
+			memcpy(&elem_tag, &data[i*(4+num_components*m_data_size)],4);
+            const int elem_idx = element_dense_index(elem_tag);
 			
 			// directly copy data into vector storage space
 			memcpy(&field[elem_idx*num_components], &data[i*(4+num_components*m_data_size) + 4], m_data_size*num_components);
         }
         delete [] data;
     } else {
-        int elem_idx;
+        int elem_tag;
         for (size_t i=0; i<num_entries; i++) {
-            fin >> elem_idx;
-            elem_idx -= 1;
+            fin >> elem_tag;
+            const int elem_idx = element_dense_index(elem_tag);
             for (size_t j=0; j<num_components; j++) {
                 fin >> field[elem_idx*num_components+j];
             }

include/igl/MshLoader.h

@@ -15,6 +15,7 @@
 #include <string>
 #include <vector>
 #include <algorithm>
+#include <unordered_map>
 
 namespace igl {
 
@@ -155,11 +156,15 @@ class MshLoader {
         void parse_element_field(std::ifstream& fin);
         void parse_unknown_field(std::ifstream& fin,
                 const std::string& fieldname);
+        int node_dense_index(int node_tag) const;
+        int element_dense_index(int elem_tag) const;
 
     private:
         bool   m_binary;
         size_t m_data_size;
-        
+        std::unordered_map<int, int> m_node_tag_to_dense;
+        std::unordered_map<int, int> m_element_tag_to_dense;
+
         FloatVector m_nodes;    // len x 3 vector 
 
         IndexVector m_elements; // linear array for nodes corresponding to each element 

include/igl/Orientation.h

@@ -7,23 +7,21 @@
 // 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/.
 #pragma once
-#ifndef IGL_PREDICATES_ORIENTATION_H
-#define IGL_PREDICATES_ORIENTATION_H
+#ifndef IGL_ORIENTATION_H
+#define IGL_ORIENTATION_H
 
-#include "../igl_inline.h"
+#include "igl_inline.h"
 #include <Eigen/Core>
 
 namespace igl {
-  namespace predicates {
-    /// Types of orientations and other predicate results.
-    ///
-    /// \fileinfo
-    enum class Orientation {
-      POSITIVE=1, INSIDE=1,
-      NEGATIVE=-1, OUTSIDE=-1,
-      COLLINEAR=0, COPLANAR=0, COCIRCULAR=0, COSPHERICAL=0, DEGENERATE=0
-    };
-  }
+  /// Types of orientations and other predicate results.
+  ///
+  /// \fileinfo
+  enum class Orientation {
+    POSITIVE=1, INSIDE=1,
+    NEGATIVE=-1, OUTSIDE=-1,
+    COLLINEAR=0, COPLANAR=0, COCIRCULAR=0, COSPHERICAL=0, DEGENERATE=0
+  };
 }
 
 

include/igl/SphereMeshWedge.cpp

@@ -8,6 +8,7 @@
 #include "SphereMeshWedge.h"
 #include "round_cone_signed_distance.h"
 #include "sign.h"
+#include <cassert>
 #include <Eigen/QR>
 #include <Eigen/Geometry>
 

include/igl/WindingNumberAntipodalScene.h

@@ -0,0 +1,259 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2026 Philip Trettner <trettner@shapedcode.com>, Cedric Martens <cedric.martens@umontreal.ca>
+//
+// 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_WINDINGNUMBERANTIPODALSCENE_H
+#define IGL_WINDINGNUMBERANTIPODALSCENE_H
+
+#include "PI.h"
+#include "parallel_for.h"
+
+#include <Eigen/Core>
+
+#include <cassert>
+#include <cmath>
+#include <cstdint>
+#include <unordered_map>
+#include <vector>
+
+namespace igl
+{
+  /// Precomputed scene for the Antipodal Method generalized winding number
+  ///
+  /// The scene stores only the weighted open-boundary edges of the input
+  /// triangle mesh and a fixed antipodal reference direction `x0`. Closed
+  /// (manifold) meshes have an empty boundary and the fractional term is
+  /// exactly zero.
+  ///
+  /// Querying the winding number additionally requires an `Intersector` that
+  /// returns the signed ray-mesh crossing count along `(p, x0)` over the
+  /// original mesh; the scene itself is intersector-agnostic. See
+  /// `igl::embree::EmbreeIntersector::signedIntersectionsRay` for an
+  /// optimized concrete implementation.
+  ///
+  /// ### Intersector concept
+  /// A type `I` satisfies the concept when it exposes:
+  ///   - `using OriginType    = ...;`     (3D row vector type)
+  ///   - `using DirectionType = ...;`     (3D row vector type)
+  ///   - `int signedIntersectionsRay(
+  ///         OriginType origin, DirectionType direction,
+  ///         /* defaulted tnear, tfar, mask */) const;`
+  ///
+  /// `winding_number` casts query point/direction to the intersector's types
+  /// at the call site, so a `double`-precision scene against a float-only
+  /// intersector works without an adaptor.
+  template <typename Scalar>
+  class WindingNumberAntipodalScene
+  {
+  public:
+    using Point     = Eigen::Matrix<Scalar, 1, 3>;
+    using Direction = Eigen::Matrix<Scalar, 1, 3>;
+
+  private:
+    struct WeightedSeg
+    {
+      Point a;
+      Point b;
+      Scalar w;
+    };
+
+  public:
+    /// Build a scene from a 3D triangle mesh.
+    ///
+    /// @param[in] V   #V by 3 list of vertex positions
+    /// @param[in] F   #F by 3 list of triangle indices
+    /// @param[in] x0  unit reference direction (defaults to a fixed non-axis-aligned vector)
+    template <typename DerivedV, typename DerivedF>
+    WindingNumberAntipodalScene(
+      const Eigen::MatrixBase<DerivedV> & V,
+      const Eigen::MatrixBase<DerivedF> & F,
+      const Direction & x0 = default_x0())
+      : m_x0(x0), m_face_count(static_cast<size_t>(F.rows()))
+    {
+      assert(V.cols() == 3 && "WindingNumberAntipodalScene: only 3D vertex positions are supported");
+      assert(F.cols() == 3 && "WindingNumberAntipodalScene: only triangle meshes are supported");
+      build_boundary_segments(V, F, m_boundary);
+    }
+
+    /// Single-point query: full generalized winding number at `p`
+    /// (fractional + signed integer crossings).
+    template <typename Intersector, typename Derivedp>
+    Scalar winding_number(
+      const Intersector & intersector,
+      const Eigen::MatrixBase<Derivedp> & p) const
+    {
+      const Point pp(static_cast<Scalar>(p(0)),
+                     static_cast<Scalar>(p(1)),
+                     static_cast<Scalar>(p(2)));
+      const Direction x1 = -m_x0;
+      Scalar area = Scalar(0);
+      for (const auto & ws : m_boundary)
+      {
+        const Point v0 = ws.a - pp;
+        const Point v1 = ws.b - pp;
+        area += ws.w * half_solid_angle_unorm(x1, v0, v1);
+      }
+      const Scalar frac = area / (Scalar(2) * Scalar(igl::PI));
+
+      using IO = typename Intersector::OriginType;
+      using ID = typename Intersector::DirectionType;
+      const IO io(static_cast<typename IO::Scalar>(pp(0)),
+                  static_cast<typename IO::Scalar>(pp(1)),
+                  static_cast<typename IO::Scalar>(pp(2)));
+      const ID id(static_cast<typename ID::Scalar>(m_x0(0)),
+                  static_cast<typename ID::Scalar>(m_x0(1)),
+                  static_cast<typename ID::Scalar>(m_x0(2)));
+      const int c = intersector.signedIntersectionsRay(io, id);
+
+      return frac + Scalar(c);
+    }
+
+    /// Batch query, parallelized via `igl::parallel_for`.
+    ///
+    /// @param[in]  intersector  Concept-compatible intersector built over the
+    ///                          same mesh used to construct the scene. Must
+    ///                          be safe to query concurrently.
+    /// @param[in]  O            #O by 3 list of query points
+    /// @param[out] W            #O by 1 list of winding numbers
+    template <typename Intersector, typename DerivedO, typename DerivedW>
+    void winding_number(
+      const Intersector & intersector,
+      const Eigen::MatrixBase<DerivedO> & O,
+      Eigen::PlainObjectBase<DerivedW> & W) const
+    {
+      W.resize(O.rows(), 1);
+
+      // Adaptive parallel-for threshold.
+      //
+      // The libigl thread pool has a roughly fixed ~1 ms TOTAL overhead per
+      // parallel_for invocation (not per iteration). So we only spawn the
+      // pool when the WHOLE batch is expected to take ≥ 1 ms.
+      //
+      // Per-query work heuristic:
+      //   t_q ≈ 50 * B + 100 * sqrt(F)  ns
+      // with B = boundary segment count, F = triangle count. Pool wins once
+      //   t_q * O > 10^6 ns  ⇒  O > 10^6 / t_q
+      // which is exactly parallel_for's `min_parallel` semantics.
+      //
+      // (This is a rough heuristic and should be revisited once parallel_for becomes lower-overhead)
+      const double t_q_ns =
+          50.0  * static_cast<double>(m_boundary.size()) +
+          100.0 * std::sqrt(static_cast<double>(m_face_count));
+      const size_t min_parallel = static_cast<size_t>(
+          std::ceil(1.0e6 / std::max(t_q_ns, 1.0)));
+
+      igl::parallel_for(O.rows(), [&](const int o)
+      {
+        W(o) = winding_number(intersector, O.row(o));
+      }, min_parallel);
+    }
+
+    /// Reference direction `x0` used to evaluate this scene.
+    const Direction & x0() const { return m_x0; }
+    /// Number of weighted boundary edge segments.
+    size_t num_boundary_segments() const { return m_boundary.size(); }
+    /// Number of triangles in the original mesh.
+    size_t num_faces() const { return m_face_count; }
+
+    /// Default reference direction: normalize(1, sqrt(2), sqrt(3)). A fixed
+    /// non-axis-aligned unit vector; any unit vector works per the paper;
+    /// this choice avoids accidental alignment with axis-aligned geometry.
+    static Direction default_x0()
+    {
+      Direction d(Scalar(1),
+                  Scalar(std::sqrt(2.0)),
+                  Scalar(std::sqrt(3.0)));
+      return d / d.norm();
+    }
+
+    /// Half the signed solid angle subtended by the spherical triangle
+    /// (x1, v0, v1) at the origin, via the unnormalized
+    /// Van Oosterom-Strackee formula. `x1` is expected to be a unit vector
+    /// (the antipodal "south pole" `-x0`); `v0`, `v1` need not be normalized.
+    /// Callers accumulate per-edge contributions and divide by 2π.
+    static Scalar half_solid_angle_unorm(
+      const Direction & x1, const Point & v0, const Point & v1)
+    {
+      const Scalar l0 = v0.norm();
+      const Scalar l1 = v1.norm();
+      const Scalar num   = x1.dot(v0.cross(v1));
+      const Scalar denom = l0 * l1
+                         + l1 * x1.dot(v0)
+                         + l0 * x1.dot(v1)
+                         + v0.dot(v1);
+      return std::atan2(num, denom);
+    }
+
+    /// Extract the open boundary as oriented, weighted edge segments.
+    ///
+    /// Each undirected edge {min, max} accumulates +1 for every triangle
+    /// that traverses it as min→max and -1 for every traversal max→min.
+    /// Interior edges of an oriented manifold cancel to zero and are
+    /// dropped. Surviving edges are emitted with positive weight; the
+    /// segment direction is flipped when the net count is negative so the
+    /// weight is always > 0. Non-manifold (≥3 incident triangles per edge)
+    /// is handled the same way. The surviving net count becomes the weight.
+    template <typename DerivedV, typename DerivedF>
+    static void build_boundary_segments(
+      const Eigen::MatrixBase<DerivedV> & V,
+      const Eigen::MatrixBase<DerivedF> & F,
+      std::vector<WeightedSeg> & out)
+    {
+      out.clear();
+      std::unordered_map<std::uint64_t, int> counts;
+      counts.reserve(static_cast<size_t>(F.rows()) * 3);
+
+      const auto pack = [](int lo, int hi) -> std::uint64_t
+      {
+        return (static_cast<std::uint64_t>(static_cast<std::uint32_t>(lo)) << 32)
+             |  static_cast<std::uint64_t>(static_cast<std::uint32_t>(hi));
+      };
+
+      for (Eigen::Index t = 0; t < F.rows(); ++t)
+      {
+        const int tri[3] = {
+          static_cast<int>(F(t, 0)),
+          static_cast<int>(F(t, 1)),
+          static_cast<int>(F(t, 2))
+        };
+        for (int e = 0; e < 3; ++e)
+        {
+          const int u = tri[e];
+          const int v = tri[(e + 1) % 3];
+          const int lo = u < v ? u : v;
+          const int hi = u < v ? v : u;
+          const std::uint64_t k = pack(lo, hi);
+          counts[k] += (u < v) ? +1 : -1;
+        }
+      }
+
+      out.reserve(counts.size());
+      for (const auto & kv : counts)
+      {
+        const int sum = kv.second;
+        if (sum == 0) continue;
+        const std::uint32_t lo = static_cast<std::uint32_t>(kv.first >> 32);
+        const std::uint32_t hi = static_cast<std::uint32_t>(kv.first & 0xFFFFFFFFu);
+        const Point a(static_cast<Scalar>(V(lo, 0)),
+                      static_cast<Scalar>(V(lo, 1)),
+                      static_cast<Scalar>(V(lo, 2)));
+        const Point b(static_cast<Scalar>(V(hi, 0)),
+                      static_cast<Scalar>(V(hi, 1)),
+                      static_cast<Scalar>(V(hi, 2)));
+        WeightedSeg ws;
+        if (sum > 0) { ws.a = a; ws.b = b; ws.w = static_cast<Scalar>(sum); }
+        else         { ws.a = b; ws.b = a; ws.w = static_cast<Scalar>(-sum); }
+        out.push_back(ws);
+      }
+    }
+
+    std::vector<WeightedSeg> m_boundary;
+    Direction m_x0;
+    size_t m_face_count = 0;
+  };
+}
+
+#endif

include/igl/WindingNumberTree.h

@@ -11,6 +11,7 @@
 #include <map>
 #include <Eigen/Dense>
 #include "WindingNumberMethod.h"
+#include <cassert>
 #include <memory>
 
 namespace igl

include/igl/active_set.cpp

@@ -14,6 +14,7 @@
 #include "placeholders.h"
 #include "PlainMatrix.h"
 
+#include <cassert>
 #include <iostream>
 #include <limits>
 #include <algorithm>

include/igl/adjacency_matrix.cpp

@@ -9,6 +9,7 @@
 
 #include "verbose.h"
 
+#include <cassert>
 #include <vector>
 
 template <typename DerivedF, typename T>

include/igl/all_pairs_distances.cpp

@@ -7,6 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "all_pairs_distances.h"
 #include <Eigen/Dense>
+#include <cassert>
 
 template <typename Mat>
 IGL_INLINE void igl::all_pairs_distances(

include/igl/angular_distance.cpp

@@ -8,6 +8,9 @@
 #include "angular_distance.h"
 #include "EPS.h"
 #include "PI.h"
+
+#include <cassert>
+
 IGL_INLINE double igl::angular_distance(
   const Eigen::Quaterniond & A,
   const Eigen::Quaterniond & B)

include/igl/arap_linear_block.cpp

@@ -9,6 +9,7 @@
 #include "verbose.h"
 #include "cotmatrix_entries.h"
 #include <Eigen/Dense>
+#include <cassert>
 
 template <typename MatV, typename MatF, typename MatK>
 IGL_INLINE void igl::arap_linear_block(

include/igl/arap_rhs.cpp

@@ -10,6 +10,7 @@
 #include "verbose.h"
 #include "repdiag.h"
 #include "cat.h"
+#include <cassert>
 #include <iostream>
 
 template<typename DerivedV, typename DerivedF, typename DerivedK>

include/igl/average_from_edges_onto_vertices.cpp

@@ -7,6 +7,8 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "average_from_edges_onto_vertices.h"
 
+#include <cassert>
+
 template<typename DerivedF,typename DerivedE,typename DerivedoE,
 typename DeriveduE,typename DeriveduV>
 IGL_INLINE void

include/igl/barycentric_coordinates.cpp

@@ -8,6 +8,8 @@
 #include "barycentric_coordinates.h"
 #include "volume.h"
 
+#include <cassert>
+
 template <
   typename DerivedP,
   typename DerivedA,

include/igl/biharmonic_coordinates.cpp

@@ -15,6 +15,7 @@
 #include "normal_derivative.h"
 #include "on_boundary.h"
 #include <Eigen/Sparse>
+#include <cassert>
 
 template <
   typename DerivedV,

include/igl/bijective_composite_harmonic_mapping.cpp

@@ -11,6 +11,7 @@
 #include "harmonic.h"
 #include "placeholders.h"
 //#include "matlab/MatlabWorkspace.h"
+#include <cassert>
 #include <iostream>
 
 template <

include/igl/blue_noise.cpp

@@ -17,6 +17,7 @@
 #include <vector>
 #include <random>
 #include <cstdint>
+#include <cassert>
 
 namespace igl
 {

include/igl/boundary_conditions.cpp

@@ -11,6 +11,7 @@
 #include "EPS.h"
 #include "project_to_line.h"
 
+#include <cassert>
 #include <vector>
 #include <map>
 #include <iostream>

include/igl/boundary_facets.cpp

@@ -15,6 +15,7 @@
 
 #include <Eigen/Core>
 
+#include <cassert>
 #include <map>
 #include <iostream>
 

include/igl/bounding_box.cpp

@@ -6,6 +6,7 @@
 // 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 "bounding_box.h"
+#include <cassert>
 #include <iostream>
 
 template <typename DerivedV, typename DerivedBV, typename DerivedBF>

include/igl/box_simplices.cpp

@@ -6,6 +6,7 @@
 // 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 "box_simplices.h"
+#include "parallel_for.h"
 
 template <
   typename DerivedV,
@@ -20,7 +21,8 @@ IGL_INLINE void igl::box_simplices(
 {
   B1.setConstant(F.rows(),V.cols(),std::numeric_limits<double>::infinity());
   B2.setConstant(F.rows(),V.cols(),-std::numeric_limits<double>::infinity());
-  for(int f = 0; f < F.rows(); f++) 
+  //for(int f = 0; f < F.rows(); f++) 
+  igl::parallel_for(F.rows(),[&](const int f)
   {
     for(int c = 0; c < F.cols(); c++) 
     {
@@ -30,7 +32,7 @@ IGL_INLINE void igl::box_simplices(
         B2(f,d) = std::max(B2(f,d),V(F(f,c),d));
       }
     }
-  }
+  },1000);
 }
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/cat.cpp

@@ -7,6 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "cat.h"
 
+#include <cassert>
 #include <cstdio>
 
 // Bug in unsupported/Eigen/SparseExtra needs iostream first

include/igl/centroid.cpp

@@ -8,6 +8,8 @@
 #include "centroid.h"
 #include <Eigen/Geometry>
 
+#include <cassert>
+
 template <
   typename DerivedV, 
   typename DerivedF, 

include/igl/circulation.h

@@ -47,7 +47,7 @@ namespace igl
   ///     F(f,:) opposite the vth corner, where EI(e,0)=v. Similarly EF(e,1) "
   ///     e=(j->i)
   /// @param[in] EI  #E by 2 list of edge flap corners (see above).
-  /// @param[out] #vN list of of faces touched by circulation (in cyclically order).
+  /// @param[out] #vN list of faces touched by circulation (in cyclically order).
   ///   
   /// \see edge_flaps
   template <typename DerivedEMAP, typename DerivedEF, typename DerivedEI, typename DerivedvN>

include/igl/circumradius.cpp

@@ -10,6 +10,8 @@
 #include "doublearea.h"
 #include "placeholders.h"
 #include <Eigen/QR>
+#include <cassert>
+
 template <
   typename DerivedV, 
   typename DerivedF,

include/igl/collapse_edge.cpp

@@ -8,6 +8,7 @@
 #include "collapse_edge.h"
 #include "circulation.h"
 #include "edge_collapse_is_valid.h"
+#include <cassert>
 #include <vector>
 
 template <

include/igl/collapse_edge.h

@@ -36,10 +36,10 @@ namespace igl
   ///     e=(j->i)
   /// @param[in,out] EI  #E by 2 list of edge flap corners (see above).
   /// [mesh inputs]
-  /// @param[out] e1  index into E of edge collpased on left
-  /// @param[out] e2  index into E of edge collpased on right
-  /// @param[out] f1  index into F of face collpased on left
-  /// @param[out] f2  index into F of face collpased on right
+  /// @param[out] e1  index into E of edge collapsed on left
+  /// @param[out] e2  index into E of edge collapsed on right
+  /// @param[out] f1  index into F of face collapsed on left
+  /// @param[out] f2  index into F of face collapsed on right
   /// @return true if edge was collapsed
   ///
   ///

include/igl/collapse_least_cost_edge.cpp

@@ -9,6 +9,8 @@
 #include "collapse_edge.h"
 #include "circulation.h"
 
+#include <cassert>
+
 IGL_INLINE bool igl::collapse_least_cost_edge(
   const decimate_cost_and_placement_callback & cost_and_placement,
   const decimate_pre_collapse_callback       & pre_collapse,

include/igl/collapse_least_cost_edge.h

@@ -47,10 +47,10 @@ namespace igl
   /// @param[in] C  #E by dim list of stored placements
   /// @param[out] e  index into E of attempted collapsed edge. Set to -1 if Q is empty or
   ///               contains only infinite cost edges.
-  /// @param[out] e1  index into E of edge collpased on left.
-  /// @param[out] e2  index into E of edge collpased on right.
-  /// @param[out] f1  index into F of face collpased on left.
-  /// @param[out] f2  index into F of face collpased on right.
+  /// @param[out] e1  index into E of edge collapsed on left.
+  /// @param[out] e2  index into E of edge collapsed on right.
+  /// @param[out] f1  index into F of face collapsed on left.
+  /// @param[out] f2  index into F of face collapsed on right.
   ///
   /// \bug This function is not templated nicely and refactoring it and its
   /// dependencies to do so is non-trivial, see

include/igl/collapse_small_triangles.cpp

@@ -13,6 +13,8 @@
 #include "colon.h"
 #include "faces_first.h"
 
+#include <cassert>
+
 #include <limits>
 
 #include <iostream>

include/igl/colormap.cpp

@@ -1567,6 +1567,266 @@ static double parula_cm[256][3] = {
   { 0.9736, 0.9752, 0.0597 },
   { 0.9763, 0.9831, 0.0538 }
 };
+static double zoe_cm[256][3] = {
+  {1.000000,1.000000,1.000000},
+{0.944653,0.947049,0.948510},
+{0.988743,0.993787,0.996863},
+{0.933958,0.941146,0.945529},
+{0.977485,0.987574,0.993725},
+{0.923264,0.935244,0.942549},
+{0.966228,0.981361,0.990588},
+{0.912569,0.929341,0.939568},
+{0.954970,0.975148,0.987451},
+{0.901874,0.923439,0.936588},
+{0.943713,0.968934,0.984313},
+{0.891180,0.917536,0.933608},
+{0.932455,0.962721,0.981176},
+{0.880485,0.911634,0.930627},
+{0.921198,0.956508,0.978039},
+{0.869791,0.905732,0.927647},
+{0.909941,0.950295,0.974901},
+{0.859096,0.899829,0.924666},
+{0.898683,0.944082,0.971764},
+{0.848402,0.893927,0.921686},
+{0.887426,0.937869,0.968627},
+{0.837707,0.888024,0.918705},
+{0.876168,0.931656,0.965490},
+{0.827013,0.882122,0.915725},
+{0.864911,0.925443,0.962352},
+{0.816318,0.876219,0.912744},
+{0.853653,0.919229,0.959215},
+{0.805623,0.870317,0.909764},
+{0.842396,0.913016,0.956078},
+{0.794929,0.864414,0.906783},
+{0.831138,0.906803,0.952940},
+{0.784234,0.858512,0.903803},
+{0.819881,0.900590,0.949803},
+{0.773540,0.852609,0.900823},
+{0.808624,0.894377,0.946666},
+{0.762845,0.846707,0.897842},
+{0.797366,0.888164,0.943528},
+{0.752151,0.840804,0.894862},
+{0.786109,0.881951,0.940391},
+{0.741456,0.834902,0.891881},
+{0.774851,0.875738,0.937254},
+{0.730761,0.829000,0.888901},
+{0.763594,0.869524,0.934116},
+{0.720067,0.823097,0.885920},
+{0.752336,0.863311,0.930979},
+{0.709372,0.817195,0.882940},
+{0.741079,0.857098,0.927842},
+{0.698678,0.811292,0.879959},
+{0.729822,0.850885,0.924704},
+{0.687983,0.805390,0.876979},
+{0.718564,0.844672,0.921567},
+{0.677289,0.799487,0.873998},
+{0.707307,0.838459,0.918430},
+{0.666594,0.793585,0.871018},
+{0.696049,0.832246,0.915292},
+{0.655900,0.787682,0.868038},
+{0.684792,0.826033,0.912155},
+{0.645205,0.781780,0.865057},
+{0.673534,0.819820,0.909018},
+{0.634510,0.775877,0.862077},
+{0.662277,0.813606,0.905880},
+{0.623816,0.769975,0.859096},
+{0.651020,0.807393,0.902743},
+{0.613121,0.764072,0.856116},
+{0.639762,0.801180,0.899606},
+{0.602427,0.758170,0.853135},
+{0.628505,0.794967,0.896469},
+{0.591732,0.752267,0.850155},
+{0.617247,0.788754,0.893331},
+{0.581038,0.746365,0.847174},
+{0.605990,0.782541,0.890194},
+{0.570343,0.740463,0.844194},
+{0.594732,0.776328,0.887057},
+{0.559648,0.734560,0.841214},
+{0.583475,0.770115,0.883919},
+{0.548954,0.728658,0.838233},
+{0.572218,0.763901,0.880782},
+{0.538259,0.722755,0.835253},
+{0.560960,0.757688,0.877645},
+{0.527565,0.716853,0.832272},
+{0.549703,0.751475,0.874507},
+{0.516870,0.710950,0.829292},
+{0.538445,0.745262,0.871370},
+{0.506176,0.705048,0.826311},
+{0.527188,0.739049,0.868233},
+{0.495481,0.699145,0.823331},
+{0.515930,0.732836,0.865095},
+{0.484787,0.693243,0.820350},
+{0.504673,0.726623,0.861958},
+{0.474092,0.687340,0.817370},
+{0.493415,0.720410,0.858821},
+{0.463397,0.681438,0.814389},
+{0.482158,0.714197,0.855683},
+{0.452703,0.675535,0.811409},
+{0.470901,0.707983,0.852546},
+{0.442008,0.669633,0.808429},
+{0.459643,0.701770,0.849409},
+{0.431314,0.663731,0.805448},
+{0.448386,0.695557,0.846271},
+{0.420619,0.657828,0.802468},
+{0.437128,0.689344,0.843134},
+{0.409925,0.651926,0.799487},
+{0.425871,0.683131,0.839997},
+{0.399230,0.646023,0.796507},
+{0.414613,0.676918,0.836859},
+{0.388535,0.640121,0.793526},
+{0.403356,0.670705,0.833722},
+{0.377841,0.634218,0.790546},
+{0.392099,0.664492,0.830585},
+{0.367146,0.628316,0.787565},
+{0.380841,0.658278,0.827448},
+{0.356452,0.622413,0.784585},
+{0.369584,0.652065,0.824310},
+{0.345757,0.616511,0.781604},
+{0.358326,0.645852,0.821173},
+{0.335063,0.610608,0.778624},
+{0.347069,0.639639,0.818036},
+{0.324368,0.604706,0.775644},
+{0.335811,0.633426,0.814898},
+{0.313674,0.598803,0.772663},
+{0.324554,0.627213,0.811761},
+{0.302979,0.592901,0.769683},
+{0.313297,0.621000,0.808624},
+{0.292284,0.586999,0.766702},
+{0.302039,0.614787,0.805486},
+{0.281590,0.581096,0.763722},
+{0.290782,0.608573,0.802349},
+{0.270895,0.575194,0.760741},
+{1.000000,0.978342,0.887274},
+{0.950000,0.925509,0.837329},
+{1.000000,0.970099,0.875524},
+{0.950000,0.917678,0.826167},
+{1.000000,0.961856,0.863775},
+{0.950000,0.909847,0.815005},
+{1.000000,0.953612,0.852025},
+{0.950000,0.902016,0.803843},
+{1.000000,0.945369,0.840275},
+{0.950000,0.894185,0.792681},
+{1.000000,0.937126,0.828526},
+{0.950000,0.886354,0.781519},
+{1.000000,0.928883,0.816776},
+{0.950000,0.878523,0.770356},
+{1.000000,0.920640,0.805027},
+{0.950000,0.870692,0.759194},
+{1.000000,0.912397,0.793277},
+{0.950000,0.862861,0.748032},
+{1.000000,0.904154,0.781528},
+{0.950000,0.855030,0.736870},
+{1.000000,0.895910,0.769778},
+{0.950000,0.847199,0.725708},
+{1.000000,0.887667,0.758028},
+{0.950000,0.839368,0.714546},
+{1.000000,0.879424,0.746279},
+{0.950000,0.831537,0.703384},
+{1.000000,0.871181,0.734529},
+{0.950000,0.823706,0.692222},
+{1.000000,0.862938,0.722780},
+{0.950000,0.815875,0.681060},
+{1.000000,0.854695,0.711030},
+{0.950000,0.808044,0.669898},
+{1.000000,0.846452,0.699281},
+{0.950000,0.800213,0.658736},
+{1.000000,0.838208,0.687531},
+{0.950000,0.792382,0.647574},
+{1.000000,0.829965,0.675782},
+{0.950000,0.784551,0.636411},
+{1.000000,0.821722,0.664032},
+{0.950000,0.776720,0.625249},
+{1.000000,0.813479,0.652282},
+{0.950000,0.768889,0.614087},
+{1.000000,0.805236,0.640533},
+{0.950000,0.761058,0.602925},
+{1.000000,0.796993,0.628783},
+{0.950000,0.753227,0.591763},
+{1.000000,0.788749,0.617034},
+{0.950000,0.745396,0.580601},
+{1.000000,0.780506,0.605284},
+{0.950000,0.737566,0.569439},
+{1.000000,0.772263,0.593535},
+{0.950000,0.729735,0.558277},
+{1.000000,0.764020,0.581785},
+{0.950000,0.721904,0.547115},
+{1.000000,0.755777,0.570035},
+{0.950000,0.714073,0.535953},
+{1.000000,0.747534,0.558286},
+{0.950000,0.706242,0.524791},
+{1.000000,0.739291,0.546536},
+{0.950000,0.698411,0.513629},
+{1.000000,0.731047,0.534787},
+{0.950000,0.690580,0.502466},
+{1.000000,0.722804,0.523037},
+{0.950000,0.682749,0.491304},
+{1.000000,0.714561,0.511288},
+{0.950000,0.674918,0.480142},
+{1.000000,0.706318,0.499538},
+{0.950000,0.667087,0.468980},
+{1.000000,0.698075,0.487789},
+{0.950000,0.659256,0.457818},
+{1.000000,0.689832,0.476039},
+{0.950000,0.651425,0.446656},
+{1.000000,0.681589,0.464289},
+{0.950000,0.643594,0.435494},
+{1.000000,0.673345,0.452540},
+{0.950000,0.635763,0.424332},
+{1.000000,0.665102,0.440790},
+{0.950000,0.627932,0.413170},
+{1.000000,0.656859,0.429041},
+{0.950000,0.620101,0.402008},
+{1.000000,0.648616,0.417291},
+{0.950000,0.612270,0.390846},
+{1.000000,0.640373,0.405542},
+{0.950000,0.604439,0.379683},
+{1.000000,0.632130,0.393792},
+{0.950000,0.596608,0.368521},
+{1.000000,0.623887,0.382042},
+{0.950000,0.588777,0.357359},
+{1.000000,0.615643,0.370293},
+{0.950000,0.580946,0.346197},
+{1.000000,0.607400,0.358543},
+{0.950000,0.573115,0.335035},
+{1.000000,0.599157,0.346794},
+{0.950000,0.565284,0.323873},
+{1.000000,0.590914,0.335044},
+{0.950000,0.557453,0.312711},
+{1.000000,0.582671,0.323295},
+{0.950000,0.549622,0.301549},
+{1.000000,0.574428,0.311545},
+{0.950000,0.541791,0.290387},
+{1.000000,0.566184,0.299796},
+{0.950000,0.533960,0.279225},
+{1.000000,0.557941,0.288046},
+{0.950000,0.526129,0.268063},
+{1.000000,0.549698,0.276296},
+{0.950000,0.518298,0.256901},
+{1.000000,0.541455,0.264547},
+{0.950000,0.510467,0.245738},
+{1.000000,0.533212,0.252797},
+{0.950000,0.502636,0.234576},
+{1.000000,0.524969,0.241048},
+{0.950000,0.494805,0.223414},
+{1.000000,0.516726,0.229298},
+{0.950000,0.486974,0.212252},
+{1.000000,0.508482,0.217549},
+{0.950000,0.479143,0.201090},
+{1.000000,0.500239,0.205799},
+{0.950000,0.471312,0.189928},
+{1.000000,0.491996,0.194050},
+{0.950000,0.463481,0.178766},
+{1.000000,0.483753,0.182300},
+{0.950000,0.455650,0.167604},
+{1.000000,0.475510,0.170550},
+{0.950000,0.447819,0.156442},
+{1.000000,0.467267,0.158801},
+{0.950000,0.439988,0.145280},
+{1.000000,0.459024,0.147051},
+{0.950000,0.432157,0.134118}
+};
+
+
 }
 
 template <typename T>
@@ -1602,6 +1862,9 @@ IGL_INLINE void igl::colormap(
     case COLOR_MAP_TYPE_VIRIDIS:
       colormap(viridis_cm, x_in, r, g, b);
       break;
+    case COLOR_MAP_TYPE_ZOE:
+      colormap(zoe_cm, x_in, r, g, b);
+      break;
     default:
       throw std::invalid_argument("igl::colormap(): Selected colormap is unsupported!");
       break;

include/igl/colormap.h

@@ -24,7 +24,8 @@ namespace igl {
     COLOR_MAP_TYPE_PLASMA = 4,
     COLOR_MAP_TYPE_VIRIDIS = 5,
     COLOR_MAP_TYPE_TURBO = 6,
-    NUM_COLOR_MAP_TYPES = 7
+    COLOR_MAP_TYPE_ZOE = 7,
+    NUM_COLOR_MAP_TYPES = 8
   };
   /// Compute [r,g,b] values of the selected colormap for
   /// a given factor f between 0 and 1

include/igl/comb_cross_field.cpp

@@ -18,6 +18,8 @@
 
 #include "triangle_triangle_adjacency.h"
 
+#include <cassert>
+
 namespace igl {
   template <typename DerivedV, typename DerivedF>
   class Comb

include/igl/comb_frame_field.cpp

@@ -16,6 +16,8 @@
 #include "PI.h"
 #include "PlainMatrix.h"
 
+#include <cassert>
+
 template <typename DerivedV, typename DerivedF, typename DerivedP>
 IGL_INLINE void igl::comb_frame_field(
   const Eigen::MatrixBase<DerivedV> &V,

include/igl/connected_components.cpp

@@ -7,6 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 
 #include "connected_components.h"
+#include <cassert>
 #include <queue>
 
 template < typename Atype, typename DerivedC, typename DerivedK>

include/igl/copyleft/cgal/BinaryWindingNumberOperations.h

@@ -13,6 +13,7 @@
 #include "../../igl_inline.h"
 #include "../../MeshBooleanType.h"
 #include <Eigen/Core>
+#include <cassert>
 
 // TODO: This is not written according to libigl style. These should be
 // function handles.

include/igl/copyleft/cgal/extract_cells_single_component.cpp

@@ -15,6 +15,7 @@
 #include <set>
 #include <map>
 #include <queue>
+#include <cassert>
 
 //#define EXTRACT_CELLS_SINGLE_COMPONENT_TIMING
 

include/igl/copyleft/cgal/is_self_intersecting.cpp

@@ -7,6 +7,8 @@
 #include "remesh_self_intersections.h"
 #include "../../collapse_edge.h"
 
+#include <cassert>
+
 template <
   typename DerivedV,
   typename DerivedF>

include/igl/copyleft/cgal/mesh_boolean_type_to_funcs.cpp

@@ -1,6 +1,8 @@
 #include "mesh_boolean_type_to_funcs.h"
 #include "BinaryWindingNumberOperations.h"
 
+#include <cassert>
+
 IGL_INLINE void igl::copyleft::cgal::mesh_boolean_type_to_funcs(
   const MeshBooleanType & type,
   std::function<int(const Eigen::Matrix<int,1,Eigen::Dynamic>) >& wind_num_op,

include/igl/copyleft/cgal/oriented_bounding_box.cpp

@@ -1,4 +1,5 @@
 #include "oriented_bounding_box.h"
+#include <boost/range/has_range_iterator.hpp>
 #include <CGAL/Aff_transformation_3.h>
 //#include <CGAL/Simple_cartesian.h>
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

include/igl/copyleft/cgal/outer_edge.cpp

@@ -8,6 +8,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "outer_edge.h"
 #include "outer_vertex.h"
+#include <cassert>
 #include <iostream>
 #include <vector>
 

include/igl/copyleft/cgal/outer_vertex.cpp

@@ -8,6 +8,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "outer_vertex.h"
 #include "outer_edge.h"
+#include <cassert>
 #include <iostream>
 #include <vector>
 

include/igl/copyleft/cgal/peel_outer_hull_layers.cpp

@@ -10,6 +10,7 @@
 #include "../../LinSpaced.h"
 #include <vector>
 #include <iostream>
+#include <cassert>
 //#define IGL_PEEL_OUTER_HULL_LAYERS_DEBUG
 #ifdef IGL_PEEL_OUTER_HULL_LAYERS_DEBUG
 #include "../../writePLY.h"

include/igl/copyleft/cgal/trim_with_solid.h

@@ -29,11 +29,11 @@ namespace igl
         /// output and strictly fewer FLOPS than CHECK_EACH_FACE. There will be
         /// many tiny patches along the intersection of A and B. 
         CHECK_EACH_PATCH = 2,
-        /// Merge A and B into the same mesh and resolve all self-itnersections.
+        /// Merge A and B into the same mesh and resolve all self-intersections.
         /// Then "undo" remeshing on faces of A not involved in intersections
-        /// with B (i.e., self-intersections in A). Then seperate into patches
+        /// with B (i.e., self-intersections in A). Then separate into patches
         /// based on connected faces --- where connected means sharing a
-        /// _manifold edge_ --- then label each aptch as inside or outside.
+        /// _manifold edge_ --- then label each patch as inside or outside.
         /// Results in fewer patches than CHECK_EACH_PATCH but finding, meshing and
         /// mesh-undoing self-intersections in A can be costly. This
         /// could result in different output from CHECK_EACH_PATCH because of

include/igl/copyleft/cgal/wire_mesh.cpp

@@ -7,6 +7,7 @@
 #include "coplanar.h"
 #include "mesh_boolean.h"
 #include <Eigen/Geometry>
+#include <cassert>
 #include <vector>
 
 template <

include/igl/copyleft/marching_cubes.cpp

@@ -26,6 +26,7 @@
 #include "marching_cubes.h"
 #include "marching_cubes_tables.h"
 
+#include <cassert>
 #include <unordered_map>
 
 

include/igl/copyleft/tetgen/cdt.cpp

@@ -10,6 +10,8 @@
 #include "../../PlainMatrix.h"
 #include "tetrahedralize.h"
 
+#include <cassert>
+
 template <
   typename DerivedV, 
   typename DerivedF, 

include/igl/copyleft/tetgen/mesh_with_skeleton.cpp

@@ -11,6 +11,7 @@
 #include "../../sample_edges.h"
 #include "../../cat.h"
 
+#include <cassert>
 #include <iostream>
 // Default settings pq2Y tell tetgen to mesh interior of triangle mesh and
 // to produce a graded tet mesh

include/igl/copyleft/tetgen/tetgenio_to_tetmesh.cpp

@@ -11,6 +11,7 @@
 #include "../../list_to_matrix.h"
 
 // STL includes
+#include <cassert>
 #include <iostream>
 
 template <

include/igl/cotmatrix.cpp

@@ -15,6 +15,8 @@
 // Bug in unsupported/Eigen/SparseExtra needs iostream first
 #include <iostream>
 
+#include <cassert>
+
 template <typename DerivedV, typename DerivedF, typename Scalar>
 IGL_INLINE void igl::cotmatrix(
   const Eigen::MatrixBase<DerivedV> & V, 

include/igl/cotmatrix_entries.cpp

@@ -17,6 +17,7 @@
 
 #include "verbose.h"
 
+#include <cassert>
 
 template <typename DerivedV, typename DerivedF, typename DerivedC>
 IGL_INLINE void igl::cotmatrix_entries(

include/igl/covariance_scatter_matrix.cpp

@@ -13,6 +13,7 @@
 #include "verbose.h"
 #include "cat.h"
 #include "PI.h"
+#include <cassert>
 
 template <
   typename DerivedV, 

include/igl/cr_vector_curvature_correction.cpp

@@ -18,6 +18,7 @@
 
 #include "PI.h"
 
+#include <cassert>
 
 template <typename DerivedV, typename DerivedF, typename DerivedE,
 typename DerivedOE, typename ScalarK>

include/igl/cr_vector_laplacian.cpp

@@ -14,6 +14,7 @@
 #include "doublearea.h"
 #include "squared_edge_lengths.h"
 
+#include <cassert>
 
 template <typename DerivedV, typename DerivedF, typename DerivedE,
 typename DerivedOE, typename ScalarL>

include/igl/cr_vector_mass.cpp

@@ -7,6 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "cr_vector_mass.h"
 
+#include <cassert>
 #include <vector>
 
 #include "orient_halfedges.h"

include/igl/cross.cpp

@@ -6,6 +6,7 @@
 // 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 "cross.h"
+#include <cassert>
 
 // http://www.antisphere.com/Wiki/tools:anttweakbar
 IGL_INLINE void igl::cross(

include/igl/crouzeix_raviart_cotmatrix.cpp

@@ -10,6 +10,7 @@
 #include "oriented_facets.h"
 #include "is_edge_manifold.h"
 #include "cotmatrix_entries.h"
+#include <cassert>
 
 template <typename DerivedV, typename DerivedF, typename LT, typename DerivedE, typename DerivedEMAP>
 void igl::crouzeix_raviart_cotmatrix(

include/igl/cubic.cpp

@@ -0,0 +1,41 @@
+#include "cubic.h"
+#include <cassert>
+
+  template
+    <
+    typename  DerivedC,
+    typename  DerivedP>
+  IGL_INLINE void igl::cubic(
+    const Eigen::MatrixBase<DerivedC>& C,
+    const typename DerivedP::Scalar & t,
+    Eigen::PlainObjectBase<DerivedP>& P)
+{
+  // static assert that C has 4 rows or Dynamic
+  static_assert(
+    DerivedC::RowsAtCompileTime == 4 ||
+    DerivedC::RowsAtCompileTime == Eigen::Dynamic,
+    "C must have 4 rows.");
+  // runtime assert that C has 4 rows
+  assert(C.rows() == 4 && "C must have 4 rows.");
+  using Scalar = typename DerivedP::Scalar;
+  // Evaluate cubic Bezier at parameter t
+  P = 
+    (Scalar(1) - t) * (Scalar(1) - t) * (Scalar(1) - t) * C.row(0)
+    + Scalar(3) * (Scalar(1) - t) * (Scalar(1) - t) * t * C.row(1)
+    + Scalar(3) * (Scalar(1) - t) * t * t * C.row(2)
+    + t * t * t * C.row(3);
+}
+
+
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::cubic<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, 1, -1, 1, 1, -1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::Matrix<double, 1, -1, 1, 1, -1>::Scalar const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1>>&);
+// generated by autoexplicit.sh
+template void igl::cubic<Eigen::Matrix<double, 4, -1, 0, 4, -1>, Eigen::Matrix<double, 1, -1, 1, 1, -1>>(Eigen::MatrixBase<Eigen::Matrix<double, 4, -1, 0, 4, -1>> const&, Eigen::Matrix<double, 1, -1, 1, 1, -1>::Scalar const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1>>&);
+// generated by autoexplicit.sh
+template void igl::cubic<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, -1, 1, true>, Eigen::Matrix<double, 1, 1, 0, 1, 1>>(Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, -1, 1, true>> const&, Eigen::Matrix<double, 1, 1, 0, 1, 1>::Scalar const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 1, 0, 1, 1>>&);
+// generated by autoexplicit.sh
+template void igl::cubic<Eigen::IndexedView<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Block<Eigen::Matrix<int, -1, -1, 0, -1, -1>, 1, -1, false>, Eigen::internal::AllRange<-1>>, Eigen::Matrix<double, 1, -1, 1, 1, -1>>(Eigen::MatrixBase<Eigen::IndexedView<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Block<Eigen::Matrix<int, -1, -1, 0, -1, -1>, 1, -1, false>, Eigen::internal::AllRange<-1>>> const&, Eigen::Matrix<double, 1, -1, 1, 1, -1>::Scalar const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1>>&);
+#endif

include/igl/cubic.h

@@ -0,0 +1,29 @@
+#ifndef IGL_CUBIC_H
+#define IGL_CUBIC_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  /// Evaluate a cubic Bézier curve defined by control points C at parameter t.
+  ///
+  /// @param[in] C  4 by dimensions matrix of control points for a cubic
+  /// Bézier curve
+  /// @param[in] t  Parameter at which to evaluate curve
+  /// @param[out] P  1 by dimensions point on curve C(t)
+  template
+    <
+    typename  DerivedC,
+    typename  DerivedP>
+  IGL_INLINE void cubic(
+    const Eigen::MatrixBase<DerivedC>& C,
+    const typename DerivedP::Scalar & t,
+    Eigen::PlainObjectBase<DerivedP>& P);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "cubic.cpp"
+#endif
+
+#endif

include/igl/cubic_is_flat.cpp

@@ -0,0 +1,23 @@
+#include "cubic_is_flat.h"
+
+template <typename DerivedC>
+IGL_INLINE bool igl::cubic_is_flat(
+  const Eigen::MatrixBase<DerivedC>& C,
+  const typename DerivedC::Scalar squared_distance_bound)
+{
+  using Scalar = typename DerivedC::Scalar;
+  const auto u = (Scalar(3) * C.row(1) - Scalar(2) * C.row(0) - C.row(3)).array().square().eval();
+  const auto v = (Scalar(3) * C.row(2) - Scalar(2) * C.row(3) - C.row(0)).array().square().eval();
+  const auto max_uv = u.cwiseMax(v).eval();
+  const Scalar tao = max_uv.sum();
+  const Scalar tolerance = Scalar(16) * squared_distance_bound;
+  return tao <= tolerance;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template bool igl::cubic_is_flat<Eigen::Matrix<double, 4, 2, 1, 4, 2>>(Eigen::MatrixBase<Eigen::Matrix<double, 4, 2, 1, 4, 2>> const&, Eigen::Matrix<double, 4, 2, 1, 4, 2>::Scalar);
+// generated by autoexplicit.sh
+template bool igl::cubic_is_flat<Eigen::Matrix<double, 4, 2, 0, 4, 2>>(Eigen::MatrixBase<Eigen::Matrix<double, 4, 2, 0, 4, 2>> const&, Eigen::Matrix<double, 4, 2, 0, 4, 2>::Scalar);
+#endif

include/igl/cubic_is_flat.h

@@ -0,0 +1,29 @@
+#ifndef IGL_CUBIC_IS_FLAT_H
+#define IGL_CUBIC_IS_FLAT_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  /// "Piecewise Linear Approximation of Bézier Curves" [Fischer 2000]. The test
+  /// computes a scalar and checks if it is less than 16 times
+  /// `squared_distance_bound`. If so, the curve's maximum squared distance to
+  /// the chord from first to last point will be less than
+  /// `squared_distance_bound`. The test is well behaved for degenerate input
+  /// (e.g., all rows in C are the same point).
+  ///
+  /// @param[in] C  4 by dimensions matrix of control points for a cubic
+  /// Bézier curve
+  /// @param[in] squared_distance_bound  Squared distance tolerance
+  /// @return  True if the curve is flat within the given tolerance
+  template <typename DerivedC>
+  IGL_INLINE bool cubic_is_flat(
+    const Eigen::MatrixBase<DerivedC>& C,
+    const typename DerivedC::Scalar squared_distance_bound);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "cubic_is_flat.cpp"
+#endif
+
+#endif

include/igl/cubic_monomial_bases.cpp

@@ -0,0 +1,58 @@
+#include "cubic_monomial_bases.h"
+#include <cassert>
+
+template <
+  typename DerivedC,
+  typename DerivedM,
+  typename DerivedD,
+  typename DerivedB>
+IGL_INLINE void igl::cubic_monomial_bases(
+  const Eigen::MatrixBase<DerivedC>& C,
+  Eigen::PlainObjectBase<DerivedM>& M,
+  Eigen::PlainObjectBase<DerivedD>& D,
+  Eigen::PlainObjectBase<DerivedB>& B)
+{
+  // static assert that C has 4 rows or Dynamic
+  static_assert(
+    DerivedC::RowsAtCompileTime == 4 ||
+    DerivedC::RowsAtCompileTime == Eigen::Dynamic,
+    "C must have 4 rows.");
+  // runtime assert that C has 4 rows
+  assert(C.rows() == 4 && "C must have 4 rows.");
+  const auto C0 = C.row(0);
+  const auto C1 = C.row(1);
+  const auto C2 = C.row(2);
+  const auto C3 = C.row(3);
+  // monomial coefficients for C(t)
+  M.resize(4,C.cols());
+  M <<
+    C.row(0).eval(),
+    3 * (C1 - C0),
+    3 * (C0 - 2 * C1 + C2),
+    C3 - C0 + 3 * (C1 - C2);
+  // Monomial coefficients for dC/dt
+  D.resize(3,C.cols());
+  D << M.row(1),
+     2 * M.row(2),
+     3 * M.row(3);
+  // Should work for row and column vectors
+  B.resize(6);
+  B.setConstant(0);
+  for (int i = 0; i < 3; ++i)
+  {
+    for (int j = 0; j < 4; ++j)
+    {
+      B(i + j) += D.row(i).dot(M.row(j));
+    }
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::cubic_monomial_bases<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, 4, -1, 0, 4, -1>, Eigen::Matrix<double, 3, -1, 0, 3, -1>, Eigen::Matrix<double, 1, 6, 1, 1, 6>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 4, -1, 0, 4, -1>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, -1, 0, 3, -1>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 6, 1, 1, 6>>&);
+// generated by autoexplicit.sh
+template void igl::cubic_monomial_bases<Eigen::Matrix<double, 4, -1, 0, 4, -1>, Eigen::Matrix<double, 4, -1, 0, 4, -1>, Eigen::Matrix<double, 3, -1, 0, 3, -1>, Eigen::Matrix<double, 6, 1, 0, 6, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, 4, -1, 0, 4, -1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 4, -1, 0, 4, -1>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, -1, 0, 3, -1>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, 6, 1, 0, 6, 1>>&);
+// generated by autoexplicit.sh
+template void igl::cubic_monomial_bases<Eigen::Matrix<double, 4, -1, 0, 4, -1>, Eigen::Matrix<double, 4, -1, 0, 4, -1>, Eigen::Matrix<double, 3, -1, 0, 3, -1>, Eigen::Matrix<double, 1, 6, 1, 1, 6>>(Eigen::MatrixBase<Eigen::Matrix<double, 4, -1, 0, 4, -1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 4, -1, 0, 4, -1>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, -1, 0, 3, -1>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 6, 1, 1, 6>>&);
+#endif

include/igl/cubic_monomial_bases.h

@@ -0,0 +1,34 @@
+#ifndef IGL_CUBIC_MONOMIAL_BASES_H
+#define IGL_CUBIC_MONOMIAL_BASES_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  /// Compute monomial basis representations for cubic Bézier curves.
+  ///
+  /// @param[in]  C  4 by dimensions matrix of control points for a cubic
+  /// Bézier curve
+  /// @param[out] M  4 by dimensions matrix of monomial coefficients for C(t)
+  /// @param[out] D  3 by dimensions matrix of monomial coefficients for dC/dt
+  /// @param[out] B  6-vector of inner products of those basis functions for C(t)
+  ///
+  /// \see igl::cycodebase::point_cubic_squared_distance
+  template <
+    typename DerivedC,
+    typename DerivedM,
+    typename DerivedD,
+    typename DerivedB>
+  IGL_INLINE void cubic_monomial_bases(
+    const Eigen::MatrixBase<DerivedC>& C,
+    Eigen::PlainObjectBase<DerivedM>& M,
+    Eigen::PlainObjectBase<DerivedD>& D,
+    Eigen::PlainObjectBase<DerivedB>& B);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "cubic_monomial_bases.cpp"
+#endif
+
+#endif

include/igl/cubic_split.cpp

@@ -0,0 +1,65 @@
+#include "cubic_split.h"
+
+template <
+  typename DerivedC,
+  typename DerivedK>
+IGL_INLINE void igl::cubic_split(
+  const Eigen::MatrixBase<DerivedC>& C,
+  const typename DerivedC::Scalar & t,
+  Eigen::PlainObjectBase<DerivedK>& C01,
+  Eigen::PlainObjectBase<DerivedK>& C012,
+  Eigen::PlainObjectBase<DerivedK>& C0123,
+  Eigen::PlainObjectBase<DerivedK>&  C123,
+  Eigen::PlainObjectBase<DerivedK>&   C23)
+{
+  const auto C0 = C.row(0);
+  const auto C1 = C.row(1);
+  const auto C2 = C.row(2);
+  const auto C3 = C.row(3);
+
+             C01 =  (C1 - C0) * t + C0;
+  const auto C12 = ((C2 - C1) * t + C1).eval();
+             C23 =  (C3 - C2) * t + C2;
+
+  C012 = (C12 - C01) * t + C01;
+  C123 = (C23 - C12) * t + C12;
+
+  C0123 = (C123 - C012) * t + C012;
+}
+
+template <
+  typename DerivedC,
+  typename DerivedK>
+IGL_INLINE void igl::cubic_split(
+  const Eigen::MatrixBase<DerivedC>& C,
+  const typename DerivedC::Scalar & t,
+  Eigen::PlainObjectBase<DerivedK>& C1,
+  Eigen::PlainObjectBase<DerivedK>& C2)
+{
+  using Scalar = typename DerivedC::Scalar;
+  typedef Eigen::Matrix<Scalar,1,DerivedC::ColsAtCompileTime> RowVectorS;
+  RowVectorS C01,C012,C0123,C123,C23;
+  igl::cubic_split(C,t,C01,C012,C0123,C123,C23);
+  C1.resize(4,C.cols());
+  C1 << C.row(0),
+        C01,
+        C012,
+        C0123;
+  C2.resize(4,C.cols());
+  C2 << C0123,
+        C123,
+        C23,
+        C.row(3);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::cubic_split<Eigen::Matrix<double, 4, -1, 1, 4, -1>, Eigen::Matrix<double, 4, 2, 1, 4, 2>>(Eigen::MatrixBase<Eigen::Matrix<double, 4, -1, 1, 4, -1>> const&, Eigen::Matrix<double, 4, -1, 1, 4, -1>::Scalar const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 4, 2, 1, 4, 2>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, 4, 2, 1, 4, 2>>&);
+// generated by autoexplicit.sh
+template void igl::cubic_split<Eigen::Matrix<double, 4, 2, 1, 4, 2>, Eigen::Matrix<double, 4, 2, 1, 4, 2>>(Eigen::MatrixBase<Eigen::Matrix<double, 4, 2, 1, 4, 2>> const&, Eigen::Matrix<double, 4, 2, 1, 4, 2>::Scalar const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 4, 2, 1, 4, 2>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, 4, 2, 1, 4, 2>>&);
+// generated by autoexplicit.sh
+template void igl::cubic_split<Eigen::Matrix<double, 4, 2, 0, 4, 2>, Eigen::Matrix<double, 4, 2, 1, 4, 2>>(Eigen::MatrixBase<Eigen::Matrix<double, 4, 2, 0, 4, 2>> const&, Eigen::Matrix<double, 4, 2, 0, 4, 2>::Scalar const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 4, 2, 1, 4, 2>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, 4, 2, 1, 4, 2>>&);
+template void igl::cubic_split<Eigen::Matrix<double, 4, 2, 0, 4, 2>, Eigen::Matrix<double, 4, 2, 0, 4, 2>>(Eigen::MatrixBase<Eigen::Matrix<double, 4, 2, 0, 4, 2>> const&, Eigen::Matrix<double, 4, 2, 0, 4, 2>::Scalar const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 4, 2, 0, 4, 2>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, 4, 2, 0, 4, 2>>&);
+template void igl::cubic_split<Eigen::Matrix<double, 4, 2, 0, 4, 2>, Eigen::Matrix<double, 1, 2, 1, 1, 2>>(Eigen::MatrixBase<Eigen::Matrix<double, 4, 2, 0, 4, 2>> const&, Eigen::Matrix<double, 4, 2, 0, 4, 2>::Scalar const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 2, 1, 1, 2>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 2, 1, 1, 2>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 2, 1, 1, 2>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 2, 1, 1, 2>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 2, 1, 1, 2>>&);
+#endif

include/igl/cubic_split.h

@@ -0,0 +1,52 @@
+#ifndef IGL_CUBIC_SPLIT_H
+#define IGL_CUBIC_SPLIT_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  /// Splits a cubic Bézier curve defined by control points C at parameter t
+  /// into two cubic Bézier curves defined by control points C1 and C2.
+  /// C1 = [C0,C01,C012,C0123], C2 = [C0123,C123,C23,C3]
+  ///
+  /// @param[in] C  4 by dimensions matrix of control points for a cubic
+  /// Bézier curve
+  /// @param[in] t  Parameter at which to split curve
+  /// @param[out] C01  1 by dim new control point
+  /// @param[out] C012  1 by dim new control point
+  /// @param[out] C0123  1 by dim new control point
+  /// @param[out] C123  1 by dim new control point
+  /// @param[out] C23  1 by dim new control point
+  ///
+  template <
+    typename DerivedC,
+    typename DerivedK>
+  IGL_INLINE void cubic_split(
+    const Eigen::MatrixBase<DerivedC>& C,
+    const typename DerivedC::Scalar & t,
+    Eigen::PlainObjectBase<DerivedK>& C01,
+    Eigen::PlainObjectBase<DerivedK>& C012,
+    Eigen::PlainObjectBase<DerivedK>& C0123,
+    Eigen::PlainObjectBase<DerivedK>&  C123,
+    Eigen::PlainObjectBase<DerivedK>&   C23);
+  ///
+  /// @param[out] C1  4 by dimensions matrix of control points for first cubic
+  /// Bézier curve from C(0) to C(t)
+  /// @param[out] C2  4 by dimensions matrix of control points for second cubic
+  /// Bézier curve from C(t) to C(1)
+  template <
+    typename DerivedC,
+    typename DerivedK>
+  IGL_INLINE void cubic_split(
+    const Eigen::MatrixBase<DerivedC>& C,
+    const typename DerivedC::Scalar & t,
+    Eigen::PlainObjectBase<DerivedK>& C1,
+    Eigen::PlainObjectBase<DerivedK>& C2);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "cubic_split.cpp"
+#endif
+
+#endif
+

include/igl/cycodebase/box_cubic.cpp

@@ -0,0 +1,83 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2026 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 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "box_cubic.h"
+#include "../cubic.h"
+#include "../parallel_for.h"
+#include "../placeholders.h"
+#include <cyPolynomial.h>
+
+template < 
+  typename DerivedC, 
+  typename DerivedB>
+IGL_INLINE void igl::cycodebase::box_cubic(
+  const Eigen::MatrixBase<DerivedC>& C,
+  Eigen::PlainObjectBase<DerivedB>& B1,
+  Eigen::PlainObjectBase<DerivedB>& B2)
+{
+  using Scalar = typename DerivedC::Scalar;
+  typedef Eigen::Matrix<Scalar,3,1> VectorS3;
+  typedef Eigen::Matrix<Scalar,DerivedC::RowsAtCompileTime,1> VectorSC;
+
+  // Using the control points is a simple (but not tight) bound
+  B1 = C({0,3},igl::placeholders::all).colwise().minCoeff();
+  B2 = C({0,3},igl::placeholders::all).colwise().maxCoeff();
+
+  // Better to find each of the t values where dC/dt = 0 and evaluate C there
+  for(int d = 0;d<C.cols();d++)
+  {
+    // [3*C1 - 3*C0, 6*C0 - 12*C1 + 6*C2, 9*C1 - 3*C0 - 9*C2 + 3*C3]
+    VectorS3 coef(
+      Scalar(3) * (C(1,d) - C(0,d)),
+      Scalar(6) * (C(0,d) - Scalar(2) * C(1,d) + C(2,d)),
+      Scalar(3) * (C(3,d) - C(0,d) + Scalar(3) * (C(1,d) - C(2,d)))
+    );
+    VectorS3 r;
+    int nr = cy::PolynomialRoots<2>(r.data(), coef.data(),0.0,1.0);
+    for(int i = 0;i<nr;i++)
+    {
+      const Scalar t = r[i];
+      // Evaluate cubic at t
+      Eigen::Matrix<Scalar,1,1> Cdt;
+      igl::cubic(C.col(d),t,Cdt);
+      B1(d) = std::min(B1(d),Cdt(0));
+      B2(d) = std::max(B2(d),Cdt(0));
+    }
+  }
+}
+
+template < 
+  typename DerivedP, 
+  typename DerivedC, 
+  typename DerivedB>
+IGL_INLINE void igl::cycodebase::box_cubic(
+  const Eigen::MatrixBase<DerivedP>& P,
+  const Eigen::MatrixBase<DerivedC>& C,
+  Eigen::PlainObjectBase<DerivedB>& B1,
+  Eigen::PlainObjectBase<DerivedB>& B2)
+{
+  B1.resize(C.rows(),P.cols());
+  B2.resize(C.rows(),P.cols());
+  typedef Eigen::Matrix<typename DerivedP::Scalar,1,DerivedP::ColsAtCompileTime> RowVectorP;
+
+  igl::parallel_for(C.rows(),[&](const int c)
+  {
+    RowVectorP B1_c, B2_c;
+    // Eval copies, but is it really good to make a template for the non copy?
+    box_cubic(P(C.row(c),igl::placeholders::all).eval(),B1_c,B2_c);
+    B1.row(c) = B1_c;
+    B2.row(c) = B2_c;
+  },1000);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::cycodebase::box_cubic<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 1, -1, -1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1>>&);
+template void igl::cycodebase::box_cubic<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, 1, -1, 1, 1, -1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1>>&);
+template void igl::cycodebase::box_cubic<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&);
+#endif

include/igl/cycodebase/box_cubic.h

@@ -0,0 +1,46 @@
+#ifndef IGL_CYCODEBASE_BOX_CUBIC_H
+#define IGL_CYCODEBASE_BOX_CUBIC_H
+
+#include "../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl {
+  namespace cycodebase {
+    /// Compute the min/max box corners tightly containing a given cubic bezier
+    /// curve.
+    ///
+    /// @param[in] C  4 by dim matrix of control points defining the cubic bezier
+    /// curve
+    /// @param[out] B1  1 by dim min corner of the bounding box
+    /// @param[out] B2  1 by dim max corner of the bounding box
+    ///
+    /// \see igl::box_simplices
+    template < 
+      typename DerivedC, 
+      typename DerivedB>
+    IGL_INLINE void box_cubic(
+      const Eigen::MatrixBase<DerivedC>& C,
+      Eigen::PlainObjectBase<DerivedB>& B1,
+      Eigen::PlainObjectBase<DerivedB>& B2);
+    /// \brief overload
+    ///
+    /// @param[in] P  #P by dim matrix of control point locations
+    /// @param[in] C  #C by 4 matrix of indices into P defining the cubics
+    /// @param[out] B1  #C by dim matrix of min corners of the bounding boxes
+    /// @param[out] B2  #C by dim matrix of max corners of the bounding boxes
+    template < 
+      typename DerivedP, 
+      typename DerivedC, 
+      typename DerivedB>
+    IGL_INLINE void box_cubic(
+      const Eigen::MatrixBase<DerivedP>& P,
+      const Eigen::MatrixBase<DerivedC>& C,
+      Eigen::PlainObjectBase<DerivedB>& B1,
+      Eigen::PlainObjectBase<DerivedB>& B2);
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+  #include "box_cubic.cpp"
+#endif
+#endif

include/igl/cycodebase/point_cubic_squared_distance.cpp

@@ -0,0 +1,196 @@
+#include "point_cubic_squared_distance.h"
+#include "../parallel_for.h"
+#include "../cubic.h"
+#include "../cubic_monomial_bases.h"
+#include <cyPolynomial.h>
+#include <limits>
+
+template <
+  typename DerivedQ,
+  typename DerivedC,
+  typename DerivedsqrD,
+  typename DerivedS,
+  typename DerivedK>
+void igl::cycodebase::point_cubic_squared_distance(
+  const Eigen::MatrixBase<DerivedQ>& Q,
+  const Eigen::MatrixBase<DerivedC>& C,
+  Eigen::PlainObjectBase<DerivedsqrD>& sqrD,
+  Eigen::PlainObjectBase<DerivedS>& S,
+  Eigen::PlainObjectBase<DerivedK>& K)
+{
+  using Scalar = typename DerivedQ::Scalar;
+  constexpr int ColsAtCompileTime = DerivedQ::ColsAtCompileTime;
+  const int dim = C.cols();
+  assert(C.rows() == 4 && "C should be 4 x dim.");
+  assert(Q.cols() == dim && "Q and C should have the same number of columns.");
+  // min_t ‖ Q - C(t) ‖^2
+  //
+  // Necessary condition for a minimum:
+  //
+  // d/dt ‖ Q - C(t) ‖^2 = 0
+  // => d/dt (Q - C(t)) · (Q - C(t)) = 0
+  // => -2 (dC/dt) · (Q - C(t)) = 0
+  // => (dC/dt) · (C(t) - Q) = f(t) = 0
+  //
+  using RowVectorSD = Eigen::RowVector<Scalar,ColsAtCompileTime>;
+  using MatrixS4D = Eigen::Matrix<Scalar,4,ColsAtCompileTime>;
+  using MatrixS3D = Eigen::Matrix<Scalar,3,ColsAtCompileTime>;
+  MatrixS4D M;
+  MatrixS3D D;
+  Eigen::RowVector<Scalar,6> B;
+  cubic_monomial_bases(C, M, D, B);
+  return point_cubic_squared_distance( Q, C, D, B, sqrD, S, K);
+}
+
+template <
+  typename DerivedQ,
+  typename DerivedC,
+  typename DerivedD,
+  typename DerivedB,
+  typename DerivedsqrD,
+  typename DerivedS,
+  typename DerivedK>
+void igl::cycodebase::point_cubic_squared_distance(
+  const Eigen::MatrixBase<DerivedQ>& Q,
+  const Eigen::MatrixBase<DerivedC>& C,
+  const Eigen::MatrixBase<DerivedD>& D,
+  const Eigen::MatrixBase<DerivedB>& B,
+  Eigen::PlainObjectBase<DerivedsqrD>& sqrD,
+  Eigen::PlainObjectBase<DerivedS>& S,
+  Eigen::PlainObjectBase<DerivedK>& K)
+{
+  // static assert that C has 4 rows or Dynamic
+  static_assert(
+    DerivedC::RowsAtCompileTime == 4 ||
+    DerivedC::RowsAtCompileTime == Eigen::Dynamic,
+    "C must have 4 rows.");
+  // runtime assert that C has 4 rows
+  assert(C.rows() == 4 && "C must have 4 rows.");
+  using Scalar = typename DerivedQ::Scalar;
+  constexpr int ColsAtCompileTime = DerivedQ::ColsAtCompileTime;
+  const int dim = C.cols();
+  using RowVectorSD = Eigen::RowVector<Scalar,ColsAtCompileTime>;
+  sqrD.setConstant(Q.rows(),1,std::numeric_limits<Scalar>::infinity());
+  S.resize(Q.rows());
+  K.resize(Q.rows(),dim);
+  const int n = Q.rows();
+  igl::parallel_for(n,[&](const int i)
+  {
+    RowVectorSD k_i;
+    const RowVectorSD q_i = Q.row(i);
+    point_cubic_squared_distance( q_i, C, D, B, sqrD(i), S(i), k_i);
+    K.row(i) = k_i;
+  }
+  ,1000);
+}
+
+template <
+  typename Derivedq,
+  typename DerivedC,
+  typename DerivedD,
+  typename DerivedB,
+  typename Derivedk
+  >
+void igl::cycodebase::point_cubic_squared_distance(
+  const Eigen::MatrixBase<Derivedq>& q,
+  const Eigen::MatrixBase<DerivedC>& C,
+  const Eigen::MatrixBase<DerivedD>& D,
+  const Eigen::MatrixBase<DerivedB>& B,
+  typename Derivedq::Scalar& sqrD,
+  typename Derivedq::Scalar& s,
+  Eigen::PlainObjectBase<Derivedk>& k)
+{
+  using cy::PolynomialRoots;
+  using Scalar = typename Derivedq::Scalar;
+  const int dim = C.cols();
+  // static assert that C has 4 rows or Dynamic
+  static_assert(
+    DerivedC::RowsAtCompileTime == 4 ||
+    DerivedC::RowsAtCompileTime == Eigen::Dynamic,
+    "C must have 4 rows.");
+  assert(C.rows() == 4 && "C should be 4 x dim.");
+  assert(q.cols() == dim && "q and C should have the same number of columns.");
+  typedef Eigen::Matrix<Scalar,1,Derivedq::ColsAtCompileTime> RowVectorSD;
+  // Fill in coefficients using precomputed data and Q
+  // => (dC/dt) · (C(t) - Q) = f(t) = 0
+  //
+  // C(t) = (1-t)^3 C0 + 3(1-t)^2 t C1 + 3(1-t) t^2 C2 + t^3 C3
+  // C(t) = C0 + t⋅3(C1 - C0) + t^2⋅3(C0 - 2C1 + C2) + t^3⋅(C3 - C0 + 3(C1 - C2))
+  // dC/dt =       3(C1 - C0) +   t⋅6(C0 - 2C1 + C2) + t^2⋅3(C3 - C0 + 3(C1 - C2))
+  Eigen::RowVector<Scalar,6> coef = B;
+  for (int j = 0; j < 3; ++j)
+  {
+    coef(j) -= D.row(j).dot(q);
+  }
+
+  sqrD = std::numeric_limits<Scalar>::infinity();
+  // f is a quintic polynomial:
+  constexpr int N = 5;
+  Scalar r[N];
+  int nr = PolynomialRoots<N>(r, coef.data(),Scalar(0),Scalar(1));
+  for(int j = 0;j<nr+2;j++)
+  {
+    Scalar t;
+    if(j==nr)
+    {
+      t = Scalar(0);
+    }
+    else if(j==nr+1)
+    {
+      t = Scalar(1);
+    }else
+    {
+      t = r[j];
+    }
+    RowVectorSD Ct;
+    igl::cubic(C,t,Ct);
+    const Scalar sqrD_j = (Ct - q).squaredNorm();
+    if(sqrD_j < sqrD)
+    {
+      sqrD = sqrD_j;
+      s = t;
+      k = Ct;
+    }
+  }
+}
+
+template <
+  typename Derivedq,
+  typename DerivedC,
+  typename Derivedk
+  >
+void igl::cycodebase::point_cubic_squared_distance(
+  const Eigen::MatrixBase<Derivedq>& q,
+  const Eigen::MatrixBase<DerivedC>& C,
+  typename Derivedq::Scalar& sqrD,
+  typename Derivedq::Scalar& s,
+  Eigen::PlainObjectBase<Derivedk>& k)
+{
+  using Scalar = typename Derivedq::Scalar;
+  constexpr int ColsAtCompileTime = Derivedq::ColsAtCompileTime;
+  static_assert(
+    DerivedC::RowsAtCompileTime == 4 ||
+    DerivedC::RowsAtCompileTime == Eigen::Dynamic,
+    "C must have 4 rows.");
+  static_assert(
+    int(Derivedq::ColsAtCompileTime) == int(DerivedC::ColsAtCompileTime),
+    "q and C must have the same number of columns.");
+  const int dim = C.cols();
+  assert(C.rows() == 4 && "C should be 4 x dim.");
+  assert(q.cols() == dim && "q and C should have the same number of columns.");
+  using RowVectorSD = Eigen::RowVector<Scalar,ColsAtCompileTime>;
+  using MatrixS4D = Eigen::Matrix<Scalar,4,ColsAtCompileTime>;
+  using MatrixS3D = Eigen::Matrix<Scalar,3,ColsAtCompileTime>;
+  MatrixS4D M;
+  MatrixS3D D;
+  Eigen::RowVector<Scalar,6> B;
+  cubic_monomial_bases(C, M, D, B);
+  return point_cubic_squared_distance( q, C, D, B, sqrD, s, k);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::cycodebase::point_cubic_squared_distance<Eigen::Matrix<double, 1, -1, 1, 1, -1>, Eigen::Matrix<double, 4, -1, 0, 4, -1>, Eigen::Matrix<double, 3, -1, 0, 3, -1>, Eigen::Matrix<double, 6, 1, 0, 6, 1>, Eigen::Matrix<double, 1, -1, 1, 1, -1>>(Eigen::MatrixBase<Eigen::Matrix<double, 1, -1, 1, 1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<double, 4, -1, 0, 4, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<double, 3, -1, 0, 3, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<double, 6, 1, 0, 6, 1>> const&, Eigen::Matrix<double, 1, -1, 1, 1, -1>::Scalar&, Eigen::Matrix<double, 1, -1, 1, 1, -1>::Scalar&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1>>&);
+template void igl::cycodebase::point_cubic_squared_distance<Eigen::Matrix<double, 1, -1, 1, 1, -1>, Eigen::Matrix<double, 4, -1, 0, 4, -1>, Eigen::Matrix<double, 1, -1, 1, 1, -1>>(Eigen::MatrixBase<Eigen::Matrix<double, 1, -1, 1, 1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<double, 4, -1, 0, 4, -1>> const&, Eigen::Matrix<double, 1, -1, 1, 1, -1>::Scalar&, Eigen::Matrix<double, 1, -1, 1, 1, -1>::Scalar&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1>>&);
+template void igl::cycodebase::point_cubic_squared_distance<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&);
+#endif

include/igl/cycodebase/point_cubic_squared_distance.h

@@ -0,0 +1,94 @@
+#ifndef IGL_CYCODEBASE_POINT_CUBIC_SQUARED_DISTANCE_H
+#define IGL_CYCODEBASE_POINT_CUBIC_SQUARED_DISTANCE_H
+
+#include "../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl {
+  namespace cycodebase {
+    /// Compute the squared distance from a set of points to a cubic Bezier
+    /// curve.
+    ///
+    /// @param[in] Q  #Q by dim matrix of query points
+    /// @param[in] C  4 by dim matrix of control points defining the cubic
+    /// bezier curve
+    /// @param[out] sqrD  #Q vector of squared distances from each query point
+    /// to the cubic bezier curve
+    /// @param[out] S  #Q vector of parametric locations on the cubic bezier
+    /// curve of the closest point to each query point
+    /// @param[out] K  #Q by dim matrix of closest point locations
+    ///
+    /// \see igl::cubic_monomial_bases, igl::cycodebase::point_spline_squared_distance
+    template <
+      typename DerivedQ,
+      typename DerivedC,
+      typename DerivedsqrD,
+      typename DerivedS,
+      typename DerivedK>
+    void point_cubic_squared_distance(
+      const Eigen::MatrixBase<DerivedQ>& Q,
+      const Eigen::MatrixBase<DerivedC>& C,
+      Eigen::PlainObjectBase<DerivedsqrD>& sqrD,
+      Eigen::PlainObjectBase<DerivedS>& S,
+      Eigen::PlainObjectBase<DerivedK>& K);
+    /// \brief overload
+    ///
+    /// @param[in] D 3 by dim matrix of monomial coefficients for dC/dt
+    /// @param[in] B 6-vector of inner products of monomial coefficients of C
+    /// and dCdt 
+    template <
+      typename DerivedQ,
+      typename DerivedC,
+      typename DerivedD,
+      typename DerivedB,
+      typename DerivedsqrD,
+      typename DerivedS,
+      typename DerivedK>
+    void point_cubic_squared_distance(
+      const Eigen::MatrixBase<DerivedQ>& Q,
+      const Eigen::MatrixBase<DerivedC>& C,
+      const Eigen::MatrixBase<DerivedD>& D,
+      const Eigen::MatrixBase<DerivedB>& B,
+      Eigen::PlainObjectBase<DerivedsqrD>& sqrD,
+      Eigen::PlainObjectBase<DerivedS>& S,
+      Eigen::PlainObjectBase<DerivedK>& K);
+    /// \brief single point overload
+    ///
+    /// @param[in] q  dim vector of a query point
+    template <
+      typename Derivedq,
+      typename DerivedC,
+      typename DerivedD,
+      typename DerivedB,
+      typename Derivedk
+      >
+    void point_cubic_squared_distance(
+      const Eigen::MatrixBase<Derivedq>& q,
+      const Eigen::MatrixBase<DerivedC>& C,
+      const Eigen::MatrixBase<DerivedD>& D,
+      const Eigen::MatrixBase<DerivedB>& B,
+      typename Derivedq::Scalar& sqrD,
+      typename Derivedq::Scalar& s,
+      Eigen::PlainObjectBase<Derivedk>& k);
+    /// \brief single point overload
+    template <
+      typename Derivedq,
+      typename DerivedC,
+      typename Derivedk
+      >
+    void point_cubic_squared_distance(
+      const Eigen::MatrixBase<Derivedq>& q,
+      const Eigen::MatrixBase<DerivedC>& C,
+      typename Derivedq::Scalar& sqrD,
+      typename Derivedq::Scalar& s,
+      Eigen::PlainObjectBase<Derivedk>& k);
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+  #include "point_cubic_squared_distance.cpp"
+#endif
+
+#endif
+
+

include/igl/cycodebase/point_spline_squared_distance.cpp

@@ -0,0 +1,120 @@
+#include "point_spline_squared_distance.h"
+#include "point_cubic_squared_distance.h"
+#include "spline_eytzinger_aabb.h"
+#include "../eytzinger_aabb_sdf.h"
+#include "../parallel_for.h"
+#include "../cubic_monomial_bases.h"
+#include "../placeholders.h"
+
+template <
+  typename DerivedQ,
+  typename DerivedP,
+  typename DerivedC,
+  typename DerivedsqrD,
+  typename DerivedI,
+  typename DerivedS,
+  typename DerivedK>
+void igl::cycodebase::point_spline_squared_distance(
+  const Eigen::MatrixBase<DerivedQ>& Q,
+  const Eigen::MatrixBase<DerivedP>& P,
+  const Eigen::MatrixBase<DerivedC>& C,
+  Eigen::PlainObjectBase<DerivedsqrD>& sqrD,
+  Eigen::PlainObjectBase<DerivedI>& I,
+  Eigen::PlainObjectBase<DerivedS>& S,
+  Eigen::PlainObjectBase<DerivedK>& K)
+{
+  static_assert(
+    DerivedQ::ColsAtCompileTime == DerivedP::ColsAtCompileTime,
+    "Q and P must have the same number of columns.");
+  using Scalar = typename DerivedQ::Scalar;
+  Eigen::Matrix<Scalar,DerivedC::RowsAtCompileTime,DerivedP::ColsAtCompileTime,Eigen::RowMajor> 
+    B1,B2;
+  Eigen::VectorXi leaf;
+  spline_eytzinger_aabb(P, C, B1, B2,leaf);
+  return point_spline_squared_distance(
+    Q,P,C,B1,B2,leaf,sqrD,I,S,K);
+}
+
+template <
+  typename DerivedQ,
+  typename DerivedP,
+  typename DerivedC,
+  typename DerivedB,
+  typename Derivedleaf,
+  typename DerivedsqrD,
+  typename DerivedI,
+  typename DerivedS,
+  typename DerivedK>
+void igl::cycodebase::point_spline_squared_distance(
+  const Eigen::MatrixBase<DerivedQ>& Q,
+  const Eigen::MatrixBase<DerivedP>& P,
+  const Eigen::MatrixBase<DerivedC>& C,
+  const Eigen::MatrixBase<DerivedB>& B1,
+  const Eigen::MatrixBase<DerivedB>& B2,
+  const Eigen::MatrixBase<Derivedleaf>& leaf,
+  Eigen::PlainObjectBase<DerivedsqrD>& sqrD,
+  Eigen::PlainObjectBase<DerivedI>& I,
+  Eigen::PlainObjectBase<DerivedS>& S,
+  Eigen::PlainObjectBase<DerivedK>& K)
+{
+  using Scalar = typename DerivedQ::Scalar;
+  using RowVectorS = Eigen::Matrix<Scalar, 1, DerivedQ::ColsAtCompileTime>;
+  using MatrixS4D = Eigen::Matrix<Scalar,4,DerivedQ::ColsAtCompileTime>;
+  using MatrixS3D = Eigen::Matrix<Scalar,3,DerivedQ::ColsAtCompileTime>;
+
+  // Cache these.
+  std::vector<MatrixS4D> C_vec(C.rows());
+  std::vector<MatrixS3D> D_vec(C.rows());
+  std::vector<Eigen::Matrix<Scalar,6,1>> B_vec(C.rows());
+  {
+    MatrixS4D M_unused;
+    for(int j = 0; j < C.rows(); j++)
+    {
+      C_vec[j] = P(C.row(j), igl::placeholders::all);
+      cubic_monomial_bases( C_vec[j], M_unused, D_vec[j], B_vec[j]);
+    }
+  }
+
+  sqrD.setConstant(Q.rows(),std::numeric_limits<Scalar>::infinity());
+  S.resize(Q.rows());
+  I.resize(Q.rows());
+  K.resize(Q.rows(),Q.cols());
+
+  //for(int i = 0; i < Q.rows(); i++)
+  igl::parallel_for(Q.rows(), [&](const int i)
+  {
+    const Eigen::Matrix<Scalar,1,DerivedP::ColsAtCompileTime> q = Q.row(i);
+    const std::function<Scalar(const int)> primitive_i = [&](const int j)
+    {
+      Scalar sqrD_ij,s_ij;
+      RowVectorS k_ij;
+      point_cubic_squared_distance( q, 
+        C_vec[j], D_vec[j], B_vec[j],
+        sqrD_ij, s_ij, k_ij);
+      if(sqrD_ij < sqrD(i))
+      {
+        sqrD(i) = sqrD_ij;
+        S(i) = s_ij;
+        K.row(i) = k_ij;
+        I(i) = j;
+      }
+      return sqrD_ij;
+    };
+    // This is expected (un)signed **unsquared** _distance. So let it compute
+    // that in there. We will square it and acculumate our results in the lambda
+    // anyway. So dist_i is ignored.
+    Scalar dist_i;
+    igl::eytzinger_aabb_sdf<true>(q,primitive_i,B1,B2,leaf,dist_i);
+  },1000);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::cycodebase::point_spline_squared_distance<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&);
+// generated by autoexplicit.sh
+template void igl::cycodebase::point_spline_squared_distance<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 2, 1, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 2, 1, -1, 2>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 2, 1, -1, 2>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&);
+// generated by autoexplicit.sh
+template void igl::cycodebase::point_spline_squared_distance<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&);
+template void igl::cycodebase::point_spline_squared_distance<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&);
+#endif

include/igl/cycodebase/point_spline_squared_distance.h

@@ -0,0 +1,80 @@
+#ifndef IGL_CYCODEBASE_POINT_SPLINE_SQUARED_DISTANCE_H
+#define IGL_CYCODEBASE_POINT_SPLINE_SQUARED_DISTANCE_H
+
+#include "../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl {
+  namespace cycodebase {
+    /// Compute the squared distance from a set of points to a spline defined by
+    /// cubic Bezier curves.
+    ///
+    /// @param[in] Q  #Q by dim matrix of query points
+    /// @param[in] P  #P by dim matrix of spline control points
+    /// @param[in] C  #C by 4 matrix of indices into P defining the cubic bezier
+    /// curves making up the spline
+    /// @param[out] sqrD  #Q vector of squared distances from each query point
+    /// to the spline
+    /// @param[out] I  #Q vector of indices of the closest cubic bezier curve in
+    /// the spline to each query point
+    /// @param[out] S  #Q vector of parametric locations on the closest cubic
+    /// bezier curve of the closest point to each query point
+    /// @param[out] K  #Q by dim matrix of closest point locations
+    ///
+    /// \see igl::cycodebase::point_cubic_squared_distance,
+    /// igl::cycodebase::spline_eytzinger_aabb
+    template <
+      typename DerivedQ,
+      typename DerivedP,
+      typename DerivedC,
+      typename DerivedsqrD,
+      typename DerivedI,
+      typename DerivedS,
+      typename DerivedK>
+    void point_spline_squared_distance(
+      const Eigen::MatrixBase<DerivedQ>& Q,
+      const Eigen::MatrixBase<DerivedP>& P,
+      const Eigen::MatrixBase<DerivedC>& C,
+      Eigen::PlainObjectBase<DerivedsqrD>& sqrD,
+      Eigen::PlainObjectBase<DerivedI>& I,
+      Eigen::PlainObjectBase<DerivedS>& S,
+      Eigen::PlainObjectBase<DerivedK>& K);
+    /// \brief overload with AABB acceleration structure
+    ///
+    /// @param[in] B1  #B by dim matrix of AABB min box corners
+    /// @param[in] B2  #B by dim matrix of AABB max box corners
+    /// @param[in] leaf  #B by 1 matrix of AABB leaf node indices/flags
+    ///
+    /// \see igl::eytzinger_aabb
+    template <
+      typename DerivedQ,
+      typename DerivedP,
+      typename DerivedC,
+      typename DerivedB,
+      typename Derivedleaf,
+      typename DerivedsqrD,
+      typename DerivedI,
+      typename DerivedS,
+      typename DerivedK>
+    void point_spline_squared_distance(
+      const Eigen::MatrixBase<DerivedQ>& Q,
+      const Eigen::MatrixBase<DerivedP>& P,
+      const Eigen::MatrixBase<DerivedC>& C,
+      const Eigen::MatrixBase<DerivedB>& B1,
+      const Eigen::MatrixBase<DerivedB>& B2,
+      const Eigen::MatrixBase<Derivedleaf>& leaf,
+      Eigen::PlainObjectBase<DerivedsqrD>& sqrD,
+      Eigen::PlainObjectBase<DerivedI>& I,
+      Eigen::PlainObjectBase<DerivedS>& S,
+      Eigen::PlainObjectBase<DerivedK>& K);
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+  #include "point_spline_squared_distance.cpp"
+#endif
+
+#endif
+
+
+

include/igl/cycodebase/roots.cpp

@@ -0,0 +1,83 @@
+#include "roots.h"
+#include <cyPolynomial.h>
+
+namespace
+{
+  template <int MAX_DEG, typename Derivedcoef, typename DerivedR>
+  int roots_helper(
+      const int deg,
+      const Eigen::MatrixBase<Derivedcoef>& coef,
+      const typename Derivedcoef::Scalar xmin,
+      const typename Derivedcoef::Scalar xmax,
+      Eigen::PlainObjectBase<DerivedR>& R)
+  {
+    using cy::PolynomialRoots;
+    if constexpr (MAX_DEG == 0)
+    {
+      // Should be impossible to reach here from roots()
+      return 0;
+    }
+    else
+    {
+      if (deg == MAX_DEG)
+      {
+        using Scalar = typename Derivedcoef::Scalar;
+        Scalar r[MAX_DEG];
+        int nr = PolynomialRoots<MAX_DEG>(r, coef.derived().data(), xmin, xmax);
+        // will resize if needed. Safe on compile-time sized matrices
+        R.setConstant(MAX_DEG, std::numeric_limits<Scalar>::quiet_NaN());
+        for (int i = 0; i < nr; ++i)
+        {
+          R[i] = r[i];
+        }
+        return nr;
+      }
+      return roots_helper<MAX_DEG - 1>(deg, coef, xmin, xmax, R);
+    }
+  }
+
+}
+
+template <
+  typename Derivedcoef,
+  typename DerivedR>
+IGL_INLINE int igl::cycodebase::roots(
+    const Eigen::MatrixBase<Derivedcoef>& coef,
+    const typename Derivedcoef::Scalar xmin,
+    const typename Derivedcoef::Scalar xmax,
+    Eigen::PlainObjectBase<DerivedR>& R)
+{
+  using Scalar = typename Derivedcoef::Scalar;
+  // Static assert that DerivedR::Scalar is same as Derivedcoef::Scalar
+  static_assert(
+      std::is_same<typename DerivedR::Scalar, Scalar>::value,
+      "DerivedR::Scalar must be the same as Derivedcoef::Scalar");
+  // Check that Derivedcoef and DerivedR are vectors        
+  static_assert(
+      (Derivedcoef::ColsAtCompileTime == 1 ||
+       Derivedcoef::RowsAtCompileTime == 1) &&
+      (DerivedR::ColsAtCompileTime == 1 ||
+       DerivedR::RowsAtCompileTime == 1),
+      "Derivedcoef and DerivedR must be vectors");
+
+  constexpr int MAX_DEG = 16;
+  const int deg = coef.size() - 1;
+
+  if (deg < 1 || deg > MAX_DEG)
+  {
+    throw std::runtime_error("Polynomial degree out of range");
+  }
+
+  return ::roots_helper<MAX_DEG>( deg, coef, xmin, xmax, R);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+template IGL_INLINE int igl::cycodebase::roots<Eigen::RowVectorXd, Eigen::RowVectorXd>(
+    const Eigen::MatrixBase<Eigen::RowVectorXd>& coef,
+    const double xmin,
+    const double xmax,
+    Eigen::PlainObjectBase<Eigen::RowVectorXd>& R);
+template int igl::cycodebase::roots<Eigen::Matrix<double, 4, 1, 0, 4, 1>, Eigen::Matrix<double, 3, 1, 0, 3, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, 4, 1, 0, 4, 1>> const&, Eigen::Matrix<double, 4, 1, 0, 4, 1>::Scalar, Eigen::Matrix<double, 4, 1, 0, 4, 1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 1, 0, 3, 1>>&);
+template int igl::cycodebase::roots<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1>::Scalar, Eigen::Matrix<double, -1, 1, 0, -1, 1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>>&);
+#endif

include/igl/cycodebase/roots.h

@@ -0,0 +1,33 @@
+#ifndef IGL_CYCODEBASE_ROOTS_H
+#define IGL_CYCODEBASE_ROOTS_H
+
+#include "../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl {
+  namespace cycodebase {
+    /// Compute the real roots of a polynomial with given coefficients within a
+    /// given interval [xmin, xmax].
+    ///
+    /// @param[in] coef  Vector of polynomial coefficients in increasing order
+    /// of degree (i.e., coef[0] + coef[1]*x + coef[2]*x^2 + ... )
+    /// @param[in] xmin  Minimum x value of the interval
+    /// @param[in] xmax  Maximum x value of the interval
+    /// @param[out] R  Vector of real roots within [xmin, xmax]
+    /// @return  Number of real roots found within [xmin, xmax]
+    template <
+      typename Derivedcoef,
+      typename DerivedR>
+    IGL_INLINE int roots(
+        const Eigen::MatrixBase<Derivedcoef>& coef,
+        const typename Derivedcoef::Scalar xmin,
+        const typename Derivedcoef::Scalar xmax,
+        Eigen::PlainObjectBase<DerivedR>& R);
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+  #include "roots.cpp"
+#endif
+
+#endif

include/igl/cycodebase/spline_eytzinger_aabb.cpp

@@ -0,0 +1,32 @@
+#include "spline_eytzinger_aabb.h"
+#include "box_cubic.h"
+#include "../eytzinger_aabb.h"
+
+template <
+  typename DerivedP,
+  typename DerivedC,
+  typename DerivedB,
+  typename Derivedleaf>
+void igl::cycodebase::spline_eytzinger_aabb(
+  const Eigen::MatrixBase<DerivedP>& P,
+  const Eigen::MatrixBase<DerivedC>& C,
+  Eigen::PlainObjectBase<DerivedB>& B1,
+  Eigen::PlainObjectBase<DerivedB>& B2,
+  Eigen::PlainObjectBase<Derivedleaf>& leaf)
+{
+  using Scalar = typename DerivedP::Scalar;
+  Eigen::Matrix<Scalar,DerivedC::RowsAtCompileTime,DerivedP::ColsAtCompileTime,Eigen::RowMajor> 
+    PB1,PB2;
+  box_cubic(P,C,PB1,PB2);
+  eytzinger_aabb( PB1, PB2, B1, B2,leaf);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::cycodebase::spline_eytzinger_aabb<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>>&);
+// generated by autoexplicit.sh
+template void igl::cycodebase::spline_eytzinger_aabb<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 2, 1, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 1, -1, 2>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 1, -1, 2>>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>>&);
+// generated by autoexplicit.sh
+template void igl::cycodebase::spline_eytzinger_aabb<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1>>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>>&);
+#endif

include/igl/cycodebase/spline_eytzinger_aabb.h

@@ -0,0 +1,43 @@
+#ifndef IGL_CYCODEBASE_SPLINE_EYTZINGER_AABB_H
+#define IGL_CYCODEBASE_SPLINE_EYTZINGER_AABB_H
+
+#include "../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl {
+  namespace cycodebase {
+    /// Compute an AABB acceleration structure for a spline of cubic Bezier
+    /// curves in Eytzinger layout.
+    ///
+    /// @param[in] P  #P by dim matrix of spline control points
+    /// @param[in] C  #C by 4 matrix of indices into P defining the cubic Bézier
+    /// curves making up the spline
+    /// @param[out] B1  #B by dim matrix of AABB min box corners
+    /// @param[out] B2  #B by dim matrix of AABB max box corners
+    /// @param[out] leaf  #B by 1 matrix of AABB leaf node indices/flags
+    ///
+    /// \see igl::cycodebase::box_cubic, igl::eytzinger_aabb
+    ///
+    template <
+      typename DerivedP,
+      typename DerivedC,
+      typename DerivedB,
+      typename Derivedleaf>
+    void spline_eytzinger_aabb(
+      const Eigen::MatrixBase<DerivedP>& P,
+      const Eigen::MatrixBase<DerivedC>& C,
+      Eigen::PlainObjectBase<DerivedB>& B1,
+      Eigen::PlainObjectBase<DerivedB>& B2,
+      Eigen::PlainObjectBase<Derivedleaf>& leaf);
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+  #include "spline_eytzinger_aabb.cpp"
+#endif
+
+#endif
+
+
+
+

include/igl/decimate_callback_types.h

@@ -64,10 +64,10 @@ namespace igl
   /// @param[in] C  #E by dim list of stored placements
   /// @param[in] e  index into E of attempted collapsed edge. Set to -1 if Q is empty or
   ///               contains only infinite cost edges.
-  /// @param[in] e1  index into E of edge collpased on left.
-  /// @param[in] e2  index into E of edge collpased on right.
-  /// @param[in] f1  index into F of face collpased on left.
-  /// @param[in] f2  index into F of face collpased on right.
+  /// @param[in] e1  index into E of edge collapsed on left.
+  /// @param[in] e2  index into E of edge collapsed on right.
+  /// @param[in] f1  index into F of face collapsed on left.
+  /// @param[in] f2  index into F of face collapsed on right.
   /// @return whether to stop
   using decimate_stopping_condition_callback = 
     std::function<bool(
@@ -139,10 +139,10 @@ namespace igl
   /// @param[in] C  #E by dim list of stored placements
   /// @param[in] e  index into E of attempted collapsed edge. Set to -1 if Q is empty or
   ///               contains only infinite cost edges.
-  /// @param[in] e1  index into E of edge collpased on left.
-  /// @param[in] e2  index into E of edge collpased on right.
-  /// @param[in] f1  index into F of face collpased on left.
-  /// @param[in] f2  index into F of face collpased on right.
+  /// @param[in] e1  index into E of edge collapsed on left.
+  /// @param[in] e2  index into E of edge collapsed on right.
+  /// @param[in] f1  index into F of face collapsed on left.
+  /// @param[in] f2  index into F of face collapsed on right.
   /// @param[in] collapsed whether collapse actual took place
   using decimate_post_collapse_callback = 
     std::function<void(

include/igl/deform_skeleton.cpp

@@ -6,6 +6,8 @@
 // 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 "deform_skeleton.h"
+#include <cassert>
+
 void igl::deform_skeleton(
   const Eigen::MatrixXd & C,
   const Eigen::MatrixXi & BE,

include/igl/delaunay_triangulation.cpp

@@ -12,6 +12,7 @@
 #include "unique_edge_map.h"
 #include "is_delaunay.h"
 
+#include <cassert>
 #include <vector>
 #include <sstream>
 

include/igl/direct_delta_mush.cpp

@@ -9,6 +9,8 @@
 #include "cotmatrix.h"
 #include "PlainMatrix.h"
 
+#include <cassert>
+
 template <
   typename DerivedV,
   typename DerivedOmega,

include/igl/edge_collapse_is_valid.cpp

@@ -11,6 +11,7 @@
 #include "intersect.h"
 #include "unique.h"
 #include "list_to_matrix.h"
+#include <cassert>
 #include <vector>
 
 template <

include/igl/edge_exists_near.cpp

@@ -6,6 +6,7 @@
 // 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 "edge_exists_near.h"
+#include <cassert>
 
 template <
   typename DeriveduE,

include/igl/edge_midpoints.cpp

@@ -7,6 +7,8 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "edge_midpoints.h"
 
+#include <cassert>
+
 template<typename DerivedV,typename DerivedF,typename DerivedE,
 typename DerivedoE, typename Derivedmps>
 IGL_INLINE void

include/igl/edge_topology.cpp

@@ -8,6 +8,7 @@
 #include "edge_topology.h"
 #include "is_edge_manifold.h"
 #include <algorithm>
+#include <cassert>
 
 template<typename DerivedV, typename DerivedF, typename DerivedE>
 IGL_INLINE void igl::edge_topology(

include/igl/edge_vectors.cpp

@@ -13,6 +13,7 @@
 
 #include "PI.h"
 
+#include <cassert>
 
 template<typename DerivedV,typename DerivedF,typename DerivedE,
 typename DerivedoE, typename Derivedvec>

include/igl/edges.cpp

@@ -7,6 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "edges.h"
 #include "adjacency_matrix.h"
+#include <cassert>
 #include <iostream>
 
 template <typename DerivedF, typename DerivedE>

include/igl/edges_to_path.cpp

@@ -5,6 +5,7 @@
 #include "unique.h"
 #include "adjacency_list.h"
 #include "PlainMatrix.h"
+#include <cassert>
 
 template <
   typename DerivedE,

include/igl/eigs.cpp

@@ -12,6 +12,7 @@
 #include "massmatrix.h"
 #include "placeholders.h"
 #include <iostream>
+#include <cassert>
 
 template <
   typename Atype,

include/igl/embree/EmbreeDevice.h

@@ -12,7 +12,9 @@
 #ifndef IGL_EMBREE_EMBREE_DEVICE_H
 #define IGL_EMBREE_EMBREE_DEVICE_H
 #include <embree4/rtcore.h>
+#include <cassert>
 #include <iostream>
+#include <cassert>
 
 namespace igl
 {

include/igl/embree/EmbreeIntersector.cpp

@@ -386,3 +386,58 @@ igl::embree::EmbreeIntersector
   ray.hit.instID[0] = RTC_INVALID_GEOMETRY_ID;
   ray.hit.primID = RTC_INVALID_GEOMETRY_ID;
 }
+
+IGL_INLINE int
+igl::embree::EmbreeIntersector
+::signedIntersectionsRay(
+  const OriginType    & origin,
+  const DirectionType & direction,
+  float tnear,
+  float tfar,
+  int   mask) const
+{
+  struct query_context
+  {
+    RTCRayQueryContext base; // MUST be first for the reinterpret_cast in the filter
+    int sum;
+  };
+
+  query_context q;
+  rtcInitRayQueryContext(&q.base);
+  q.sum = 0;
+
+  RTCRay ray{};
+  ray.org_x = origin[0];
+  ray.org_y = origin[1];
+  ray.org_z = origin[2];
+  ray.dir_x = direction[0];
+  ray.dir_y = direction[1];
+  ray.dir_z = direction[2];
+  ray.tnear = tnear;
+  ray.tfar  = tfar;
+  ray.mask  = static_cast<unsigned int>(mask);
+  ray.flags = 0;
+
+  RTCOccludedArguments rargs;
+  rtcInitOccludedArguments(&rargs);
+  rargs.flags = (RTCRayQueryFlags)(
+      RTC_RAY_QUERY_FLAG_COHERENT | RTC_RAY_QUERY_FLAG_INVOKE_ARGUMENT_FILTER);
+  rargs.feature_mask = RTC_FEATURE_FLAG_ALL;
+  rargs.context = &q.base;
+  rargs.filter = +[](RTCFilterFunctionNArguments const* fargs)
+  {
+    assert(fargs->N == 1 && fargs->valid[0]);
+    auto const& fray = reinterpret_cast<RTCRay&>(*fargs->ray);
+    auto const& fhit = reinterpret_cast<RTCHit&>(*fargs->hit);
+    auto const d = fray.dir_x * fhit.Ng_x
+                 + fray.dir_y * fhit.Ng_y
+                 + fray.dir_z * fhit.Ng_z;
+    reinterpret_cast<query_context*>(fargs->context)->sum += d > 0.0f ? +1 : -1;
+    // Reject the hit so traversal continues and we visit every crossing.
+    fargs->valid[0] = 0;
+  };
+  rargs.occluded = nullptr;
+
+  rtcOccluded1(scene, &ray, &rargs);
+  return q.sum;
+}

include/igl/embree/EmbreeIntersector.h

@@ -37,6 +37,8 @@ namespace igl
     public:
       typedef Eigen::Matrix<float,Eigen::Dynamic,3> PointMatrixType;
       typedef Eigen::Matrix<int,Eigen::Dynamic,3> FaceMatrixType;
+      typedef Eigen::RowVector3f OriginType;
+      typedef Eigen::RowVector3f DirectionType;
     public:
       EmbreeIntersector();
     private:
@@ -150,6 +152,23 @@ namespace igl
         Hit<float> &hit,
         int mask = 0xFFFFFFFF) const;
 
+      /// Signed ray-mesh crossing count along `(origin, direction)`. Visits
+      /// every hit (occluded ray + argument filter rejecting each) and
+      /// accumulates `sign(direction · Ng)` per hit.
+      ///
+      /// @param[in] origin     ray origin
+      /// @param[in] direction  ray direction (need not be normalized)
+      /// @param[in] tnear      start of ray segment
+      /// @param[in] tfar       end of ray segment
+      /// @param[in] mask       a 32 bit mask to identify active geometries
+      /// @return signed crossing count
+      int signedIntersectionsRay(
+        const OriginType   & origin,
+        const DirectionType& direction,
+        float tnear = 0.0f,
+        float tfar  = std::numeric_limits<float>::infinity(),
+        int   mask  = 0xFFFFFFFF) const;
+
     private:
 
       struct Vertex   {float x,y,z,a;};