fast_winding_number: point-cloud and pre-built tree support (#292)

src/fast_winding_number.cpp

@@ -1,38 +1,181 @@
 #include <nanobind/nanobind.h>
 #include <nanobind/eigen/dense.h>
+#include <nanobind/stl/tuple.h>
+#include <nanobind/stl/vector.h>
 #include <igl/fast_winding_number.h>
 #include "default_types.h"
+#include <vector>
 
 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<const Eigen::MatrixXN> &V,
     const nb::DRef<const Eigen::MatrixXI> &F,
     const nb::DRef<const Eigen::MatrixXN> &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<const Eigen::MatrixXN> &P,
+    const nb::DRef<const Eigen::MatrixXN> &N,
+    const nb::DRef<const Eigen::VectorXN> &A,
+    const nb::DRef<const Eigen::MatrixXN> &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<const Eigen::MatrixXN> &P,
+    const nb::DRef<const Eigen::MatrixXN> &N,
+    const nb::DRef<const Eigen::VectorXN> &A,
+    const std::vector<std::vector<Integer>> &point_indices,
+    const nb::DRef<const Eigen::MatrixXI> &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<const Eigen::MatrixXN> &P,
+    const nb::DRef<const Eigen::MatrixXN> &N,
+    const nb::DRef<const Eigen::VectorXN> &A,
+    const std::vector<std::vector<Integer>> &point_indices,
+    const nb::DRef<const Eigen::MatrixXI> &CH,
+    const nb::DRef<const Eigen::MatrixXN> &CM,
+    const nb::DRef<const Eigen::VectorXN> &R,
+    const nb::DRef<const Eigen::MatrixXN> &EC,
+    const nb::DRef<const Eigen::MatrixXN> &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<const Eigen::MatrixXN> &V,
+    const nb::DRef<const Eigen::MatrixXI> &F,
+    const int order)
+  {
+    igl::fast_winding_number(V, F, order, bvh);
+  }
+
+  auto fwn_bvh_query(
+    const igl::FastWindingNumberBVH &bvh,
+    const nb::DRef<const Eigen::MatrixXN> &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_<igl::FastWindingNumberBVH>(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)");
 }

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)