From dc669ba84ccc5b7b9d6caff3e686464090cecb13 Mon Sep 17 00:00:00 2001 From: Horde Date: Mon, 10 Aug 2026 21:39:05 +0000 Subject: [PATCH] fast_winding_number: add point-cloud and pre-built tree support (#292) The binding previously only supported the (V,F,Q) triangle-soup one-shot. This adds the rest of libigl's fast_winding_number API so point clouds and reusable acceleration structures are available: - Point-cloud one-shot: fast_winding_number(P, N, A, Q, expansion_order=2, beta=2.0). - Reusable triangle-soup BVH: FastWindingNumberBVH class with init(V, F, order) and winding_number(Q, accuracy_scale) so many query sets reuse one build. - Point-cloud octree precompute: fast_winding_number_precompute(P, N, A, point_indices, CH, order) -> (CM, R, EC), plus a cached-evaluation overload fast_winding_number(P, N, A, point_indices, CH, CM, R, EC, Q, beta). Build the octree with the existing igl.octree(P). Adds test_fast_winding_number covering all four paths on a unit-sphere point set / mesh (inside ~1, outside ~0; cached paths match the one-shots). Full suite: 92 passed. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/fast_winding_number.cpp | 161 ++++++++++++++++++++++++++++++++++-- tests/test_all.py | 44 ++++++++++ 2 files changed, 196 insertions(+), 9 deletions(-) diff --git a/src/fast_winding_number.cpp b/src/fast_winding_number.cpp index 421bce2a..54fe3ccc 100644 --- a/src/fast_winding_number.cpp +++ b/src/fast_winding_number.cpp @@ -1,38 +1,181 @@ #include #include +#include +#include #include #include "default_types.h" +#include namespace nb = nanobind; using namespace nb::literals; namespace pyigl { - auto fast_winding_number( + // ---- Triangle soup, one-shot ------------------------------------------ + auto fast_winding_number_mesh( const nb::DRef &V, const nb::DRef &F, const nb::DRef &Q) { - Eigen::VectorXN W; // Winding number values for each query point + Eigen::VectorXN W; igl::fast_winding_number(V, F, Q, W); return W; } + + // ---- Point cloud, one-shot -------------------------------------------- + auto fast_winding_number_points( + const nb::DRef &P, + const nb::DRef &N, + const nb::DRef &A, + const nb::DRef &Q, + const int expansion_order, + const Numeric beta) + { + Eigen::VectorXN WN; + igl::fast_winding_number(P, N, A, Q, expansion_order, beta, WN); + return WN; + } + + // ---- Point cloud, cache the octree expansion -------------------------- + auto fast_winding_number_precompute( + const nb::DRef &P, + const nb::DRef &N, + const nb::DRef &A, + const std::vector> &point_indices, + const nb::DRef &CH, + const int expansion_order) + { + Eigen::MatrixXN CM; + Eigen::VectorXN R; + Eigen::MatrixXN EC; + igl::fast_winding_number( + P, N, A, point_indices, CH, expansion_order, CM, R, EC); + return std::make_tuple(CM, R, EC); + } + + // ---- Point cloud, evaluate using cached octree expansion -------------- + auto fast_winding_number_points_cached( + const nb::DRef &P, + const nb::DRef &N, + const nb::DRef &A, + const std::vector> &point_indices, + const nb::DRef &CH, + const nb::DRef &CM, + const nb::DRef &R, + const nb::DRef &EC, + const nb::DRef &Q, + const Numeric beta) + { + Eigen::VectorXN WN; + igl::fast_winding_number( + P, N, A, point_indices, CH, CM, R, EC, Q, beta, WN); + return WN; + } + + // ---- Triangle soup, cached BVH for repeated queries ------------------- + void fwn_bvh_init( + igl::FastWindingNumberBVH &bvh, + const nb::DRef &V, + const nb::DRef &F, + const int order) + { + igl::fast_winding_number(V, F, order, bvh); + } + + auto fwn_bvh_query( + const igl::FastWindingNumberBVH &bvh, + const nb::DRef &Q, + const Numeric accuracy_scale) + { + Eigen::VectorXN W; + igl::fast_winding_number(bvh, (float)accuracy_scale, Q, W); + return W; + } } -// Bind the wrapper to the Python module void bind_fast_winding_number(nb::module_ &m) { m.def( "fast_winding_number", - &pyigl::fast_winding_number, - "V"_a, - "F"_a, - "Q"_a, + &pyigl::fast_winding_number_mesh, + "V"_a, "F"_a, "Q"_a, R"(Compute approximate winding number for each query point based on a triangle soup mesh. @param[in] V #V by 3 matrix of mesh vertex positions @param[in] F #F by 3 matrix of triangle indices @param[in] Q #Q by 3 matrix of query positions -@return W #Q vector of winding number values for each query point)" - ); +@return W #Q vector of winding number values for each query point)"); + + m.def( + "fast_winding_number", + &pyigl::fast_winding_number_points, + "P"_a, "N"_a, "A"_a, "Q"_a, "expansion_order"_a = 2, "beta"_a = 2.0, +R"(Compute approximate winding number for each query point based on an oriented point cloud. + +@param[in] P #P by 3 matrix of point locations +@param[in] N #P by 3 matrix of point normals +@param[in] A #P vector of point areas +@param[in] Q #Q by 3 matrix of query positions +@param[in] expansion_order Taylor series expansion order (1, 2, or 3) +@param[in] beta Barnes-Hut accuracy term separating near from far field; + larger is more accurate and slower (2 is recommended). beta <= 0 forces exact + evaluation. +@return WN #Q vector of winding number values for each query point)"); + + m.def( + "fast_winding_number_precompute", + &pyigl::fast_winding_number_precompute, + "P"_a, "N"_a, "A"_a, "point_indices"_a, "CH"_a, "expansion_order"_a = 2, +R"(Precompute the octree Taylor expansion for point-cloud fast winding number, +so many query sets can reuse it. Build the octree with igl.octree(P). + +@param[in] P #P by 3 matrix of point locations +@param[in] N #P by 3 matrix of point normals +@param[in] A #P vector of point areas +@param[in] point_indices list of lists of point indices per octree cell (from igl.octree) +@param[in] CH #cells by 8 matrix of octree children (from igl.octree) +@param[in] expansion_order Taylor series expansion order (1, 2, or 3) +@return Tuple (CM, R, EC) where + CM #cells by 3 matrix of each cell's center of mass + R #cells vector of each cell's radius + EC #cells by #coefficients matrix of expansion coefficients)"); + + m.def( + "fast_winding_number", + &pyigl::fast_winding_number_points_cached, + "P"_a, "N"_a, "A"_a, "point_indices"_a, "CH"_a, "CM"_a, "R"_a, "EC"_a, + "Q"_a, "beta"_a = 2.0, +R"(Evaluate point-cloud fast winding number using a precomputed octree expansion +(see fast_winding_number_precompute). + +@param[in] P #P by 3 matrix of point locations +@param[in] N #P by 3 matrix of point normals +@param[in] A #P vector of point areas +@param[in] point_indices list of lists of point indices per octree cell +@param[in] CH #cells by 8 matrix of octree children +@param[in] CM #cells by 3 matrix of cell centers of mass +@param[in] R #cells vector of cell radii +@param[in] EC #cells by #coefficients matrix of expansion coefficients +@param[in] Q #Q by 3 matrix of query positions +@param[in] beta Barnes-Hut accuracy term (2 recommended; <= 0 forces exact) +@return WN #Q vector of winding number values for each query point)"); + + nb::class_(m, "FastWindingNumberBVH", +R"(Cached bounding-volume hierarchy for triangle-soup fast winding number, +enabling many query sets against the same mesh without rebuilding. + +Construct empty, then call init(V, F). Query with winding_number(Q).)") + .def(nb::init<>()) + .def("init", &pyigl::fwn_bvh_init, "V"_a, "F"_a, "order"_a = 2, +R"(Build the hierarchy for mesh (V, F). + +@param[in] V #V by 3 matrix of mesh vertex positions +@param[in] F #F by 3 matrix of triangle indices +@param[in] order Taylor series expansion order (e.g. 2))") + .def("winding_number", &pyigl::fwn_bvh_query, "Q"_a, "accuracy_scale"_a = 2.0, +R"(Compute winding numbers at query points using the cached hierarchy. + +@param[in] Q #Q by 3 matrix of query positions +@param[in] accuracy_scale accuracy parameter (e.g. 2) +@return W #Q vector of winding number values for each query point)"); } diff --git a/tests/test_all.py b/tests/test_all.py index 370c7431..0f6ce171 100644 --- a/tests/test_all.py +++ b/tests/test_all.py @@ -1791,3 +1791,47 @@ def test_resolve_duplicated_faces(): F = np.array([[0, 1, 2], [0, 1, 2], [0, 1, 2], [3, 4, 5]], dtype=np.int64) F2, J = igl.resolve_duplicated_faces(F) assert set(map(tuple, F2.tolist())) == {(3, 4, 5)} + + +def test_fast_winding_number(): + # Unit-sphere mesh from a subdivided icosahedron. + V, F = igl.icosahedron() + for _ in range(3): + V, F = igl.upsample(V, F) + V = V / np.linalg.norm(V, axis=1, keepdims=True) + + Q = np.array([[0.0, 0.0, 0.0], # inside -> ~1 + [0.5, 0.0, 0.0], # inside -> ~1 + [2.0, 0.0, 0.0], # outside -> ~0 + [0.0, 0.0, 3.0]]) # outside -> ~0 + + def inside_outside(W): + assert W.shape == (Q.shape[0],) + assert W[0] > 0.9 and W[1] > 0.9 + assert abs(W[2]) < 0.1 and abs(W[3]) < 0.1 + + # Mesh one-shot + Wm = igl.fast_winding_number(V, F, Q) + inside_outside(Wm) + + # Point cloud: vertices as oriented points with voronoi areas + N = V.copy() # outward unit normals on a unit sphere + A = np.asarray(igl.massmatrix(V, F, igl.MASSMATRIX_TYPE_VORONOI).diagonal()).ravel() + + # Point-cloud one-shot + Wp = igl.fast_winding_number(V, N, A, Q) + inside_outside(Wp) + + # Reusable triangle-soup BVH: same answer as the one-shot, reusable across + # multiple query sets without rebuilding. + bvh = igl.FastWindingNumberBVH() + bvh.init(V, F, 2) + Wb = bvh.winding_number(Q) + inside_outside(Wb) + assert np.allclose(Wb, bvh.winding_number(Q)) + + # Point-cloud octree precompute + cached evaluation matches the one-shot. + point_indices, CH, CN, W = igl.octree(V) + CM, R, EC = igl.fast_winding_number_precompute(V, N, A, point_indices, CH, 2) + Wc = igl.fast_winding_number(V, N, A, point_indices, CH, CM, R, EC, Q) + assert np.allclose(Wc, Wp, atol=1e-6)