cycodebase: reuse AABB / precomputed bases for fast repeated distance queries

src/cycodebase/point_cubic_squared_distance.cpp

@@ -19,6 +19,20 @@ namespace pyigl
     igl::cycodebase::point_cubic_squared_distance(Q, C, sqrD, S, K);
     return std::make_tuple(sqrD, S, K);
   }
+
+  // Precomputed overload: reuse the monomial bases (D, B) from
+  // igl.cubic_monomial_bases so repeated queries don't recompute them.
+  auto point_cubic_squared_distance_bases(
+    const nb::DRef<const Eigen::MatrixXN> &Q,
+    const nb::DRef<const Eigen::MatrixXN> &C,
+    const nb::DRef<const Eigen::MatrixXN> &D,
+    const nb::DRef<const Eigen::VectorXN> &B)
+  {
+    Eigen::VectorXN sqrD, S;
+    Eigen::MatrixXN K;
+    igl::cycodebase::point_cubic_squared_distance(Q, C, D, B, sqrD, S, K);
+    return std::make_tuple(sqrD, S, K);
+  }
 }
 
 void bind_point_cubic_squared_distance(nb::module_ &m)
@@ -29,6 +43,23 @@ void bind_point_cubic_squared_distance(nb::module_ &m)
 
 @param[in] Q  #Q by dim matrix of query points
 @param[in] C  4 by dim matrix of control points for the cubic Bézier curve
+@return Tuple (sqrD, S, K) where
+  sqrD  #Q vector of smallest squared distances
+  S     #Q vector of parameters of the closest points on the curve
+  K     #Q by dim matrix of closest points on the curve)");
+
+  m.def("point_cubic_squared_distance", &pyigl::point_cubic_squared_distance_bases,
+    "Q"_a, "C"_a, "D"_a, "B"_a,
+    R"(Squared distance to a cubic Bézier curve, reusing precomputed monomial bases.
+
+Compute the bases once with igl.cubic_monomial_bases(C) -> (M, D, B) and pass
+(D, B) here so repeated queries against the same curve skip recomputing them.
+The result matches the two-argument overload.
+
+@param[in] Q  #Q by dim matrix of query points
+@param[in] C  4 by dim matrix of control points for the cubic Bézier curve
+@param[in] D  3 by dim matrix of monomial coefficients for dC/dt (from cubic_monomial_bases)
+@param[in] B  6-vector of inner products of the monomial bases (from cubic_monomial_bases)
 @return Tuple (sqrD, S, K) where
   sqrD  #Q vector of smallest squared distances
   S     #Q vector of parameters of the closest points on the curve

src/cycodebase/point_spline_squared_distance.cpp

@@ -21,6 +21,24 @@ namespace pyigl
     igl::cycodebase::point_spline_squared_distance(Q, P, C, sqrD, I, S, K);
     return std::make_tuple(sqrD, I, S, K);
   }
+
+  // Accelerated overload: reuse a prebuilt Eytzinger AABB (B1, B2, leaf) from
+  // igl.cycodebase.spline_eytzinger_aabb so repeated queries don't rebuild it.
+  auto point_spline_squared_distance_aabb(
+    const nb::DRef<const Eigen::MatrixXN> &Q,
+    const nb::DRef<const Eigen::MatrixXN> &P,
+    const nb::DRef<const Eigen::MatrixXI> &C,
+    const nb::DRef<const Eigen::MatrixXN> &B1,
+    const nb::DRef<const Eigen::MatrixXN> &B2,
+    const nb::DRef<const Eigen::VectorXI> &leaf)
+  {
+    Eigen::VectorXN sqrD, S;
+    Eigen::VectorXI I;
+    Eigen::MatrixXN K;
+    igl::cycodebase::point_spline_squared_distance(
+      Q, P, C, B1, B2, leaf, sqrD, I, S, K);
+    return std::make_tuple(sqrD, I, S, K);
+  }
 }
 
 void bind_point_spline_squared_distance(nb::module_ &m)
@@ -32,6 +50,26 @@ void bind_point_spline_squared_distance(nb::module_ &m)
 @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 Bézier curves
+@return Tuple (sqrD, I, S, K) where
+  sqrD  #Q vector of smallest squared distances
+  I     #Q vector of indices of the closest cubic (row of C)
+  S     #Q vector of parameters of the closest points on that cubic
+  K     #Q by dim matrix of closest points on the spline)");
+
+  m.def("point_spline_squared_distance", &pyigl::point_spline_squared_distance_aabb,
+    "Q"_a, "P"_a, "C"_a, "B1"_a, "B2"_a, "leaf"_a,
+    R"(Squared distance to a spline, reusing a prebuilt Eytzinger AABB.
+
+Build the AABB once with igl.cycodebase.spline_eytzinger_aabb(P, C) and pass its
+(B1, B2, leaf) here so repeated queries against the same spline skip rebuilding
+the acceleration structure. The result matches the unaccelerated overload.
+
+@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 Bézier curves
+@param[in] B1    #B by dim matrix of AABB min box corners (from spline_eytzinger_aabb)
+@param[in] B2    #B by dim matrix of AABB max box corners (from spline_eytzinger_aabb)
+@param[in] leaf  #B vector of AABB leaf node indices/flags (from spline_eytzinger_aabb)
 @return Tuple (sqrD, I, S, K) where
   sqrD  #Q vector of smallest squared distances
   I     #Q vector of indices of the closest cubic (row of C)

tests/test_all.py

@@ -1715,6 +1715,37 @@ def test_cycodebase_spline_eytzinger_aabb_and_distance():
     assert np.isclose(sqrD[1], 0.0, atol=1e-12)
 
 
+def test_cycodebase_point_spline_squared_distance_aabb_reuse():
+    # Build the Eytzinger AABB once, then reuse it across query sets without
+    # rebuilding. Accelerated results must match the unaccelerated overload.
+    P, C = _unit_square_spline()
+    B1, B2, leaf = igl.cycodebase.spline_eytzinger_aabb(P, C)
+
+    Q = np.array([[0.5, 0.5], [0.5, 0.0], [0.2, 0.9], [1.1, 0.3]])
+    ref = igl.cycodebase.point_spline_squared_distance(Q, P, C)
+    fast = igl.cycodebase.point_spline_squared_distance(Q, P, C, B1, B2, leaf)
+    for a, b in zip(ref, fast):
+        assert np.allclose(a, b)
+
+    # Reuse the same tree for a different query set.
+    Q2 = np.array([[0.9, 0.9], [0.5, 0.5]])
+    sqrD2, I2, S2, K2 = igl.cycodebase.point_spline_squared_distance(
+        Q2, P, C, B1, B2, leaf)
+    assert sqrD2.shape == (2,)
+    assert np.isclose(sqrD2[1], 0.25, atol=1e-12)
+
+
+def test_cycodebase_point_cubic_squared_distance_bases_reuse():
+    # Precompute the monomial bases once, then reuse across query sets.
+    C = np.array([[0.0, 0.0], [1.0, 2.0], [2.0, -2.0], [3.0, 0.0]])
+    Q = np.array([[1.5, 0.0], [2.0, 0.5], [2.5, 1.0]])
+    ref = igl.cycodebase.point_cubic_squared_distance(Q, C)
+    M, D, B = igl.cubic_monomial_bases(C)
+    fast = igl.cycodebase.point_cubic_squared_distance(Q, C, D, B)
+    for a, b in zip(ref, fast):
+        assert np.allclose(a, b)
+
+
 # --------------------------------------------------------------------------
 # New predicates for cubic Bézier curves / splines
 # --------------------------------------------------------------------------