n3xd_ocp.measure: every face's area and centroid in one call

BRepGProp.SurfaceProperties is the hottest kernel call the app makes — 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.

face_surface_props(shape, *, parallel=True, eps=None) runs the whole scan
C++-side with the GIL released, optionally over OSD_Parallel, and returns
(areas[F], centroids[F,3]) in MapShapes(FACE) order — the face identity the
topology layer already keys on, so a caller indexes straight into it. A face
OCCT cannot integrate reports zeros, matching what the app's own try/except
substitutes.

Tested against the same stock-recorded per-face reference the Inc 1 gate uses,
and against the one-call-per-face loop it replaces. Parallel and serial are
compared with array_equal rather than approx: the parallel path shares one
TopoDS_Shape across threads, so an exact match is the evidence that nothing
reachable from it gets mutated while measuring.

Backend adoption comes after the cutover, so this changes nothing yet.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01DfriM8XUkn7uYf5Dwe2xo6

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"]

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

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();
 }

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 <nanobind/ndarray.h>
+
+#include <BRepGProp.hxx>
+#include <GProp_GProps.hxx>
+#include <OSD_Parallel.hxx>
+#include <TopAbs_ShapeEnum.hxx>
+#include <TopExp.hxx>
+#include <TopTools_IndexedMapOfShape.hxx>
+#include <TopoDS.hxx>
+#include <TopoDS_Face.hxx>
+#include <TopoDS_Shape.hxx>
+
+#include <vector>
+
+namespace {
+
+/// Own a heap buffer for as long as the numpy array that wraps it lives.
+template <typename T> 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<double>(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<nb::numpy, double, nb::ndim<1>>(areas, 1, area_shape,
+                                                            owner_of(areas)),
+                nb::ndarray<nb::numpy, double, nb::ndim<2>>(
+                    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");
+}

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"]

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)