better cm

doc
working spline winding number and demo
2026-01-26 19:43:58 -05:00 · 2026-01-26 18:21:23 -05:00 · 2026-01-26 17:11:02 -05:00 · 2026-01-26 12:13:43 -05:00 · 2026-01-24 23:27:28 -05:00
67 changed files with 2945 additions and 119 deletions
@@ -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})
@@ -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})
+
@@ -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)
@@ -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)
@@ -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)

@@ -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
+  };
 }


@@ -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
@@ -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;
@@ -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
@@ -0,0 +1,39 @@
+#include "cubic.h"
+  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
@@ -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
@@ -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
@@ -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
@@ -0,0 +1,57 @@
+#include "cubic_monomial_bases.h"
+
+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
@@ -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
@@ -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
@@ -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
+
@@ -0,0 +1,82 @@
+// 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 <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},Eigen::all).colwise().minCoeff();
+  B2 = C({0,3},Eigen::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),Eigen::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
@@ -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
@@ -0,0 +1,197 @@
+#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();
+  for(int i = 0;i<n;i++)
+  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
@@ -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
+
+
@@ -0,0 +1,117 @@
+#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"
+
+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), Eigen::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, 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
@@ -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
+
+
+
@@ -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
@@ -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
@@ -0,0 +1,30 @@
+#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, 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
@@ -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
+
+
+
+
@@ -85,6 +85,12 @@ IGL_INLINE void igl::eytzinger_aabb(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::eytzinger_aabb<Eigen::Matrix<double, -1, -1, 1, -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, 1, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -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::eytzinger_aabb<Eigen::Matrix<double, -1, -1, 1, -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, 1, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -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>>&);
+// generated by autoexplicit.sh
+template void igl::eytzinger_aabb<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::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<int, -1, 1, 0, -1, 1>>&);
+// generated by autoexplicit.sh
 template void igl::eytzinger_aabb<Eigen::Matrix<double, -1, 2, 1, -1, 2>, Eigen::Matrix<double, -1, 2, 1, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 2, 1, -1, 2>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 2, 1, -1, 2>> 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::eytzinger_aabb<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>>&);
@@ -2,6 +2,7 @@
 #include "parallel_for.h"

 template <
+  bool Squared,
  typename Derivedp,
  typename Func,
  typename DerivedB,
@@ -25,9 +26,16 @@ IGL_INLINE void igl::eytzinger_aabb_sdf(
    // relative position
    const RowVectorS r = p - 0.5*(b1 + b2);
    const RowVectorS q = r.cwiseAbs() - b;
-    const Scalar t1 = q.cwiseMax(0.0).squaredNorm();
-    const Scalar t2 = std::min(q.maxCoeff(),0.0);
-    return std::sqrt(t1 + t2);
+    if constexpr (Squared)
+    {
+      return q.cwiseMax(Scalar(0)).squaredNorm();
+    }else
+    {
+      const Scalar t1 = q.cwiseMax(0.0).norm();
+      const Scalar t2 = std::min(q.maxCoeff(),0.0);
+      const Scalar sdf = t1 + t2;
+      return sdf;
+    }
  };

  // List of current active nodes
@@ -63,6 +71,7 @@ IGL_INLINE void igl::eytzinger_aabb_sdf(
 }

 template <
+  bool Squared,
  typename DerivedP,
  typename Func,
  typename DerivedB,
@@ -87,16 +96,18 @@ IGL_INLINE void igl::eytzinger_aabb_sdf(
    {
      return primitive(p,j);
    };
-    igl::eytzinger_aabb_sdf(p,primitive_i,B1,B2,leaf,S(i));
+    igl::eytzinger_aabb_sdf<Squared>(p,primitive_i,B1,B2,leaf,S(i));
  },1000);
 }

 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
-template void igl::eytzinger_aabb_sdf<Eigen::Matrix<double, 1, 2, 1, 1, 2>, std::function<double (int)>, Eigen::Matrix<double, -1, 2, 1, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, 1, 2, 1, 1, 2>> const&, std::function<double (int)> 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::Matrix<double, 1, 2, 1, 1, 2>::Scalar&);
+template void igl::eytzinger_aabb_sdf<true, Eigen::Matrix<double, 1, -1, 1, 1, -1>, std::function<double (int)>, Eigen::Matrix<double, -1, 2, 1, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, 1, -1, 1, 1, -1>> const&, std::function<double (int)> 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::Matrix<double, 1, -1, 1, 1, -1>::Scalar&);
 // generated by autoexplicit.sh
-template void igl::eytzinger_aabb_sdf<Eigen::Matrix<double, -1, 3, 1, -1, 3>, std::function<double (Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, int)>, Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>> const&, std::function<double (Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, int)> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>>&);
-template void igl::eytzinger_aabb_sdf<Eigen::Matrix<double, 1, 3, 1, 1, 3>, std::function<double (int)>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3>> const&, std::function<double (int)> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const&, Eigen::Matrix<double, 1, 3, 1, 1, 3>::Scalar&);
-template void igl::eytzinger_aabb_sdf<Eigen::Matrix<double, -1, 3, 0, -1, 3>, std::function<double (Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, int)>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>> const&, std::function<double (Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, int)> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>>&);
+template void igl::eytzinger_aabb_sdf<false, Eigen::Matrix<double, 1, 2, 1, 1, 2>, std::function<double (int)>, Eigen::Matrix<double, -1, 2, 1, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, 1, 2, 1, 1, 2>> const&, std::function<double (int)> 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::Matrix<double, 1, 2, 1, 1, 2>::Scalar&);
+// generated by autoexplicit.sh
+template void igl::eytzinger_aabb_sdf<true, Eigen::Matrix<double, 1, -1, 1, 1, -1>, std::function<double (int)>, Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, 1, -1, 1, 1, -1>> const&, std::function<double (int)> 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::Matrix<double, 1, -1, 1, 1, -1>::Scalar&);
+template void igl::eytzinger_aabb_sdf<false, Eigen::Matrix<double, -1, 3, 1, -1, 3>, std::function<double (Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, int)>, Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>> const&, std::function<double (Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, int)> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>>&);
+template void igl::eytzinger_aabb_sdf<false, Eigen::Matrix<double, 1, 3, 1, 1, 3>, std::function<double (int)>, Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3>> const&, std::function<double (int)> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const&, Eigen::Matrix<double, 1, 3, 1, 1, 3>::Scalar&);
 #endif
@@ -21,7 +21,13 @@ namespace igl
  /// and returns the SDF to that primitive `primitive(i)`
  /// @param[out] f  SDF value at query point x
  ///
+  /// eytzinger_aabb_sdf<false> computes the signed distance function and
+  /// assumes that primitive(i) returns the signed distance to primitive i
+  ///
+  /// eytzinger_aabb_sdf<true> computes the squared (unsigned) distance function
+  /// and assumes that primitive(i) returns the squared distance to primitive i
  template <
+    bool Squared,
    typename Derivedp,
    typename Func,
    typename DerivedB,
@@ -45,6 +51,7 @@ namespace igl
  ///   `primitive(p,i)`
  /// @param[out] S  #P list of SDF values at query points P
  template <
+    bool Squared,
    typename DerivedP,
    typename Func,
    typename DerivedB,
@@ -1,4 +1,26 @@
 #include "eytzinger_aabb_winding_number.h"
+#include "PI.h"
+
+// signed_angle helper. Not sure why igl::signed_angle is doing all that
+// normalization (divisions and square roots) when atan2 doesn't need it.
+namespace 
+{
+  template <typename MatV, typename VecP>
+  inline typename MatV::Scalar signed_angle(
+    const MatV & V,
+    const VecP & p,
+    const int e0,
+    const int e1)
+  {
+    const Eigen::Matrix<typename MatV::Scalar,1,2> v0 = V.row(e0) - p;
+    const Eigen::Matrix<typename MatV::Scalar,1,2> v1 = V.row(e1) - p;
+    const typename MatV::Scalar angle = std::atan2(
+      v0(0)*v1(1) - v0(1)*v1(0),
+      v0(0)*v1(0) + v0(1)*v1(1));
+    return angle;
+  }
+}
+
 template <
  typename Derivedp,
  typename DerivedV,
@@ -18,16 +40,34 @@ IGL_INLINE void igl::eytzinger_aabb_winding_number(
  const Eigen::MatrixBase<DerivedC> & C,
  typename DerivedV::Scalar & wn)
 {
-  const auto signed_angle = [&V,&p]( const int e0, const int e1)->typename DerivedV::Scalar
+  using Scalar = typename DerivedV::Scalar;
+  const auto primitive = [&p,&V,&E](const int e)->typename DerivedV::Scalar
  {
-    const Eigen::Matrix<typename DerivedV::Scalar,1,2> v0 = V.row(e0) - p;
-    const Eigen::Matrix<typename DerivedV::Scalar,1,2> v1 = V.row(e1) - p;
-    const typename DerivedV::Scalar angle = std::atan2(
-      v0(0)*v1(1) - v0(1)*v1(0),
-      v0(0)*v1(0) + v0(1)*v1(1));
-    return angle;
+    return ::signed_angle(V,p,E(e,0),E(e,1)) / Scalar(2.0 * igl::PI);
  };
+  return igl::eytzinger_aabb_winding_number(
+    p,V,primitive,B1,B2,leaf,I,C,wn);
+}

+template <
+  typename Derivedp,
+  typename DerivedV,
+  typename DerivedB,
+  typename Derivedleaf,
+  typename DerivedI,
+  typename DerivedC>
+IGL_INLINE void igl::eytzinger_aabb_winding_number(
+  const Eigen::MatrixBase<Derivedp> & p,
+  const Eigen::MatrixBase<DerivedV> & V,
+  const std::function<typename DerivedV::Scalar(const int)> & primitive,
+  const Eigen::MatrixBase<DerivedB> & B1,
+  const Eigen::MatrixBase<DerivedB> & B2,
+  const Eigen::MatrixBase<Derivedleaf> & leaf,
+  const Eigen::MatrixBase<DerivedI> & I,
+  const Eigen::MatrixBase<DerivedC> & C,
+  typename DerivedV::Scalar & wn)
+{
+  using Scalar = typename DerivedV::Scalar;
  wn = 0;
  if(leaf.size() == 0) { wn = 0; return; }
  // I don't think stack or queue matters
@@ -39,7 +79,8 @@ IGL_INLINE void igl::eytzinger_aabb_winding_number(
    stack.pop_back();
    if(leaf(r) >= 0)
    {
-      wn += signed_angle(E(leaf(r),0),E(leaf(r),1));
+      // Default behavior
+      wn += primitive(leaf(r));
      continue;
    }
    // otherwise is p outside the box B1.row(r), B2.row(r)?
@@ -50,7 +91,7 @@ IGL_INLINE void igl::eytzinger_aabb_winding_number(
      assert(C(r+1) - C(r) % 2 == 0);
      for(int i = C(r); i < C(r+1); i+=2)
      {
-        wn += signed_angle(I(i),I(i+1));
+        wn += ::signed_angle(V,p,I(i),I(i+1)) / Scalar(2.0 * igl::PI);
      }
      continue;
    }
@@ -64,6 +105,7 @@ IGL_INLINE void igl::eytzinger_aabb_winding_number(

 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::eytzinger_aabb_winding_number<Eigen::Matrix<double, 1, -1, 1, 1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 2, 1, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, 1, -1, 1, 1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, std::function<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar (int)> 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::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar&);
 template void igl::eytzinger_aabb_winding_number<Eigen::Matrix<double, 1, 2, 1, 1, 2>, 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<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, 1, 2, 1, 1, 2>> 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::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar&);
 #endif
-
@@ -41,6 +41,23 @@ namespace igl
    const Eigen::MatrixBase<DerivedI> & I,
    const Eigen::MatrixBase<DerivedC> & C,
    typename DerivedV::Scalar & wn);
+  template <
+    typename Derivedp,
+    typename DerivedV,
+    typename DerivedB,
+    typename Derivedleaf,
+    typename DerivedI,
+    typename DerivedC>
+  IGL_INLINE void eytzinger_aabb_winding_number(
+    const Eigen::MatrixBase<Derivedp> & p,
+    const Eigen::MatrixBase<DerivedV> & V,
+    const std::function<typename DerivedV::Scalar(const int)> & primitive,
+    const Eigen::MatrixBase<DerivedB> & B1,
+    const Eigen::MatrixBase<DerivedB> & B2,
+    const Eigen::MatrixBase<Derivedleaf> & leaf,
+    const Eigen::MatrixBase<DerivedI> & I,
+    const Eigen::MatrixBase<DerivedC> & C,
+    typename DerivedV::Scalar & wn);
 }

 #ifndef IGL_STATIC_LIBRARY
@@ -114,5 +114,7 @@ IGL_INLINE void igl::eytzinger_aabb_winding_number_tree(

 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::eytzinger_aabb_winding_number_tree<Eigen::Matrix<int, -1, 2, 1, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<int, -1, 2, 1, -1, 2>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>>&);
 template void igl::eytzinger_aabb_winding_number_tree<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>>&);
 #endif
@@ -0,0 +1,83 @@
+#include "cubic_winding_number.h"
+#include "point_in_convex_hull.h"
+#include "orient2d.h"
+#include "../Orientation.h"
+#include "../cubic_is_flat.h"
+#include "../cubic_split.h"
+#include "../PI.h"
+#include "../matlab_format.h"
+#include <iostream>
+#include <cstdio>
+
+namespace 
+{
+  // helper where we assume that the input is not degenerate 
+  template <typename DerivedC, typename Derivedq>
+  inline typename DerivedC::Scalar cubic_winding_number_helper(
+    const Eigen::MatrixBase<DerivedC>& C,
+    const Eigen::MatrixBase<Derivedq>& q)
+  {
+    using Scalar = typename DerivedC::Scalar;
+
+    if(igl::predicates::point_in_convex_hull(
+        q,
+        C.row(0).eval(),
+        C.row(1).eval(),
+        C.row(2).eval(),
+        C.row(3).eval()) == igl::Orientation::NEGATIVE 
+      )
+    {
+      const auto v0 = (C.row(0) - q).eval();
+      const auto v3 = (C.row(3) - q).eval();
+      const auto w = std::atan2(
+        v0(0)*v3(1) - v0(1)*v3(0),
+        v0(0)*v3(0) + v0(1)*v3(1)) /
+        Scalar(2.0 * igl::PI);
+      return w;
+    }
+    else
+    {
+      Eigen::Matrix<Scalar,4,2,Eigen::RowMajor> C1,C2;
+      igl::cubic_split(C,Scalar(0.5),C1,C2);
+      const auto w1 = cubic_winding_number_helper(C1,q);
+      const auto w2 = cubic_winding_number_helper(C2,q);
+      return w1 + w2;
+    }
+  }
+}
+
+template <typename DerivedC, typename Derivedq>
+IGL_INLINE typename DerivedC::Scalar igl::predicates::cubic_winding_number(
+  const Eigen::MatrixBase<DerivedC>& C,
+  const Eigen::MatrixBase<Derivedq>& q)
+{
+  // Degenerate cases.
+  using Scalar = typename DerivedC::Scalar;
+  // check if C are all the same point
+  if( (C.row(0)-C.row(1)).squaredNorm() == Scalar(0) &&
+      (C.row(0)-C.row(2)).squaredNorm() == Scalar(0) &&
+      (C.row(0)-C.row(3)).squaredNorm() == Scalar(0) )
+  {
+    return Scalar(0);
+  }
+  // Check if C are all COLLINEAR
+  if( 
+      (orient2d(C.row(0),C.row(3),C.row(1)) == Orientation::COLLINEAR) &&
+      (orient2d(C.row(0),C.row(3),C.row(2)) == Orientation::COLLINEAR) )
+  {
+    const auto v0 = (C.row(0) - q).eval();
+    const auto v3 = (C.row(3) - q).eval();
+    const auto w = std::atan2(
+      v0(0)*v3(1) - v0(1)*v3(0),
+      v0(0)*v3(0) + v0(1)*v3(1)) /
+      Scalar(2.0 * igl::PI);
+    return w;
+  }
+  return cubic_winding_number_helper(C,q);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template Eigen::Matrix<double, 4, -1, 1, 4, -1>::Scalar igl::predicates::cubic_winding_number<Eigen::Matrix<double, 4, -1, 1, 4, -1>, Eigen::Matrix<double, 1, -1, 1, 1, -1>>(Eigen::MatrixBase<Eigen::Matrix<double, 4, -1, 1, 4, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, -1, 1, 1, -1>> const&);
+template Eigen::Matrix<double, 4, 2, 0, 4, 2>::Scalar igl::predicates::cubic_winding_number<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::MatrixBase<Eigen::Matrix<double, 1, 2, 1, 1, 2>> const&);
+#endif
@@ -0,0 +1,30 @@
+#ifndef IGL_PREDICATES_CUBIC_WINDING_NUMBER_H
+#define IGL_PREDICATES_CUBIC_WINDING_NUMBER_H
+#include "../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  namespace predicates
+  {
+    /// Computes the (generalized) winding number of a cubic Bézier curve around a query.
+    /// This implementation is similar to "Robust Containment Queries Over
+    /// Collections of Rational Parametric Curves via Generalized Winding
+    /// Numbers" [Spainhour et al. 2024].
+    ///
+    /// @param[in] C  4 by dimensions matrix of control points for a cubic
+    /// Bézier curve
+    /// @param[in] q  1 by dimensions query point
+    template <typename DerivedC, typename Derivedq>
+    IGL_INLINE typename DerivedC::Scalar cubic_winding_number(
+      const Eigen::MatrixBase<DerivedC>& C,
+      const Eigen::MatrixBase<Derivedq>& q);
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "cubic_winding_number.cpp"
+#endif
+
+#endif
+
@@ -46,8 +46,8 @@ IGL_INLINE void igl::predicates::ear_clipping(
    Eigen::Matrix<Scalar,1,2> pb = P.row(b);
    Eigen::Matrix<Scalar,1,2> pi = P.row(i);
    auto r = igl::predicates::orient2d(pa, pi, pb);
-    if(r == igl::predicates::Orientation::NEGATIVE ||
-       r == igl::predicates::Orientation::COLLINEAR) return false;
+    if(r == igl::Orientation::NEGATIVE ||
+       r == igl::Orientation::COLLINEAR) return false;

    // check edge (b, R(b))
    Index k = b;
@@ -57,8 +57,8 @@ IGL_INLINE void igl::predicates::ear_clipping(
    auto r1 = igl::predicates::orient2d(pb, pl, pa);
    auto r2 = igl::predicates::orient2d(pi, pl, pb);
    if (l == a) return true;  // single triangle
-    if (r1 != igl::predicates::Orientation::POSITIVE &&
-        r2 != igl::predicates::Orientation::POSITIVE) {
+    if (r1 != igl::Orientation::POSITIVE &&
+        r2 != igl::Orientation::POSITIVE) {
      return false;
    }
    
@@ -83,8 +83,8 @@ IGL_INLINE void igl::predicates::ear_clipping(
    pl = P.row(l);
    r1 = igl::predicates::orient2d(pb, pk, pl);
    r2 = igl::predicates::orient2d(pi, pa, pk);
-    if (r1 != igl::predicates::Orientation::POSITIVE &&
-        r2 != igl::predicates::Orientation::POSITIVE
+    if (r1 != igl::Orientation::POSITIVE &&
+        r2 != igl::Orientation::POSITIVE
    ){
      return false;
    }
@@ -41,7 +41,7 @@ IGL_INLINE Orientation incircle(
 }

 #ifdef IGL_STATIC_LIBRARY
-#define IGL_INCIRCLE(Vector) template igl::predicates::Orientation igl::predicates::incircle<Vector>(const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&)
+#define IGL_INCIRCLE(Vector) template igl::Orientation igl::predicates::incircle<Vector>(const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&)
 #define IGL_MATRIX(T, R, C) Eigen::Matrix<T, R, C>
 IGL_INCIRCLE(IGL_MATRIX(float, 1, 2));
 IGL_INCIRCLE(IGL_MATRIX(float, 2, 1));
@@ -11,7 +11,7 @@

 #include "../igl_inline.h"
 #include <Eigen/Core>
-#include "Orientation.h"
+#include "../Orientation.h"

 namespace igl {
  namespace predicates {
@@ -44,7 +44,7 @@ IGL_INLINE Orientation insphere(
 }

 #ifdef IGL_STATIC_LIBRARY
-#define IGL_INSPHERE(Vector) template igl::predicates::Orientation igl::predicates::insphere<Vector>(const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&)
+#define IGL_INSPHERE(Vector) template igl::Orientation igl::predicates::insphere<Vector>(const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&)
 #define IGL_MATRIX(T, R, C) Eigen::Matrix<T, R, C>
 IGL_INSPHERE(IGL_MATRIX(float, 1, 3));
 IGL_INSPHERE(IGL_MATRIX(float, 3, 1));
@@ -11,7 +11,7 @@

 #include "../igl_inline.h"
 #include <Eigen/Core>
-#include "Orientation.h"
+#include "../Orientation.h"

 namespace igl {
  namespace predicates {
@@ -16,14 +16,38 @@ namespace predicates {
 using REAL = IGL_PREDICATES_REAL;
 #include "IGL_PREDICATES_ASSERT_SCALAR.h"

-template<typename Vector2D>
+template <
+      typename Derivedpa,
+      typename Derivedpb,
+      typename Derivedpc>
 IGL_INLINE Orientation orient2d(
-    const Eigen::MatrixBase<Vector2D>& pa,
-    const Eigen::MatrixBase<Vector2D>& pb,
-    const Eigen::MatrixBase<Vector2D>& pc)
+    const Eigen::MatrixBase<Derivedpa>& pa,
+    const Eigen::MatrixBase<Derivedpb>& pb,
+    const Eigen::MatrixBase<Derivedpc>& pc)
 {
-  EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Vector2D, 2);
-  IGL_PREDICATES_ASSERT_SCALAR(Vector2D);
+  static_assert(
+      (Derivedpa::RowsAtCompileTime == 2 && Derivedpa::ColsAtCompileTime == 1) ||
+      (Derivedpa::RowsAtCompileTime == 1 && Derivedpa::ColsAtCompileTime == 2) ||
+      (Derivedpa::RowsAtCompileTime == Eigen::Dynamic && Derivedpa::ColsAtCompileTime == 1 ) ||
+      (Derivedpa::RowsAtCompileTime == 1 && Derivedpa::ColsAtCompileTime == Eigen::Dynamic ),
+      "pa must be a 2D point");
+  assert(pa.size() == 2 && "pa must be a 2D point");
+  static_assert(
+      (Derivedpb::RowsAtCompileTime == 2 && Derivedpb::ColsAtCompileTime == 1) ||
+      (Derivedpb::RowsAtCompileTime == 1 && Derivedpb::ColsAtCompileTime == 2) ||
+      (Derivedpb::RowsAtCompileTime == Eigen::Dynamic && Derivedpb::ColsAtCompileTime == 1 ) ||
+      (Derivedpb::RowsAtCompileTime == 1 && Derivedpb::ColsAtCompileTime == Eigen::Dynamic ),
+      "pb must be a 2D point");
+  assert(pb.size() == 2 && "pb must be a 2D point");
+  static_assert(
+      (Derivedpc::RowsAtCompileTime == 2 && Derivedpc::ColsAtCompileTime == 1) ||
+      (Derivedpc::RowsAtCompileTime == 1 && Derivedpc::ColsAtCompileTime == 2) ||
+      (Derivedpc::RowsAtCompileTime == Eigen::Dynamic && Derivedpc::ColsAtCompileTime == 1 ) ||
+      (Derivedpc::RowsAtCompileTime == 1 && Derivedpc::ColsAtCompileTime == Eigen::Dynamic ),
+      "pc must be a 2D point");
+  assert(pc.size() == 2 && "pc must be a 2D point");
+
+

  using Point = Eigen::Matrix<REAL, 2, 1>;
  Point a{pa[0], pa[1]};
@@ -72,14 +96,11 @@ IGL_INLINE void orient2d(
 }

 #ifdef IGL_STATIC_LIBRARY
-#define IGL_ORIENT2D(Vector) template igl::predicates::Orientation igl::predicates::orient2d<Vector>(const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&)
-#define IGL_MATRIX(T, R, C) Eigen::Matrix<T, R, C>
-IGL_ORIENT2D(IGL_MATRIX(float, 1, 2));
-IGL_ORIENT2D(IGL_MATRIX(float, 2, 1));
-#ifndef LIBIGL_PREDICATES_USE_FLOAT
-IGL_ORIENT2D(IGL_MATRIX(double, 1, 2));
-IGL_ORIENT2D(IGL_MATRIX(double, 2, 1));
-#endif
-#undef IGL_MATRIX
-#undef IGL_ORIENT2D
+// Explicit template instantiation
+template igl::Orientation igl::predicates::orient2d<Eigen::Block<Eigen::Matrix<double, 4, -1, 1, 4, -1> const, 1, -1, true>, Eigen::Block<Eigen::Matrix<double, 4, -1, 1, 4, -1> const, 1, -1, true>, Eigen::Block<Eigen::Matrix<double, 4, -1, 1, 4, -1> const, 1, -1, true>>(Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, 4, -1, 1, 4, -1> const, 1, -1, true>> const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, 4, -1, 1, 4, -1> const, 1, -1, true>> const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, 4, -1, 1, 4, -1> const, 1, -1, true>> const&);
+template igl::Orientation igl::predicates::orient2d<Eigen::Block<Eigen::Matrix<double, 4, 2, 0, 4, 2> const, 1, 2, false>, Eigen::Block<Eigen::Matrix<double, 4, 2, 0, 4, 2> const, 1, 2, false>, Eigen::Block<Eigen::Matrix<double, 4, 2, 0, 4, 2> const, 1, 2, false>>(Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, 4, 2, 0, 4, 2> const, 1, 2, false>> const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, 4, 2, 0, 4, 2> const, 1, 2, false>> const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, 4, 2, 0, 4, 2> const, 1, 2, false>> const&);
+template igl::Orientation igl::predicates::orient2d<Eigen::Matrix<double, 1, -1, 1, 1, -1>, Eigen::Matrix<double, 1, -1, 1, 1, -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, 1, -1, 1, 1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, -1, 1, 1, -1>> const&);
+template igl::Orientation igl::predicates::orient2d<Eigen::Matrix<double, 1, 2, 1, 1, 2>, Eigen::Matrix<double, 1, 2, 1, 1, 2>, Eigen::Matrix<double, 1, -1, 1, 1, -1>>(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<double, 1, -1, 1, 1, -1>> const&);
+template igl::Orientation igl::predicates::orient2d<Eigen::Matrix<double, 1, 2, 1, 1, 2>, Eigen::Matrix<double, 1, 2, 1, 1, 2>, Eigen::Matrix<double, 1, 2, 1, 1, 2>>(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<double, 1, 2, 1, 1, 2>> const&);
+template igl::Orientation igl::predicates::orient2d<Eigen::Matrix<double, 2, 1, 0, 2, 1>, Eigen::Matrix<double, 2, 1, 0, 2, 1>, Eigen::Matrix<double, 2, 1, 0, 2, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, 2, 1, 0, 2, 1>> const&, Eigen::MatrixBase<Eigen::Matrix<double, 2, 1, 0, 2, 1>> const&, Eigen::MatrixBase<Eigen::Matrix<double, 2, 1, 0, 2, 1>> const&);
 #endif
@@ -11,7 +11,7 @@

 #include "../igl_inline.h"
 #include <Eigen/Core>
-#include "Orientation.h"
+#include "../Orientation.h"

 namespace igl {
  namespace predicates {
@@ -25,11 +25,14 @@ namespace igl {
    ///          COLLINEAR if they are collinear.
    ///
    /// \fileinfo
-    template<typename Vector2D>
+    template<
+      typename Derivedpa,
+      typename Derivedpb,
+      typename Derivedpc>
    IGL_INLINE Orientation orient2d(
-        const Eigen::MatrixBase<Vector2D>& pa,
-        const Eigen::MatrixBase<Vector2D>& pb,
-        const Eigen::MatrixBase<Vector2D>& pc);
+        const Eigen::MatrixBase<Derivedpa>& pa,
+        const Eigen::MatrixBase<Derivedpb>& pb,
+        const Eigen::MatrixBase<Derivedpc>& pc);
    /// Compute the orientation of the tetrahedron formed by each 4-tuple of
    /// points
    ///
@@ -81,7 +81,7 @@ IGL_INLINE void orient3d(
 // Explicit template instantiation
 template void igl::predicates::orient3d<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<int, -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<double, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>>&);

-#define IGL_ORIENT3D(Vector) template igl::predicates::Orientation igl::predicates::orient3d<Vector>(const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&)
+#define IGL_ORIENT3D(Vector) template igl::Orientation igl::predicates::orient3d<Vector>(const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&)
 #define IGL_MATRIX(T, R, C) Eigen::Matrix<T, R, C>
 IGL_ORIENT3D(IGL_MATRIX(float, 1, 3));
 IGL_ORIENT3D(IGL_MATRIX(float, 3, 1));
@@ -11,7 +11,7 @@

 #include "../igl_inline.h"
 #include <Eigen/Core>
-#include "Orientation.h"
+#include "../Orientation.h"

 namespace igl {
  namespace predicates {
@@ -0,0 +1,204 @@
+#include "point_in_convex_hull.h"
+#include "orient2d.h"
+#include <algorithm>
+#include <cassert>
+
+// Lightly edited LLM code.
+
+namespace 
+{
+  using igl::Orientation;
+template <typename Scalar>
+struct SH
+{
+  static inline int osign(const Orientation o)
+  {
+    using namespace igl::predicates;
+    switch(o)
+    {
+      case Orientation::POSITIVE:  return +1;
+      case Orientation::NEGATIVE:  return -1;
+      case Orientation::COLLINEAR: return  0;
+    }
+    return 0;
+  }
+  
+  static inline bool between_1d(const Scalar x, Scalar lo, Scalar hi)
+  {
+    if(lo > hi) std::swap(lo, hi);
+    return (lo <= x) && (x <= hi);
+  }
+
+  // Assumes q is collinear with segment endpoints (by exact orient2d).
+  static inline bool on_segment_collinear(
+    const Scalar qx, const Scalar qy,
+    const Scalar ax, const Scalar ay,
+    const Scalar bx, const Scalar by)
+  {
+    // Use dominant axis to avoid issues with vertical-ish segments.
+    const Scalar dx = std::abs(bx - ax);
+    const Scalar dy = std::abs(by - ay);
+    if(dx >= dy) return between_1d(qx, ax, bx);
+    else         return between_1d(qy, ay, by);
+  }
+
+  // Closed triangle membership using only cached signs; handles degenerate triples.
+  static inline bool in_triangle_cached(
+    const int o_ab, const int o_bc, const int o_ca,
+    const Scalar qx, const Scalar qy,
+    const Scalar ax, const Scalar ay,
+    const Scalar bx, const Scalar by,
+    const Scalar cx, const Scalar cy)
+  {
+    const bool all_nonneg = (o_ab >= 0 && o_bc >= 0 && o_ca >= 0);
+    const bool all_nonpos = (o_ab <= 0 && o_bc <= 0 && o_ca <= 0);
+    if(!(all_nonneg || all_nonpos))
+      return false;
+
+    // Fully collinear (or coincident) triple: triangle degenerates to a segment/point.
+    if(o_ab == 0 && o_bc == 0 && o_ca == 0)
+    {
+      const Scalar xmin = std::min(ax, std::min(bx, cx));
+      const Scalar xmax = std::max(ax, std::max(bx, cx));
+      const Scalar ymin = std::min(ay, std::min(by, cy));
+      const Scalar ymax = std::max(ay, std::max(by, cy));
+      const Scalar sx = xmax - xmin;
+      const Scalar sy = ymax - ymin;
+
+      if(sx == 0.0 && sy == 0.0)
+        return (qx == ax) && (qy == ay);
+
+      if(sx >= sy) return between_1d(qx, xmin, xmax);
+      else         return between_1d(qy, ymin, ymax);
+    }
+
+    return true;
+  }
+
+};
+}
+
+template <
+  typename Derivedq,
+  typename Deriveda,
+  typename Derivedb,
+  typename Derivedc,
+  typename Derivedd>
+IGL_INLINE igl::Orientation igl::predicates::point_in_convex_hull(
+  const Eigen::MatrixBase<Derivedq> & q,
+  const Eigen::MatrixBase<Deriveda> & a,
+  const Eigen::MatrixBase<Derivedb> & b,
+  const Eigen::MatrixBase<Derivedc> & c,
+  const Eigen::MatrixBase<Derivedd> & d)
+{
+  assert(q.size()==2 && a.size()==2 && b.size()==2 && c.size()==2 && d.size()==2);
+  using Scalar = typename Derivedq::Scalar;
+  using SH = ::SH<Scalar>;
+  using igl::Orientation;
+
+  // For cheap coordinate comparisons (between tests)
+  const Scalar qx = (Scalar)q(0), qy = (Scalar)q(1);
+  const Scalar px[4] = {(Scalar)a(0), (Scalar)b(0), (Scalar)c(0), (Scalar)d(0)};
+  const Scalar py[4] = {(Scalar)a(1), (Scalar)b(1), (Scalar)c(1), (Scalar)d(1)};
+
+  // Cache the 6 unordered edge orientations with q: sign(orient2d(p_i,p_j,q)) for i<j.
+  const int s01 = SH::osign(orient2d(a,b,q));
+  const int s02 = SH::osign(orient2d(a,c,q));
+  const int s03 = SH::osign(orient2d(a,d,q));
+  const int s12 = SH::osign(orient2d(b,c,q));
+  const int s13 = SH::osign(orient2d(b,d,q));
+  const int s23 = SH::osign(orient2d(c,d,q));
+
+
+  const std::function<int(int,int)> Oq = [&](int i, int j)->int
+  {
+    // directed sign(orient2d(p_i,p_j,q)) via antisymmetry
+    if(i == j) return 0;
+    if(i > j) return -Oq(j,i);
+    // i<j
+    if(i==0 && j==1) return s01;
+    if(i==0 && j==2) return s02;
+    if(i==0 && j==3) return s03;
+    if(i==1 && j==2) return s12;
+    if(i==1 && j==3) return s13;
+    /*i==2 && j==3*/ return s23;
+  };
+
+  // Carathéodory: q ∈ hull(4 pts) iff q ∈ some hull(3 pts)
+  auto tri = [&](int ia,int ib,int ic)->bool
+  {
+    const int o_ab = Oq(ia,ib);
+    const int o_bc = Oq(ib,ic);
+    const int o_ca = Oq(ic,ia);
+    return SH::in_triangle_cached(
+      o_ab, o_bc, o_ca,
+      qx, qy,
+      px[ia], py[ia],
+      px[ib], py[ib],
+      px[ic], py[ic]);
+  };
+
+  const bool inside =
+    tri(0,1,2) || tri(0,1,3) || tri(0,2,3) || tri(1,2,3);
+
+  if(!inside)
+    return Orientation::NEGATIVE;
+
+  // --- Boundary classification ---
+  // A segment (i,j) is a supporting segment of conv(P) iff the other two points
+  // lie in the same closed half-plane of the line through (i,j):
+  // sign(orient(p_i,p_j,p_k)) and sign(orient(p_i,p_j,p_l)) are NOT opposite.
+  auto supporting_pair = [&](int i, int j)->bool
+  {
+    int other[2], t=0;
+    for(int k=0;k<4;k++) if(k!=i && k!=j) other[t++] = k;
+
+    // Choose Eigen refs for orient2d calls
+    auto Pi = [&](int idx)->auto const& {
+      switch(idx){
+        case 0: return a;
+        case 1: return b;
+        case 2: return c;
+        default:return d;
+      }
+    };
+
+    const int s1 = SH::osign(orient2d(Pi(i), Pi(j), Pi(other[0])));
+    const int s2 = SH::osign(orient2d(Pi(i), Pi(j), Pi(other[1])));
+
+    // Opposite strict signs => not supporting
+    return !((s1 > 0 && s2 < 0) || (s1 < 0 && s2 > 0));
+  };
+
+  // If hull is fully collinear (degenerate), any in-hull point is "COLLINEAR".
+  const bool nondegenerate =
+    (orient2d(a,b,c) != Orientation::COLLINEAR) ||
+    (orient2d(a,b,d) != Orientation::COLLINEAR) ||
+    (orient2d(a,c,d) != Orientation::COLLINEAR) ||
+    (orient2d(b,c,d) != Orientation::COLLINEAR);
+
+  if(!nondegenerate)
+    return Orientation::COLLINEAR;
+
+  // For non-degenerate hull: q is on boundary iff it lies on some supporting segment.
+  for(int i=0;i<4;i++)
+  for(int j=i+1;j<4;j++)
+  {
+    if(!supporting_pair(i,j)) continue;
+
+    if(Oq(i,j) == 0) // q collinear with (p_i, p_j) by exact orient2d
+    {
+      if(SH::on_segment_collinear(qx,qy, px[i],py[i], px[j],py[j]))
+        return Orientation::COLLINEAR;
+    }
+  }
+
+  return Orientation::POSITIVE; // inside and not on boundary
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template igl::Orientation igl::predicates::point_in_convex_hull<Eigen::Matrix<double, 1, 2, 1, 1, 2>, Eigen::Matrix<double, 1, 2, 1, 1, 2>, Eigen::Matrix<double, 1, 2, 1, 1, 2>, Eigen::Matrix<double, 1, 2, 1, 1, 2>, Eigen::Matrix<double, 1, 2, 1, 1, 2>>(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<double, 1, 2, 1, 1, 2>> const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 2, 1, 1, 2>> const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 2, 1, 1, 2>> const&);
+template igl::Orientation igl::predicates::point_in_convex_hull<Eigen::Matrix<double, 1, -1, 1, 1, -1>, Eigen::Matrix<double, 1, 2, 1, 1, 2>, Eigen::Matrix<double, 1, 2, 1, 1, 2>, Eigen::Matrix<double, 1, 2, 1, 1, 2>, Eigen::Matrix<double, 1, 2, 1, 1, 2>>(Eigen::MatrixBase<Eigen::Matrix<double, 1, -1, 1, 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<double, 1, 2, 1, 1, 2>> const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 2, 1, 1, 2>> const&);
+template igl::Orientation igl::predicates::point_in_convex_hull<Eigen::Matrix<double, 1, -1, 1, 1, -1>, Eigen::Matrix<double, 1, -1, 1, 1, -1>, Eigen::Matrix<double, 1, -1, 1, 1, -1>, Eigen::Matrix<double, 1, -1, 1, 1, -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, 1, -1, 1, 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<double, 1, -1, 1, 1, -1>> const&);
+#endif
@@ -0,0 +1,43 @@
+#ifndef IGL_PREDICATES_POINT_IN_CONVEX_POLYGON_H
+#define IGL_PREDICATES_POINT_IN_CONVEX_POLYGON_H
+
+#include "../igl_inline.h"
+#include <Eigen/Core>
+#include "../Orientation.h"
+
+namespace igl 
+{
+  namespace predicates 
+  {
+    /// Determine if the 2D point q is inside, outside, or on the boundary of the
+    /// convex hull of the 2D points a,b,c,d. The points a,b,c,d may be given in
+    /// any order and the convex hull may be a point, segment, triangle or
+    /// quadrilateral.
+    ///
+    /// @param[in] q  2D query point
+    /// @param[in] a  2D point possibly on the convex hull
+    /// @param[in] b  2D point possibly on the convex hull
+    /// @param[in] c  2D point possibly on the convex hull
+    /// @param[in] d  2D point possibly on the convex hull
+    /// @return Returns +1 if q is inside the convex hull, -1 if outside, and 0
+    /// if on the boundary.
+    template <
+      typename Derivedq,
+      typename Deriveda,
+      typename Derivedb,
+      typename Derivedc,
+      typename Derivedd>
+    IGL_INLINE Orientation point_in_convex_hull(
+      const Eigen::MatrixBase<Derivedq> & q,
+      const Eigen::MatrixBase<Deriveda> & a,
+      const Eigen::MatrixBase<Derivedb> & b,
+      const Eigen::MatrixBase<Derivedc> & c,
+      const Eigen::MatrixBase<Derivedd> & d);
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "point_in_convex_hull.cpp"
+#endif
+
+#endif
@@ -22,8 +22,8 @@ IGL_INLINE bool igl::predicates::point_inside_convex_polygon(
    Eigen::Matrix<Scalar,1,2> a = P.row(i);
    Eigen::Matrix<Scalar,1,2> b = P.row(i_1);
    auto r = igl::predicates::orient2d(a,b,q);
-    if(r == igl::predicates::Orientation::COLLINEAR || 
-       r == igl::predicates::Orientation::NEGATIVE)
+    if(r == igl::Orientation::COLLINEAR || 
+       r == igl::Orientation::NEGATIVE)
      return false;
  }
  return true;
@@ -36,10 +36,10 @@ IGL_INLINE bool igl::predicates::segment_segment_intersect(
  };
  
  // colinear case
-  if((t1 == igl::predicates::Orientation::COLLINEAR && on_segment(a,b,c)) ||
-     (t2 == igl::predicates::Orientation::COLLINEAR && on_segment(c,d,b)) ||
-     (t3 == igl::predicates::Orientation::COLLINEAR && on_segment(a,b,d)) ||
-     (t4 == igl::predicates::Orientation::COLLINEAR && on_segment(c,d,a))) 
+  if((t1 == igl::Orientation::COLLINEAR && on_segment(a,b,c)) ||
+     (t2 == igl::Orientation::COLLINEAR && on_segment(c,d,b)) ||
+     (t3 == igl::Orientation::COLLINEAR && on_segment(a,b,d)) ||
+     (t4 == igl::Orientation::COLLINEAR && on_segment(c,d,a))) 
     return true;
  
  // ordinary case
@@ -0,0 +1,51 @@
+#include "spline_winding_number.h"
+#include "cubic_winding_number.h"
+#include "../parallel_for.h"
+#include "../eytzinger_aabb.h"
+#include "../eytzinger_aabb.h"
+#include "../eytzinger_aabb_winding_number_tree.h"
+#include "../eytzinger_aabb_winding_number.h"
+
+template <
+  typename DerivedP, 
+  typename DerivedC, 
+  typename DerivedB,
+  typename Derivedleaf,
+  typename DerivedQ,
+  typename DerivedW>
+IGL_INLINE void igl::predicates::spline_winding_number(
+  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,
+  const Eigen::MatrixBase<DerivedQ> & Q,
+  Eigen::PlainObjectBase<DerivedW>& W)
+{
+  using Scalar = typename DerivedP::Scalar;
+  Eigen::Matrix<typename DerivedC::Scalar,DerivedC::RowsAtCompileTime,2,Eigen::RowMajor> E(C.rows(),2);
+  E << C.col(0), C.col(3);
+  Eigen::VectorXi tI,tC;
+  eytzinger_aabb_winding_number_tree( E, leaf, tI, tC);
+
+  W.resize(Q.rows(),1);
+  //for(int i = 0;i<Q.rows();i++)
+  igl::parallel_for(Q.rows(),[&](const int i)
+  {
+    const auto qi = Q.row(i).eval();
+    const auto primitive = [&](const int c)
+    {
+      typedef Eigen::Matrix<Scalar,4,DerivedP::ColsAtCompileTime,Eigen::RowMajor> Mat4;
+      // point_spline_squared_distance is caching this slice.
+      const Mat4 Cc = P(C.row(c),Eigen::all);
+      return igl::predicates::cubic_winding_number(Cc, qi);
+    };
+    igl::eytzinger_aabb_winding_number(qi,P,primitive,B1,B2,leaf,tI,tC,W(i));
+  },1000);
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::predicates::spline_winding_number<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<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::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::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>>&);
+#endif
@@ -0,0 +1,50 @@
+#ifndef IGL_PREDICATES_SPLINE_WINDING_NUMBER_H
+#define IGL_PREDICATES_SPLINE_WINDING_NUMBER_H
+#include "../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  namespace predicates
+  {
+    /// Computes the (generalized) winding number of a spline of cubic Bézier
+    /// curves around a set of query points.
+    ///
+    /// @param[in] q  1 by dim query point
+    /// @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[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
+    /// @param[in] Q  #Q by dim matrix of query points
+    /// @param[out] W  #Q by 1 matrix of winding numbers for each query point
+    ///
+    ///
+    /// \see igl::cycodebase::box_cubic, igl::predicates::cubic_winding_number,
+    /// igl::eytzinger_aabb, igl::cycodebase::spline_eytzinger_aabb
+    template <
+      typename DerivedP, 
+      typename DerivedC, 
+      typename DerivedB,
+      typename Derivedleaf,
+      typename DerivedQ,
+      typename DerivedW>
+    IGL_INLINE void spline_winding_number(
+      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,
+      const Eigen::MatrixBase<DerivedQ> & Q,
+      Eigen::PlainObjectBase<DerivedW>& W);
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "spline_winding_number.cpp"
+#endif
+
+#endif
+
+
@@ -21,9 +21,9 @@ IGL_INLINE  bool igl::predicates::triangle_triangle_intersect(
  exactinit();
  using Vector2D = Eigen::Matrix<typename Vector3D::Scalar,2,1>;

-  constexpr Orientation COPLANAR = igl::predicates::Orientation::COPLANAR;
-  constexpr Orientation NEGATIVE = igl::predicates::Orientation::NEGATIVE;
-  constexpr Orientation POSITIVE = igl::predicates::Orientation::POSITIVE;
+  constexpr Orientation COPLANAR = igl::Orientation::COPLANAR;
+  constexpr Orientation NEGATIVE = igl::Orientation::NEGATIVE;
+  constexpr Orientation POSITIVE = igl::Orientation::POSITIVE;

  // Determine for each vertex of one triangle if it's above, below or on the
  // plane of the other triangle.
@@ -6,6 +6,17 @@ namespace igl
 {
  namespace predicates
  {
+    /// Triangle-triangle intersection test using exact predicates.
+    ///
+    /// @param[in] a1 First vertex of triangle A.
+    /// @param[in] a2 Second vertex of triangle A.
+    /// @param[in] a3 Third vertex of triangle A.
+    /// @param[in] b1 First vertex of triangle B.
+    /// @param[in] b2 Second vertex of triangle B.
+    /// @param[in] b3 Third vertex of triangle B.
+    /// @param[out] coplanar True if the triangles are coplanar.
+    /// @return True if the triangles intersect.
+    ///
    template <typename Vector3D>
    IGL_INLINE  bool triangle_triangle_intersect(
      const Vector3D & a1,
@@ -113,12 +113,9 @@ void igl::variable_radius_offset(
  Eigen::PlainObjectBase<DerivedmF> & mF)
 {
  using RowVector3S = Eigen::Matrix<Scalar,1,3>;
-  IGL_TICTOC_LAMBDA;
-  tictoc();
  Eigen::Matrix<Scalar,Eigen::Dynamic,3,Eigen::RowMajor> B1,B2;
  Eigen::VectorXi leaf;
  igl::eytzinger_aabb(PB1,PB2,B1,B2,leaf);
-  printf("%-20s: %g secs\n","eytzinger_aabb",tictoc());

  const std::function<Scalar(const RowVector3S &)>
    sdf = [&](const RowVector3S & p) -> Scalar
@@ -128,7 +125,7 @@ void igl::variable_radius_offset(
      return primitive(p,j);
    };
    Scalar f;
-    igl::eytzinger_aabb_sdf(p,primitive_p,B1,B2,leaf,f);
+    igl::eytzinger_aabb_sdf<false>(p,primitive_p,B1,B2,leaf,f);
    return f;
  };
  const std::function<Scalar(const RowVector3S &)>
@@ -157,18 +154,14 @@ void igl::variable_radius_offset(
  };
  Eigen::Matrix<int,Eigen::Dynamic,3,Eigen::RowMajor> ijk;
  igl::lipschitz_octree( origin,h0,max_depth,udf,ijk);
-  printf("%-20s: %g secs\n","lipschitz_octree",tictoc());

  {
-    tictoc();
    // Gather the corners of those leaf cells
    const Scalar h = h0 / (1 << max_depth);
    Eigen::Matrix<int,Eigen::Dynamic,8,Eigen::RowMajor> J;
    Eigen::Matrix<int,Eigen::Dynamic,3,Eigen::RowMajor> unique_ijk;
    Eigen::Matrix<Scalar,Eigen::Dynamic,3,Eigen::RowMajor> unique_corner_positions;
    igl::unique_sparse_voxel_corners(origin,h0,max_depth,ijk,unique_ijk,J,unique_corner_positions);
-    //printf("unique_sparse_voxel_corners: %0.7f seconds\n",tictoc());
-    printf("%-20s: %g secs\n","unique_sparse_vo...",tictoc());
    /// Evaluate the signed distance function at the corners
    Eigen::VectorXd S(unique_corner_positions.rows());
    //for(int u = 0;u<unique_corner_positions.rows();u++)
@@ -179,12 +172,8 @@ void igl::variable_radius_offset(
        // evaluate the function at the corner
        S(u) = sdf(unique_corner_positions.row(u));
      },1000);
-      //printf("                        sdf: %0.7f seconds\n",tictoc());
-      printf("%-20s: %g secs\n","sdf",tictoc());
    // Run marching cubes on the sparse set of leaf cells
    igl::marching_cubes( S,unique_corner_positions,J, 0, mV,mF);
-    //printf("             marching_cubes: %0.7f seconds\n",tictoc());
-    printf("%-20s: %g secs\n","marching_cubes",tictoc());
  }
 }

@@ -0,0 +1,41 @@
+#include <test_common.h>
+#include <igl/cubic_is_flat.h>
+
+TEST_CASE("cubic_is_flat: simple", "[igl]" )
+{
+  {
+    Eigen::Matrix<double,4,2> C;
+    C << 0,0,
+         1,1,
+         2,-1,
+         3,0;
+    REQUIRE( igl::cubic_is_flat(C,1e-2) == false);
+  }
+  {
+    Eigen::Matrix<double,4,2> C;
+    C << 0,0,
+         1,0.1,
+         2,-0.1,
+         3,0;
+    REQUIRE( igl::cubic_is_flat(C,1e-2) == true);
+  }
+  {
+    Eigen::Matrix<double,4,2> C;
+    C << 0,0,
+         1,0.1,
+         2,-0.1,
+         3,0;
+    REQUIRE( igl::cubic_is_flat(C,1e-4) == false);
+  }
+}
+
+TEST_CASE("cubic_is_flat: degenerate", "[igl]" )
+{
+  Eigen::Matrix<double,4,2> C;
+  C << 0,0,
+    0,0,
+    0,0,
+    0,0;
+  REQUIRE( igl::cubic_is_flat(C,1) == true);
+}
+
@@ -0,0 +1,36 @@
+#include <test_common.h>
+#include <igl/cubic_split.h>
+
+TEST_CASE("cubic_split: simple", "[igl]" )
+{
+  {
+    Eigen::Matrix<double,4,2> C;
+    C << 0,0,
+         1,1,
+         2,-1,
+         3,0;
+    double t = 0.5;
+    Eigen::RowVector2d C01,C012,C0123,C123,C23;
+    igl::cubic_split(C,t,C01,C012,C0123,C123,C23);
+    // C01 = (0.5,0.5)
+    // C012 = (1.0,0.25)
+    // C0123 = (1.5,0.0)
+    // C123 = (2.0,-0.25)
+    // C23 = (2.5,-0.5)
+    REQUIRE( C01.isApprox( Eigen::RowVector2d(0.5,0.5) ) );
+    REQUIRE( C012.isApprox( Eigen::RowVector2d(1.0,0.25) ) );
+    REQUIRE( C0123.isApprox( Eigen::RowVector2d(1.5,0.0) ) );
+    REQUIRE( C123.isApprox( Eigen::RowVector2d(2.0,-0.25) ) );
+    REQUIRE( C23.isApprox( Eigen::RowVector2d(2.5,-0.5) ) );
+    Eigen::Matrix<double,4,2> C1,C2;
+    igl::cubic_split(C,t,C1,C2);
+    REQUIRE( C1.row(0).isApprox( C.row(0) ) );
+    REQUIRE( C1.row(1).isApprox( C01 ) );
+    REQUIRE( C1.row(2).isApprox( C012 ) );
+    REQUIRE( C1.row(3).isApprox( C0123 ) );
+    REQUIRE( C2.row(0).isApprox( C0123 ) );
+    REQUIRE( C2.row(1).isApprox( C123 ) );
+    REQUIRE( C2.row(2).isApprox( C23 ) );
+    REQUIRE( C2.row(3).isApprox( C.row(3) ) );
+  }
+}
@@ -0,0 +1,20 @@
+#include <test_common.h>
+#include <igl/cycodebase/box_cubic.h>
+
+TEST_CASE("box_cubic: simple", "[igl/cycodebase]" )
+{
+  Eigen::MatrixXd C(4,2);
+  C<<0,0,
+    1,2,
+    2,-2,
+    3,0;
+  Eigen::RowVectorXd B1,B2;
+  igl::cycodebase::box_cubic(C,B1,B2);
+  REQUIRE(B1.size() == 2);
+  REQUIRE(B2.size() == 2);
+  REQUIRE(B1[0] == Approx(0.0).margin(1e-12));
+  REQUIRE(B1[1] == Approx(-0.57735026918962584).margin(1e-12));
+  REQUIRE(B2[0] == Approx(3.0).margin(1e-12));
+  REQUIRE(B2[1] == Approx(0.57735026918962584).margin(1e-12));
+}
+
@@ -0,0 +1,39 @@
+#include <test_common.h>
+#include <igl/cycodebase/point_cubic_squared_distance.h>
+
+TEST_CASE("point_cubic_squared_distance: simple2d", "[igl/cycodebase]" )
+{
+  Eigen::MatrixXd C(4,2);
+  C<<0,0,
+    1,2,
+    2,-2,
+    3,0;
+  Eigen::MatrixXd Q(3,2);
+  //[1.5 0;2 0.5;2.5 1]
+  Q<<
+    1.5,0,
+    2,0.5,
+    2.5,1;
+  Eigen::VectorXd sqrD,S;
+  Eigen::MatrixXd K;
+  igl::cycodebase::point_cubic_squared_distance(Q,C,sqrD,S,K);
+  REQUIRE(sqrD.size() == 3);
+  REQUIRE(S.size() == 3);
+  REQUIRE(K.rows() == 3);
+  REQUIRE(K.cols() == 2);
+  // sqrd = [0;0.5;1.25]
+  REQUIRE(sqrD[0] == Approx(0.0).margin(1e-12));
+  REQUIRE(sqrD[1] == Approx(0.5).margin(1e-12));
+  REQUIRE(sqrD[2] == Approx(1.25).margin(1e-12));
+  // S = [0.5;0.5;1.0]
+  REQUIRE(S[0] == Approx(0.5).margin(1e-12));
+  REQUIRE(S[1] == Approx(0.5).margin(1e-12));
+  REQUIRE(S[2] == Approx(1.0).margin(1e-12));
+  // [1.5 0;1.5 0;3 0]
+  REQUIRE(K(0,0) == Approx(1.5).margin(1e-12));
+  REQUIRE(K(0,1) == Approx(0.0).margin(1e-12));
+  REQUIRE(K(1,0) == Approx(1.5).margin(1e-12));
+  REQUIRE(K(1,1) == Approx(0.0).margin(1e-12));
+  REQUIRE(K(2,0) == Approx(3.0).margin(1e-12));
+  REQUIRE(K(2,1) == Approx(0.0).margin(1e-12));
+}
@@ -0,0 +1,43 @@
+#include <test_common.h>
+#include <igl/cycodebase/roots.h>
+
+TEST_CASE("roots: cubic", "[igl/cycodebase]" )
+{
+  // t^3 - 6*t^2 + 11*t - 6
+  {
+    Eigen::Vector<double,4> coef(4);
+    coef << -6,11,-6,1;
+    Eigen::Vector<double,3> R;
+    const int nr = igl::cycodebase::roots(coef, 0.0, 4.0, R);
+    REQUIRE(nr == 3);
+    REQUIRE(R[0] == Approx(1.0).margin(1e-12));
+    REQUIRE(R[1] == Approx(2.0).margin(1e-12));
+    REQUIRE(R[2] == Approx(3.0).margin(1e-12));
+  }
+
+  // With dynamic size
+  {
+    Eigen::VectorXd coef(4);
+    coef << -6,11,-6,1;
+    {
+      Eigen::VectorXd R;
+      const int nr = igl::cycodebase::roots(coef, 0.0, 4.0, R);
+      REQUIRE(nr == 3);
+      REQUIRE(R.size() == 3);
+      REQUIRE(R[0] == Approx(1.0).margin(1e-12));
+      REQUIRE(R[1] == Approx(2.0).margin(1e-12));
+      REQUIRE(R[2] == Approx(3.0).margin(1e-12));
+    }
+    // Only pick the first root
+    {
+      Eigen::VectorXd R;
+      const int nr = igl::cycodebase::roots(coef, 0.0, 1.5, R);
+      REQUIRE(nr == 1);
+      REQUIRE(R.size() == 3);
+      REQUIRE(R[0] == Approx(1.0).margin(1e-12));
+      REQUIRE(std::isnan(R[1]));
+      REQUIRE(std::isnan(R[2]));
+    }
+  }
+}
+
@@ -0,0 +1,77 @@
+#include <test_common.h>
+#include <igl/predicates/cubic_winding_number.h>
+#include <iostream>
+
+TEST_CASE("cubic_winding_number: simple", "[igl/predicates]")
+{
+  Eigen::Matrix<double,4,2> C;
+  C<<0,0,
+    1,1,
+    2,-1,
+    3,0;
+  {
+    Eigen::RowVector2d q(1.1,1.1);
+    double w = igl::predicates::cubic_winding_number(C,q);
+    REQUIRE(0.29147615882815 == Approx(w).epsilon(1e-12));
+  }
+  {
+    Eigen::RowVector2d q(1.45,0.01);
+    double w = igl::predicates::cubic_winding_number(C,q);
+    REQUIRE(-0.502124394734 == Approx(w).epsilon(1e-12));
+  }
+  {
+    Eigen::RowVector2d q(1.45,0);
+    double w = igl::predicates::cubic_winding_number(C,q);
+    REQUIRE(-0.5 == Approx(w).epsilon(1e-12));
+  }
+  // Stress test the numerics
+  {
+    Eigen::RowVector2d q(1.5+1e-15,0);
+    double w = igl::predicates::cubic_winding_number(C,q);
+    REQUIRE(0.5 == Approx(w).epsilon(1e-12));
+  }
+  {
+    Eigen::RowVector2d q(1.5-1e-15,0);
+    double w = igl::predicates::cubic_winding_number(C,q);
+    REQUIRE(-0.5 == Approx(w).epsilon(1e-12));
+  }
+  {
+    Eigen::RowVector2d q(1.5+1e-16,1e-16);
+    double w = igl::predicates::cubic_winding_number(C,q);
+    REQUIRE(0.5 == Approx(w).epsilon(1e-12));
+  }
+  {
+    Eigen::RowVector2d q(1.5-1e-16,-1e-16);
+    double w = igl::predicates::cubic_winding_number(C,q);
+    REQUIRE(-0.5 == Approx(w).epsilon(1e-12));
+  }
+}
+
+TEST_CASE("cubic_winding_number: degenerate", "[igl/predicates]")
+{
+  Eigen::RowVector2d q(1.2,1.0);
+  {
+    Eigen::Matrix<double,4,2> C;
+    C<<0,0,
+      0,0,
+      0,0,
+      0,0;
+    REQUIRE( igl::predicates::cubic_winding_number(C,q) == 0.0);
+  }
+  {
+    Eigen::Matrix<double,4,2> C;
+    C<<0,0,
+      1,0,
+      2,0,
+      3,0;
+    REQUIRE(0.3087217355796 == Approx(igl::predicates::cubic_winding_number(C,q)).epsilon(1e-12));
+  }
+  {
+    Eigen::Matrix<double,4,2> C;
+    C<<0,0,
+      1,1,
+      3,3,
+      9,9;
+    REQUIRE(-0.4835565188700 == Approx(igl::predicates::cubic_winding_number(C,q)).epsilon(1e-12));
+  }
+}
@@ -35,9 +35,9 @@ TEST_CASE("orient3d", "[igl][predicates]")
    Eigen::VectorXi R;
    igl::predicates::orient3d(A,B,C,D,R);
    REQUIRE(R.size() == 3);
-    REQUIRE(R(0) == int(igl::predicates::Orientation::NEGATIVE));
-    REQUIRE(R(1) == int(igl::predicates::Orientation::POSITIVE));
-    REQUIRE(R(2) == int(igl::predicates::Orientation::COPLANAR));
+    REQUIRE(R(0) == int(igl::Orientation::NEGATIVE));
+    REQUIRE(R(1) == int(igl::Orientation::POSITIVE));
+    REQUIRE(R(2) == int(igl::Orientation::COPLANAR));
  }
  {
    Eigen::MatrixXd A(1,3);
@@ -57,9 +57,9 @@ TEST_CASE("orient3d", "[igl][predicates]")
    Eigen::VectorXi R;
    igl::predicates::orient3d(A,B,C,D,R);
    REQUIRE(R.size() == 3);
-    REQUIRE(R(0) == int(igl::predicates::Orientation::NEGATIVE));
-    REQUIRE(R(1) == int(igl::predicates::Orientation::POSITIVE));
-    REQUIRE(R(2) == int(igl::predicates::Orientation::COPLANAR));
+    REQUIRE(R(0) == int(igl::Orientation::NEGATIVE));
+    REQUIRE(R(1) == int(igl::Orientation::POSITIVE));
+    REQUIRE(R(2) == int(igl::Orientation::COPLANAR));
  }
 }

@@ -0,0 +1,83 @@
+#include <test_common.h>
+#include <igl/predicates/point_in_convex_hull.h>
+
+TEST_CASE("point_in_convex_hull: simple", "[igl/predicates]")
+{
+  Eigen::RowVector2d a(0,0);
+  Eigen::RowVector2d b(1,1);
+  Eigen::RowVector2d c(2,-1);
+  Eigen::RowVector2d d(3,0);
+
+  {
+    Eigen::RowVector2d q(1.5,0);
+    auto r = igl::predicates::point_in_convex_hull(q,a,b,c,d);
+    REQUIRE(r == igl::Orientation::POSITIVE);
+  }
+  {
+    Eigen::RowVector2d q(-1,0);
+    auto r = igl::predicates::point_in_convex_hull(q,a,b,c,d);
+    REQUIRE(r == igl::Orientation::NEGATIVE);
+  }
+  {
+    Eigen::RowVector2d q(0,0);
+    auto r = igl::predicates::point_in_convex_hull(q,a,b,c,d);
+    REQUIRE(r == igl::Orientation::COLLINEAR);
+  }
+  {
+    Eigen::RowVector2d q(0.5,0.5);
+    auto r = igl::predicates::point_in_convex_hull(q,a,b,c,d);
+    REQUIRE(r == igl::Orientation::COLLINEAR);
+  }
+  {
+    Eigen::RowVector2d q(1.0,-0.4);
+    auto r = igl::predicates::point_in_convex_hull(q,a,b,c,d);
+    REQUIRE(r == igl::Orientation::POSITIVE);
+  }
+  {
+    Eigen::RowVector2d q(2.0,0.6);
+    auto r = igl::predicates::point_in_convex_hull(q,a,b,c,d);
+    REQUIRE(r == igl::Orientation::NEGATIVE);
+  }
+}
+
+TEST_CASE("point_in_convex_hull: degenerate", "[igl/predicates]")
+{
+  // Triangle
+  {
+    Eigen::RowVector2d a(0,0);
+    Eigen::RowVector2d b(1,0);
+    Eigen::RowVector2d c(0,1);
+    Eigen::RowVector2d d(0.1,0.1);
+    REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(0.1,0.1),a,b,c,d) == igl::Orientation::POSITIVE);
+    REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(0.2,0.2),a,b,c,d) == igl::Orientation::POSITIVE);
+    REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(-0.2,-0.2),a,b,c,d) == igl::Orientation::NEGATIVE);
+    REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(0.5,0.5),a,b,c,d) == igl::Orientation::COLLINEAR);
+  }
+
+  // Segment
+  {
+    Eigen::RowVector2d a(0,0);
+    Eigen::RowVector2d b(1,0);
+    Eigen::RowVector2d c(0.5,0);
+    Eigen::RowVector2d d(0.25,0);
+    REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(0.25,0),a,b,c,d) == igl::Orientation::COLLINEAR);
+    REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(0.5,0),a,b,c,d) == igl::Orientation::COLLINEAR);
+    REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(1.5,0),a,b,c,d) == igl::Orientation::NEGATIVE);
+    REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(-0.5,0),a,b,c,d) == igl::Orientation::NEGATIVE);
+    REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(0.5,1),a,b,c,d) == igl::Orientation::NEGATIVE);
+  }
+
+  // Point
+  {
+    Eigen::RowVector2d a(0,0);
+    Eigen::RowVector2d b(0,0);
+    Eigen::RowVector2d c(0,0);
+    Eigen::RowVector2d d(0,0);
+    REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(0,0),a,b,c,d) == igl::Orientation::COLLINEAR);
+    REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(1,0),a,b,c,d) == igl::Orientation::NEGATIVE);
+    REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(0,1),a,b,c,d) == igl::Orientation::NEGATIVE);
+    REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(-1,0),a,b,c,d) == igl::Orientation::NEGATIVE);
+    REQUIRE(igl::predicates::point_in_convex_hull(Eigen::RowVector2d(0,-1),a,b,c,d) == igl::Orientation::NEGATIVE);
+  }
+
+}
@@ -23,17 +23,17 @@ TEST_CASE("predicates", "[igl][predicates]") {
        e << 0.0, 1.0;
        f << 0.5, 0.5;

-        REQUIRE(orient2d(a, b, c) == Orientation::POSITIVE);
-        REQUIRE(orient2d(a, c, b) == Orientation::NEGATIVE);
-        REQUIRE(orient2d(a, b, b) == Orientation::COLLINEAR);
-        REQUIRE(orient2d(a, a, a) == Orientation::COLLINEAR);
-        REQUIRE(orient2d(a, b, d) == Orientation::COLLINEAR);
-        REQUIRE(orient2d(a, f, c) == Orientation::COLLINEAR);
+        REQUIRE(orient2d(a, b, c) == igl::Orientation::POSITIVE);
+        REQUIRE(orient2d(a, c, b) == igl::Orientation::NEGATIVE);
+        REQUIRE(orient2d(a, b, b) == igl::Orientation::COLLINEAR);
+        REQUIRE(orient2d(a, a, a) == igl::Orientation::COLLINEAR);
+        REQUIRE(orient2d(a, b, d) == igl::Orientation::COLLINEAR);
+        REQUIRE(orient2d(a, f, c) == igl::Orientation::COLLINEAR);

-        REQUIRE(incircle(a,b,c,e) == Orientation::COCIRCULAR);
-        REQUIRE(incircle(a,b,c,a) == Orientation::COCIRCULAR);
-        REQUIRE(incircle(a,b,c,d) == Orientation::OUTSIDE);
-        REQUIRE(incircle(a,b,c,f) == Orientation::INSIDE);
+        REQUIRE(incircle(a,b,c,e) == igl::Orientation::COCIRCULAR);
+        REQUIRE(incircle(a,b,c,a) == igl::Orientation::COCIRCULAR);
+        REQUIRE(incircle(a,b,c,d) == igl::Orientation::OUTSIDE);
+        REQUIRE(incircle(a,b,c,f) == igl::Orientation::INSIDE);
    }

    SECTION("3D") {
@@ -46,16 +46,16 @@ TEST_CASE("predicates", "[igl][predicates]") {
        e << 1.0, 1.0, 1.0;
        f << std::numeric_limits<Scalar>::epsilon(), 0.0, 0.0;

-        REQUIRE(orient3d(a, b, c, d) == Orientation::NEGATIVE);
-        REQUIRE(orient3d(a, b, d, c) == Orientation::POSITIVE);
-        REQUIRE(orient3d(a, b, d, d) == Orientation::COPLANAR);
-        REQUIRE(orient3d(a, a, a, a) == Orientation::COPLANAR);
-        REQUIRE(orient3d(a, b, f, c) == Orientation::COPLANAR);
+        REQUIRE(orient3d(a, b, c, d) == igl::Orientation::NEGATIVE);
+        REQUIRE(orient3d(a, b, d, c) == igl::Orientation::POSITIVE);
+        REQUIRE(orient3d(a, b, d, d) == igl::Orientation::COPLANAR);
+        REQUIRE(orient3d(a, a, a, a) == igl::Orientation::COPLANAR);
+        REQUIRE(orient3d(a, b, f, c) == igl::Orientation::COPLANAR);

-        REQUIRE(insphere(a, b, c, d, e) == Orientation::COSPHERICAL);
-        REQUIRE(insphere(a, b, d, e, c) == Orientation::COSPHERICAL);
-        REQUIRE(insphere(b, c, e, d, ((a+b)*0.5).eval()) == Orientation::INSIDE);
-        REQUIRE(insphere(b, c, e, d, (-f).eval()) == Orientation::OUTSIDE);
-        REQUIRE(insphere(f, b, d, c, e) == Orientation::INSIDE);
+        REQUIRE(insphere(a, b, c, d, e)                  == igl::Orientation::COSPHERICAL);
+        REQUIRE(insphere(a, b, d, e, c)                  == igl::Orientation::COSPHERICAL);
+        REQUIRE(insphere(b, c, e, d, ((a+b)*0.5).eval()) == igl::Orientation::INSIDE);
+        REQUIRE(insphere(b, c, e, d, (-f).eval())        == igl::Orientation::OUTSIDE);
+        REQUIRE(insphere(f, b, d, c, e)                  == igl::Orientation::INSIDE);
    }
 }
@@ -178,7 +178,7 @@ int main(int argc, char * argv[])

  if(offset == OffsetType::VARIABLE_RADIUS_OFFSET)
  {
-    viewer.data().set_colors(R);
+    viewer.data().set_data(R);
  }
  viewer.append_mesh();
  const int m_index = viewer.selected_data_index;
@@ -207,7 +207,7 @@ int main(int argc, char * argv[])
      GV *= 2;
      GV.array() -= 1;
      Eigen::VectorXd S;
-      igl::eytzinger_aabb_sdf(GV,primitive,B1,B2,leaf,S);
+      igl::eytzinger_aabb_sdf<false>(GV,primitive,B1,B2,leaf,S);
      printf("%-20s: %g secs\n","eytzinger_aabb_sdf",tictoc());

      tictoc();
@@ -231,7 +231,7 @@ int main(int argc, char * argv[])
        return primitive(p,j);
      };
      double f;
-      igl::eytzinger_aabb_sdf(p,primitive_p,B1,B2,leaf,f);
+      igl::eytzinger_aabb_sdf<false>(p,primitive_p,B1,B2,leaf,f);
      return f - isovalue;
    };
    const std::function<double(const Eigen::RowVector3d &)>
@@ -0,0 +1,287 @@
+#include <igl/read_triangle_mesh.h>
+#include <igl/box_simplices.h>
+#include <igl/eytzinger_aabb.h>
+#include <igl/lipschitz_octree.h>
+#include <igl/unique_sparse_voxel_corners.h>
+#include <igl/eytzinger_aabb_sdf.h>
+#include <igl/unique.h>
+#include <igl/readDMAT.h>
+#include <igl/grid.h>
+#include <igl/edges.h>
+#include <igl/marching_cubes.h>
+#include <igl/parallel_for.h>
+#include <igl/opengl/glfw/Viewer.h>
+#include <igl/matlab_format.h>
+#include <igl/variable_radius_offset.h>
+#include <igl/get_seconds.h>
+#include <igl/SphereMeshWedge.h>
+
+
+#include <limits>
+
+#include <Eigen/Geometry>
+#define sign(x) ((x>0) - (x<0))
+
+
+inline double dot2(const Eigen::RowVector3d & v)
+{
+  return v.dot(v);
+}
+
+inline double clamp(const double & x, const double & a, const double & b)
+{
+  return std::max(a, std::min(b, x));
+}
+
+// https://iquilezles.org/articles/triangledistance/
+template <typename T1, typename T2, typename T3, typename Tp>
+typename Tp::Scalar udTriangle(
+    const T1 & v1, const T2 & v2, const T3 & v3, const Tp & p)
+{
+  using vec3 = Eigen::Matrix<typename Tp::Scalar, 3, 1>;
+  // prepare data    
+  vec3 v21 = v2 - v1; vec3 p1 = p - v1;
+  vec3 v32 = v3 - v2; vec3 p2 = p - v2;
+  vec3 v13 = v1 - v3; vec3 p3 = p - v3;
+  vec3 nor =  v21.cross( v13 );
+
+    return sqrt( // inside/outside test    
+                 (sign(v21.cross(nor).dot(p1)) + 
+                  sign(v32.cross(nor).dot(p2)) + 
+                  sign(v13.cross(nor).dot(p3))<2.0) 
+                  ?
+                  // 3 edges    
+                  std::min( std::min( 
+                  dot2(v21*clamp(v21.dot(p1)/dot2(v21),0.0,1.0)-p1), 
+                  dot2(v32*clamp(v32.dot(p2)/dot2(v32),0.0,1.0)-p2) ), 
+                  dot2(v13*clamp(v13.dot(p3)/dot2(v13),0.0,1.0)-p3) )
+                  :
+                  // 1 face    
+                  nor.dot(p1)*nor.dot(p1)/dot2(nor) );
+}
+    
+
+
+
+template <typename Tp, typename Ta, typename Tb>
+typename Tp::Scalar udLineSegment(const Tp & p, const Ta & a, const Tb & b)
+{
+  using vec3 = Eigen::Matrix<typename Tp::Scalar, 3, 1>;
+  vec3 pa = p - a;
+  vec3 ba = b - a;
+  const auto h = clamp( pa.dot(ba)/ba.dot(ba), 0.0, 1.0 );
+  return ( pa - ba*h ).norm();
+}
+
+
+
+int main(int argc, char * argv[])
+{
+  const int ns = 256+1;
+  Eigen::RowVector3i res(ns,ns,ns);
+
+  IGL_TICTOC_LAMBDA;
+  // Read mesh
+  Eigen::Matrix<double,Eigen::Dynamic,3,Eigen::RowMajor> V;
+  Eigen::Matrix<int,Eigen::Dynamic,3,Eigen::RowMajor> F;
+  if(!igl::read_triangle_mesh
+     (argc>1?argv[1]: TUTORIAL_SHARED_PATH "/decimated-knight.off",V,F)) {
+    std::cout << "Failed to load mesh." << std::endl;
+  }
+  // Normalize
+  V.rowwise() -= 0.5* (V.colwise().maxCoeff() + V.colwise().minCoeff());
+  const double scale = 0.8/V.array().abs().maxCoeff();
+  V *= scale;
+
+  // if second arg exists read it to dmat
+  Eigen::VectorXd R;
+  if(argc>2)
+  {
+    igl::readDMAT( argv[2],R);
+    R *= scale;
+  }else
+  {
+    //Eigen::VectorXd R = (0.1+0.9*(0.5+(Eigen::VectorXd::Random(V.rows())*0.5).array()))*((2.0/(res(0)-1))*2);
+    R = V.col(1);
+    R.array() -= R.minCoeff();
+    R /= R.maxCoeff();
+    R.array() += 1.0;
+    R = (R.array().pow(4)).eval();
+    R *= 0.01;
+  }
+
+
+  tictoc();
+  Eigen::Matrix<double,Eigen::Dynamic,3,Eigen::RowMajor> PB1,PB2;
+  Eigen::Matrix<int,Eigen::Dynamic,2,Eigen::RowMajor> E;
+  igl::edges(F,E);
+
+  enum class OffsetType
+  {
+    EDGE_OFFSET = 0,
+    TRIANGLE_OFFSET = 1,
+    VARIABLE_RADIUS_OFFSET = 2
+  } offset = OffsetType::VARIABLE_RADIUS_OFFSET;
+
+  double isovalue = 0.16;(2.0/(res(0)-1))*2;
+  std::function<double(const Eigen::RowVector3d &,const int i)> primitive;
+  std::vector<igl::SphereMeshWedge<double>> data;
+  switch(offset)
+  {
+    case OffsetType::EDGE_OFFSET:
+    {
+      igl::box_simplices(V,E,PB1,PB2);
+      primitive = [&](const Eigen::RowVector3d & p, const int i)
+      {
+        return udLineSegment(p, V.row(E(i,0)), V.row(E(i,1)));
+      };
+      break;
+    }
+    case OffsetType::TRIANGLE_OFFSET:
+    {
+      igl::box_simplices(V,F,PB1,PB2);
+      primitive = [&](const Eigen::RowVector3d & p, const int i)
+      {
+        return udTriangle( V.row(F(i,0)), V.row(F(i,1)), V.row(F(i,2)), p);
+      };
+      break;
+    }
+    case OffsetType::VARIABLE_RADIUS_OFFSET:
+    {
+      isovalue = 0;
+      igl::variable_radius_offset(V,F,R,PB1,PB2,data,primitive);
+      break;
+    }
+  }
+  printf("%-20s: %g secs\n","PB",tictoc());
+
+  tictoc();
+  Eigen::Matrix<double,Eigen::Dynamic,3,Eigen::RowMajor> B1,B2;
+  Eigen::VectorXi leaf;
+  igl::eytzinger_aabb(PB1,PB2,B1,B2,leaf);
+  printf("%-20s: %g secs\n","eytzinger_aabb",tictoc());
+  //{
+  //  Eigen::VectorXi U;
+  //  igl::unique(leaf,U);
+  //  printf("%d → %d\n",PB1.rows(),U.size()-2);
+  //}
+
+  tictoc();
+  Eigen::Matrix<double,Eigen::Dynamic,3> mV;
+  Eigen::Matrix<int,Eigen::Dynamic,3> mF;
+  if(ns <= 64+1)
+  {
+    Eigen::Matrix<double,Eigen::Dynamic,3,Eigen::RowMajor> GV;
+    igl::grid(res,GV);
+    GV *= 2;
+    GV.array() -= 1;
+    Eigen::VectorXd S;
+    igl::eytzinger_aabb_sdf(GV,primitive,B1,B2,leaf,S);
+    printf("%-20s: %g secs\n","eytzinger_aabb_sdf",tictoc());
+
+    tictoc();
+    igl::marching_cubes(S,GV,res(0),res(1),res(2),isovalue,mV,mF);
+    printf("%-20s: %g secs\n","marching_cubes",tictoc());
+  }
+
+  tictoc();
+  const std::function<double(const Eigen::RowVector3d &)>
+    sdf = [&](const Eigen::RowVector3d & p) -> double
+  {
+    const std::function<double(const int)> primitive_p = [&](const int j)
+    {
+      return primitive(p,j);
+    };
+    double f;
+    igl::eytzinger_aabb_sdf(p,primitive_p,B1,B2,leaf,f);
+    return f;
+  };
+  const std::function<double(const Eigen::RowVector3d &)>
+    udf = [&](const Eigen::RowVector3d & p) -> double
+  {
+    return std::abs(sdf(p));
+    //// This made performance worse.
+    //const std::function<double(const int)> primitive_p = [&](const int j)
+    //{
+    //  const double d = udTriangle( V.row(F(j,0)), V.row(F(j,1)), V.row(F(j,2)), p);
+    //  const Eigen::RowVector3d r(R(F(j,0)),R(F(j,1)),R(F(j,2)));
+    //  const double min_r = r.minCoeff();
+    //  const double max_r = r.maxCoeff();
+    //  if(d > max_r)
+    //  {
+    //    return d - max_r;
+    //  }else if(d < min_r)
+    //  {
+    //    return d - min_r;
+    //  }
+    //  return 0.0;
+    //};
+    //double f;
+    //igl::eytzinger_aabb_sdf(p,primitive_p,B1,B2,leaf,f);
+    //return f;
+  };
+  Eigen::RowVector3d origin(-1,-1,-1);
+  const double h0 = 2;
+  const int max_depth = floor(log2(res(0)));
+  Eigen::Matrix<int,Eigen::Dynamic,3,Eigen::RowMajor> ijk;
+  igl::lipschitz_octree( origin,h0,max_depth,udf,ijk);
+  printf("%-20s: %g secs\n","lipschitz_octree",tictoc());
+
+  Eigen::Matrix<double,Eigen::Dynamic,3> oV;
+  Eigen::Matrix<int,Eigen::Dynamic,3> oF;
+  {
+    tictoc();
+    // Gather the corners of those leaf cells
+    const double h = h0 / (1 << max_depth);
+    Eigen::Matrix<int,Eigen::Dynamic,8,Eigen::RowMajor> J;
+    Eigen::Matrix<int,Eigen::Dynamic,3,Eigen::RowMajor> unique_ijk;
+    Eigen::Matrix<double,Eigen::Dynamic,3,Eigen::RowMajor> unique_corner_positions;
+    igl::unique_sparse_voxel_corners(origin,h0,max_depth,ijk,unique_ijk,J,unique_corner_positions);
+    //printf("unique_sparse_voxel_corners: %0.7f seconds\n",tictoc());
+    printf("%-20s: %g secs\n","unique_sparse_vo...",tictoc());
+    /// Evaluate the signed distance function at the corners
+    Eigen::VectorXd S(unique_corner_positions.rows());
+    //for(int u = 0;u<unique_corner_positions.rows();u++)
+    igl::parallel_for(
+      unique_corner_positions.rows(),
+      [&](const int u)
+      {
+        // evaluate the function at the corner
+        S(u) = sdf(unique_corner_positions.row(u));
+      },1000);
+      //printf("                        sdf: %0.7f seconds\n",tictoc());
+      printf("%-20s: %g secs\n","sdf",tictoc());
+    // Run marching cubes on the sparse set of leaf cells
+    igl::marching_cubes( S,unique_corner_positions,J, 0, oV,oF);
+    //printf("             marching_cubes: %0.7f seconds\n",tictoc());
+    printf("%-20s: %g secs\n","marching_cubes",tictoc());
+  }
+
+  //tictoc();
+  //printf("------------------\n");
+  //igl::variable_radius_offset(V,F,R,origin,h0,max_depth,oV,oF);
+  //printf("%-20s: %g secs\n","variable_radius_offset",tictoc());
+
+
+  igl::opengl::glfw::Viewer viewer;
+  viewer.data().set_mesh(V,F);
+  viewer.data().set_face_based(true);
+  viewer.data().show_lines = false;
+
+  viewer.append_mesh();
+  viewer.data().set_mesh(mV,mF);
+  viewer.data().set_face_based(true);
+  viewer.data().show_lines = true;
+  viewer.data().show_faces = true;
+  viewer.data().set_colors(Eigen::RowVector3d(0.8,0.2,0.2));
+
+  viewer.append_mesh();
+  viewer.data().set_mesh(oV,oF);
+  viewer.data().set_face_based(true);
+  viewer.data().show_lines = true;
+  viewer.data().show_faces = true;
+  viewer.data().set_colors(Eigen::RowVector3d(0.2,0.8,0.2));
+
+  viewer.launch();
+  return 0;
+}
@@ -118,7 +118,7 @@ int main(int argc, char * argv[])
        };
        if(quantity == UNSIGNED_DISTANCE || quantity == SIGNED_DISTANCE)
        {
-          igl::eytzinger_aabb_sdf(p,primitive_p,B1,B2,leaf,di);
+          igl::eytzinger_aabb_sdf<false>(p,primitive_p,B1,B2,leaf,di);
        }
        if(quantity == WINDING_NUMBER || quantity == SIGNED_DISTANCE)
        {
@@ -0,0 +1,186 @@
+#include <igl/opengl/glfw/Viewer.h>
+#include <igl/get_seconds.h>
+#include <igl/find.h>
+#include <igl/readDMAT.h>
+#include <igl/writeDMAT.h>
+#include <igl/triangulated_grid.h>
+#include <igl/cubic.h>
+#include <igl/cycodebase/box_cubic.h>
+#include <igl/cycodebase/point_spline_squared_distance.h>
+#include <igl/cycodebase/spline_eytzinger_aabb.h>
+#include <igl/predicates/spline_winding_number.h>
+
+int main(int argc, char * argv[])
+{
+  IGL_TICTOC_LAMBDA;
+  Eigen::MatrixXd P;
+  Eigen::MatrixXi C;
+  igl::readDMAT(argc>1?argv[1]: TUTORIAL_SHARED_PATH "/libigl-cursive-P.dmat",P);
+  igl::readDMAT(argc>1?argv[2]: TUTORIAL_SHARED_PATH "/libigl-cursive-C.dmat",C);
+
+  // Cheezy conversion of (P,C) to polyline using 50 samples per curve
+  const int ns = 50;
+  Eigen::MatrixXd V(ns*C.rows(),P.cols());
+  Eigen::MatrixXi E((ns-1)*C.rows(),2);
+  for(int c = 0;c<C.rows();c++)
+  {
+    for(int i = 0;i<ns;i++)
+    {
+      const double t = double(i)/(ns-1);
+      // Evaluate cubic Bezier at parameter t
+      Eigen::RowVectorXd Vct;
+      igl::cubic(P(C.row(c),Eigen::all),t,Vct);
+      V.row(c*ns + i) = Vct;
+      if(i>0)
+      {
+        E.row(c*(ns-1) + (i-1))<< c*ns + (i-1), c*ns + i;
+      }
+    }
+  }
+
+
+  Eigen::MatrixXd GV;
+  Eigen::MatrixXi GF;
+  {
+    Eigen::RowVector2d min_corner = P.colwise().minCoeff();
+    Eigen::RowVector2d max_corner = P.colwise().maxCoeff();
+    // samples on x-axis
+    const int sx = 512;
+    const int sy = sx * (max_corner(1)-min_corner(1)) / (max_corner(0)-min_corner(0));
+    Eigen::RowVector2d diag = (max_corner - min_corner).eval();
+    const double bbd = diag.norm();
+    diag /= diag.norm();
+    min_corner -= 0.1*bbd*diag;
+    max_corner += 0.1*bbd*diag;
+
+    igl::triangulated_grid(sx,sy,GV,GF);
+    // Scale and translate grid to fit bounding box
+    GV.col(0) = (GV.col(0).array() * (max_corner(0)-min_corner(0))) + min_corner(0);
+    GV.col(1) = (GV.col(1).array() * (max_corner(1)-min_corner(1))) + min_corner(1);
+  }
+
+
+  Eigen::VectorXd sqrD,S;
+  Eigen::VectorXi I;
+  Eigen::MatrixXd K;
+  Eigen::VectorXd W;
+
+  Eigen::Matrix<double,Eigen::Dynamic,2,Eigen::RowMajor> B1,B2;
+  Eigen::VectorXi leaf;
+  tictoc();
+  igl::cycodebase::spline_eytzinger_aabb(P, C, B1, B2,leaf);
+  printf("AABB build time: %g secs\n",tictoc());
+  igl::writeDMAT("P.dmat", Eigen::MatrixXd(P));
+  igl::writeDMAT("C.dmat", Eigen::MatrixXi(C));
+  igl::writeDMAT("B1.dmat",Eigen::MatrixXd(B1));
+  igl::writeDMAT("B2.dmat",Eigen::MatrixXd(B2));
+
+  {
+    tictoc();
+    igl::cycodebase::point_spline_squared_distance(
+        GV, P, C, B1,B2,leaf, sqrD, I, S, K);
+    printf("sqrD: %d points in %g secs\n",GV.rows(),tictoc());
+    tictoc();
+    igl::predicates::spline_winding_number(P,C,B1,B2,leaf,GV,W);
+    printf("Wind: %d points in %g secs\n",GV.rows(),tictoc());
+  }
+  Eigen::VectorXd unsigned_D = sqrD.array().sqrt();
+  Eigen::VectorXd signed_D = unsigned_D.array() * (0.5 - W.array().abs()) * 2.0;
+
+  // Just leaves
+  //Eigen::MatrixXd B1,B2;
+  //igl::cycodebase::box_cubic(P,C,B1,B2);
+
+  const auto nonempty = igl::find((leaf.array() != -2).eval());
+  const int nb = nonempty.size();
+  Eigen::MatrixXd BV(nb*4,2);
+  BV<<
+    B1(nonempty,0), B1(nonempty,1),
+    B1(nonempty,0), B2(nonempty,1),
+    B2(nonempty,0), B2(nonempty,1),
+    B2(nonempty,0), B1(nonempty,1);
+  Eigen::MatrixXi BE(nb*4,2);
+  for(int c = 0;c<nb;c++)
+  {
+    BE.row(c*4 + 0)<< 0*nb + c, 1*nb + c;
+    BE.row(c*4 + 1)<< 1*nb + c, 2*nb + c;
+    BE.row(c*4 + 2)<< 2*nb + c, 3*nb + c;
+    BE.row(c*4 + 3)<< 3*nb + c, 0*nb + c;
+  }
+
+  igl::opengl::glfw::Viewer vr;
+  vr.data().set_mesh(GV,GF);
+  vr.data().show_lines = false;
+  const int g_index = vr.selected_data_index;
+  vr.append_mesh();
+  vr.data().set_mesh(V,E(Eigen::all,{0,1,1}).eval());
+  vr.data().set_points(P,Eigen::RowVector3d(1,0.7,0.2));
+  Eigen::MatrixXi CE(C.rows()*2,2);
+  CE<< C(Eigen::all,{0,1}),
+    C(Eigen::all,{2,3});
+  vr.data().set_edges(P,CE,Eigen::RowVector3d(1,0.7,0.2));
+  vr.data().point_size = 10;
+  vr.append_mesh();
+  vr.data().set_edges(BV,BE,Eigen::RowVector3d(0,0,0));
+
+  enum Quantity
+  {
+    SIGNED_DISTANCE,
+    UNSIGNED_DISTANCE,
+    WINDING_NUMBER
+  } quantity = SIGNED_DISTANCE;
+
+  const auto update = [&]()
+  {
+    vr.core().lighting_factor = 0;
+    const double dmax = unsigned_D.maxCoeff();
+    const double dmin = -dmax;
+    vr.data_list[g_index].show_texture = false;
+    switch(quantity)
+    {
+      case SIGNED_DISTANCE:
+        vr.data_list[g_index].set_data(signed_D,dmin,dmax,igl::COLOR_MAP_TYPE_ZOE,32);
+        break;
+      case UNSIGNED_DISTANCE:
+        vr.data_list[g_index].set_data(unsigned_D,dmin,dmax,igl::COLOR_MAP_TYPE_ZOE,32);
+        break;
+      case WINDING_NUMBER:
+      {
+        vr.data_list[g_index].set_data(W,W.minCoeff(),W.maxCoeff(),igl::COLOR_MAP_TYPE_TURBO,std::round(W.maxCoeff()-W.minCoeff())+1);
+        break;
+      }
+    }
+  };
+  update();
+
+  vr.callback_key_pressed = [&](igl::opengl::glfw::Viewer &, unsigned int key, int mod)->bool
+  {
+    switch(key)
+    {
+      default:
+        return false;
+      case 'W':
+      case 'w':
+        quantity = WINDING_NUMBER;
+        break;
+      case 'U':
+      case 'u':
+        quantity = UNSIGNED_DISTANCE;
+        break;
+      case 'S':
+      case 's':
+        quantity = SIGNED_DISTANCE;
+        break;
+    }
+    update();
+    return true;
+  };
+  std::cout<<R"(
+  S,s     Signed distance
+  W,w     Winding number
+  U,u     Unsigned distance
+  N,n     Toggle naive / accelerated
+  )"<<std::endl;
+  vr.launch();
+  return 0;
+}
@@ -137,4 +137,5 @@ if(LIBIGL_TUTORIALS_CHAPTER10)
  igl_add_tutorial(1001_LipschitzOctree igl::glfw)
  igl_add_tutorial(1002_EytzingerAABB igl::glfw)
  igl_add_tutorial(1003_WindingNumber2D igl::glfw)
+  igl_add_tutorial(1004_SplineDistance igl::glfw igl::cycodebase igl::predicates)
 endif()
Author	SHA1	Message	Date
Alec Jacobson	e0f3fe6688	better cm	2026-01-26 19:43:58 -05:00
Alec Jacobson	0e8f09a4f8	doc	2026-01-26 18:21:23 -05:00
Alec Jacobson	a439f89d22	working spline winding number and demo	2026-01-26 17:11:02 -05:00
Alec Jacobson	2ecacc3c76	Orientation to igl::, point_in_convex_hull, eyt_winding_number with func handle, fix bug in eyt_sdf, tutorial, tests	2026-01-26 12:13:43 -05:00
Alec Jacobson	0496a74ed7	roots, cubics, tests and tutorial	2026-01-24 23:27:28 -05:00