diff --git a/python/n3xd_ocp/__init__.py b/python/n3xd_ocp/__init__.py index d6a3ade..da8d817 100644 --- a/python/n3xd_ocp/__init__.py +++ b/python/n3xd_ocp/__init__.py @@ -16,6 +16,7 @@ import sys import OCP as _OCP_pkg # noqa: F401 (loads the extension, registering us) bintools = sys.modules["n3xd_ocp.bintools"] +measure = sys.modules["n3xd_ocp.measure"] _debug = sys.modules["n3xd_ocp._debug"] -__all__ = ["bintools"] +__all__ = ["bintools", "measure"] diff --git a/scripts/_env.sh b/scripts/_env.sh index 4d6f097..4bfbd4e 100755 --- a/scripts/_env.sh +++ b/scripts/_env.sh @@ -22,7 +22,9 @@ if [ ! -x "$VENV/bin/python" ]; then echo "--- creating build venv at $VENV ---" /opt/python/cp312-cp312/bin/python -m venv "$VENV" "$VENV/bin/pip" install -q --upgrade pip - "$VENV/bin/pip" install -q build nanobind scikit-build-core pytest + # numpy is a test dependency only: the n3xd_ocp bulk APIs hand back + # ndarrays, so their tests need it to check the results. + "$VENV/bin/pip" install -q build nanobind scikit-build-core pytest numpy fi PY="$VENV/bin/python" export PY VENV CACHE diff --git a/src/ext/ext_core.cpp b/src/ext/ext_core.cpp index eb6b843..914b243 100644 --- a/src/ext/ext_core.cpp +++ b/src/ext/ext_core.cpp @@ -18,8 +18,10 @@ namespace nb = nanobind; void register_ext_bintools(); void register_ext_debug(); +void register_ext_measure(); void register_ext(nb::module_ &) { register_ext_bintools(); register_ext_debug(); + register_ext_measure(); } diff --git a/src/ext/ext_measure.cpp b/src/ext/ext_measure.cpp new file mode 100644 index 0000000..7b66dc4 --- /dev/null +++ b/src/ext/ext_measure.cpp @@ -0,0 +1,113 @@ +/* + n3xd_ocp.measure — every face's area and centroid in one call. + + BRepGProp.SurfaceProperties is the single hottest kernel call the app + makes: profiled on a thread-heavy benchmark it was 94 % of + face_candidate_anchors, 0.99 s of 1.05 s for 690 faces. Not because the + kernel is slow, but because it is reached once per face from Python, so a + rebuild pays the round-trip hundreds of times per feature over a growing + shape. + + Here the whole scan runs C++-side with the GIL released, optionally spread + over OSD_Parallel. Results come back in TopExp::MapShapes(FACE) order, + which is the face identity the entire topology layer keys on, so a caller + can index straight into it. + + A face OCCT cannot integrate yields area 0 and a zero centroid rather than + raising, matching what the app's own try/except substitutes today. +*/ + +#include "../common/occt_module.h" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace { + +/// Own a heap buffer for as long as the numpy array that wraps it lives. +template nb::capsule owner_of(T *data) { + return nb::capsule(data, [](void *p) noexcept { delete[] (T *) p; }); +} + +void measure_range(const TopTools_IndexedMapOfShape &faces, double eps, + double *areas, double *centroids, int begin, int end) { + for (int i = begin; i < end; ++i) { + double area = 0.0, cx = 0.0, cy = 0.0, cz = 0.0; + try { + GProp_GProps props; + const TopoDS_Face &face = TopoDS::Face(faces.FindKey(i + 1)); + if (eps > 0.0) + BRepGProp::SurfaceProperties(face, props, eps); + else + BRepGProp::SurfaceProperties(face, props); + area = props.Mass(); + const gp_Pnt centre = props.CentreOfMass(); + cx = centre.X(); + cy = centre.Y(); + cz = centre.Z(); + } catch (const Standard_Failure &) { + // Leave the zeros: the app substitutes 0.0 for an unmeasurable + // face rather than failing the scan. + } + areas[i] = area; + centroids[3 * i + 0] = cx; + centroids[3 * i + 1] = cy; + centroids[3 * i + 2] = cz; + } +} + +} // namespace + +void register_ext_measure() { + nb::module_ m = ocp_named_module("n3xd_ocp.measure"); + + m.def( + "face_surface_props", + [](const TopoDS_Shape &shape, bool parallel, nb::object eps_obj) { + const double eps = eps_obj.is_none() ? 0.0 : nb::cast(eps_obj); + + TopTools_IndexedMapOfShape faces; + TopExp::MapShapes(shape, TopAbs_FACE, faces); + const int count = faces.Extent(); + + double *areas = new double[count > 0 ? count : 1]; + double *centroids = new double[count > 0 ? 3 * count : 1]; + + { + nb::gil_scoped_release nogil; + if (parallel && count > 1) { + OSD_Parallel::For(0, count, [&](int i) { + measure_range(faces, eps, areas, centroids, i, i + 1); + }); + } else { + measure_range(faces, eps, areas, centroids, 0, count); + } + } + + size_t area_shape[1] = {(size_t) count}; + size_t centroid_shape[2] = {(size_t) count, 3}; + return nb::make_tuple( + nb::ndarray>(areas, 1, area_shape, + owner_of(areas)), + nb::ndarray>( + centroids, 2, centroid_shape, owner_of(centroids))); + }, + "shape"_a, nb::kw_only(), "parallel"_a = true, "eps"_a = nb::none(), + R"doc(Area and centre of mass of every face, in MapShapes(FACE) order. + +Returns ``(areas[F], centroids[F, 3])`` as float64 arrays. ``eps`` selects +BRepGProp's tolerance-driven quadrature, which is more accurate on +spline-bounded faces; leave it None for OCCT's default. A face that cannot be +integrated reports area 0 and a zero centroid.)doc"); +} diff --git a/tests/conftest.py b/tests/conftest.py index 1bb5b55..c68c37e 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -49,3 +49,11 @@ def a_face(fixture_shapes): faces = TopTools_IndexedMapOfShape() TopExp.MapShapes_s(fixture_shapes["box_meshed"], TopAbs_FACE, faces) return faces.FindKey(1) + + +@pytest.fixture(scope="session") +def inc1_or_skip(manifest): + """The Inc 1 reference block, for tests outside test_inc1_modeling.py.""" + if "inc1" not in manifest: + pytest.skip("manifest predates the Inc 1 reference block") + return manifest["inc1"] diff --git a/tests/test_ext_measure.py b/tests/test_ext_measure.py new file mode 100644 index 0000000..cb5a823 --- /dev/null +++ b/tests/test_ext_measure.py @@ -0,0 +1,102 @@ +"""n3xd_ocp.measure — the batched face measurement. + +Checked against the same per-face reference the Inc 1 gate uses, which is +recorded from the stock wheel: the batched scan has to agree with calling +BRepGProp.SurfaceProperties once per face, in the same MapShapes(FACE) order, +or the anchor path would resolve against different numbers. +""" + +from __future__ import annotations + +import pytest + +np = pytest.importorskip("numpy") + +import n3xd_ocp +from OCP.BRepGProp import BRepGProp +from OCP.BRepPrimAPI import BRepPrimAPI_MakeBox +from OCP.GProp import GProp_GProps +from OCP.TopAbs import TopAbs_FACE + +from .test_inc1_modeling import sub_shapes + + +def _per_face(shape, eps=None): + """The one-call-per-face loop this API exists to replace.""" + rows = [] + for face in sub_shapes(shape, TopAbs_FACE): + props = GProp_GProps() + if eps is None: + BRepGProp.SurfaceProperties_s(face, props) + else: + BRepGProp.SurfaceProperties_s(face, props, eps) + centre = props.CentreOfMass() + rows.append((props.Mass(), (centre.X(), centre.Y(), centre.Z()))) + return rows + + +def test_matches_the_stock_reference(inc1_or_skip, fixture_shapes): + for name, expected in inc1_or_skip["face_props"].items(): + areas, centroids = n3xd_ocp.measure.face_surface_props(fixture_shapes[name]) + assert areas.shape == (len(expected),), name + assert centroids.shape == (len(expected), 3), name + for i, row in enumerate(expected): + assert areas[i] == pytest.approx(row[0], rel=1e-9), f"{name}[{i}].area" + assert centroids[i] == pytest.approx(row[1:], rel=1e-9, abs=1e-12), ( + f"{name}[{i}].centroid" + ) + + +def test_matches_the_per_face_loop(fixture_shapes): + shape = fixture_shapes["cut_cylinder"] + areas, centroids = n3xd_ocp.measure.face_surface_props(shape) + for i, (area, centre) in enumerate(_per_face(shape)): + assert areas[i] == pytest.approx(area, rel=1e-12) + assert centroids[i] == pytest.approx(centre, rel=1e-12, abs=1e-12) + + +def test_parallel_and_serial_agree_exactly(fixture_shapes): + """The parallel path shares one TopoDS_Shape across threads. If OCCT + mutated anything reachable from it while measuring, this is where it would + show — so the comparison is exact, not approximate.""" + shape = fixture_shapes["fused"] + par = n3xd_ocp.measure.face_surface_props(shape, parallel=True) + ser = n3xd_ocp.measure.face_surface_props(shape, parallel=False) + assert np.array_equal(par[0], ser[0]) + assert np.array_equal(par[1], ser[1]) + + +def test_eps_overload_is_reachable(fixture_shapes): + shape = fixture_shapes["box"] + areas, _ = n3xd_ocp.measure.face_surface_props(shape, eps=1e-6) + for i, (area, _centre) in enumerate(_per_face(shape, eps=1e-6)): + assert areas[i] == pytest.approx(area, rel=1e-12) + + +def test_empty_shape_returns_empty_arrays(fixture_shapes): + areas, centroids = n3xd_ocp.measure.face_surface_props( + fixture_shapes["empty_compound"] + ) + assert areas.shape == (0,) + assert centroids.shape == (0, 3) + + +def test_arrays_own_their_memory(fixture_shapes): + """The buffers are allocated C++-side and handed over with a capsule + deleter; the arrays must stay valid once every other reference is gone.""" + import gc + + areas, centroids = n3xd_ocp.measure.face_surface_props(fixture_shapes["box"]) + expected = float(areas.sum()) + gc.collect() + assert float(areas.sum()) == expected + assert centroids.base is not None # the capsule, not the shape + + +def test_survives_the_shape_being_dropped(): + """Nothing in the result points back into the shape.""" + areas, centroids = n3xd_ocp.measure.face_surface_props( + BRepPrimAPI_MakeBox(10.0, 20.0, 30.0).Shape() + ) + assert areas.sum() == pytest.approx(2200.0) + assert centroids.shape == (6, 3)