n3xd_ocp.tess: a shape's triangulation without the per-node round trip

cad/tessellation.py walks every node and every triangle of every face in
Python — three interpreter loops per face — to turn OCCT's mesh into flat
arrays the kernel already holds. extract_meshes hands them over in one call
with the GIL released; extract_edge_polylines does the same for the per-edge
GCPnts discretisation behind the viewport's edge overlay.

The output reproduces the Python loops exactly, including the parts that look
like quirks: unmeshed faces are omitted so face_index keeps its gaps (it is the
app-wide face identity, not a list position), triangles are rewound and normals
negated for reversed faces, and nodes come back float64 — the binary transport
narrows to float32 on its way out but the JSON one does not, so narrowing here
would quietly change what the REST payload says.

Plain dicts of arrays rather than a bound class, because the cad pool sends
meshes across a pipe and dicts of ndarrays pickle natively. There is a test for
that, and one asserting the arrays outlive the shape.

The equality tests re-implement the Python loops in full and compare with
array_equal, not approx: a difference here moves what the viewport draws or
what a face pick resolves to.

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

python/n3xd_ocp/__init__.py

@@ -17,6 +17,7 @@ import OCP as _OCP_pkg  # noqa: F401  (loads the extension, registering us)
 
 bintools = sys.modules["n3xd_ocp.bintools"]
 measure = sys.modules["n3xd_ocp.measure"]
+tess = sys.modules["n3xd_ocp.tess"]
 _debug = sys.modules["n3xd_ocp._debug"]
 
-__all__ = ["bintools", "measure"]
+__all__ = ["bintools", "measure", "tess"]

src/ext/ext_core.cpp

@@ -19,9 +19,11 @@ namespace nb = nanobind;
 void register_ext_bintools();
 void register_ext_debug();
 void register_ext_measure();
+void register_ext_tess();
 
 void register_ext(nb::module_ &) {
     register_ext_bintools();
     register_ext_debug();
     register_ext_measure();
+    register_ext_tess();
 }

src/ext/ext_tess.cpp

@@ -0,0 +1,222 @@
+/*
+    n3xd_ocp.tess — a shape's triangulation, without the per-node round trip.
+
+    cad/tessellation.py walks every node and every triangle of every face in
+    Python, three interpreter loops per face, to turn OCCT's mesh into flat
+    arrays.  The kernel already holds exactly those arrays; this hands them
+    over in one call, with the GIL released for the copy.
+
+    The output deliberately reproduces what the Python loops produce, down to
+    the details that look like quirks and are not:
+
+      * faces with no triangulation are **omitted**, leaving a gap in the id
+        sequence.  face_index is the MapShapes(FACE) ordinal, which is the face
+        identity the whole app keys on — picking, anchors, entity maps — so the
+        gap has to stay.
+      * triangles are 0-based and **rewound** for a reversed face, so the
+        winding a renderer sees is the outward one.
+      * normals are **negated** for a reversed face, for the same reason, and
+        are None when OCCT computed none (the frontend then derives its own).
+      * nodes come back float64.  The app's binary transport casts to float32
+        on its way out, but its JSON transport does not, so narrowing here
+        would quietly change what the REST payload says.
+
+    Edges get the same treatment: one call replaces a per-edge
+    GCPnts_TangentialDeflection plus a per-point Python loop.
+*/
+
+#include "../common/occt_module.h"
+
+#include <nanobind/ndarray.h>
+
+#include <BRepAdaptor_Curve.hxx>
+#include <BRep_Tool.hxx>
+#include <GCPnts_TangentialDeflection.hxx>
+#include <Poly_Triangulation.hxx>
+#include <TopAbs_ShapeEnum.hxx>
+#include <TopExp.hxx>
+#include <TopLoc_Location.hxx>
+#include <TopTools_IndexedMapOfShape.hxx>
+#include <TopoDS.hxx>
+#include <TopoDS_Edge.hxx>
+#include <TopoDS_Face.hxx>
+#include <TopoDS_Shape.hxx>
+#include <gp_Pnt.hxx>
+#include <gp_Trsf.hxx>
+
+#include <vector>
+
+namespace {
+
+template <typename T> nb::capsule owner_of(T *data) {
+    return nb::capsule(data, [](void *p) noexcept { delete[] (T *) p; });
+}
+
+/// Wrap an owned buffer as an (n, cols) numpy array.
+template <typename T>
+nb::object as_matrix(T *data, size_t rows, size_t cols) {
+    size_t shape[2] = {rows, cols};
+    return nb::cast(
+        nb::ndarray<nb::numpy, T, nb::ndim<2>>(data, 2, shape, owner_of(data)));
+}
+
+/// One face's mesh, already flattened and oriented, ready to hand to Python.
+struct FaceMesh {
+    int index = 0;
+    std::vector<double> nodes;      // 3N
+    std::vector<double> normals;    // 3N, empty when OCCT computed none
+    std::vector<uint32_t> tris;     // 3M
+};
+
+}  // namespace
+
+void register_ext_tess() {
+    nb::module_ m = ocp_named_module("n3xd_ocp.tess");
+
+    m.def(
+        "extract_meshes",
+        [](const TopoDS_Shape &shape, bool want_normals, bool apply_location,
+           bool flip_reversed) {
+            TopTools_IndexedMapOfShape faces;
+            TopExp::MapShapes(shape, TopAbs_FACE, faces);
+
+            std::vector<FaceMesh> meshes;
+            {
+                nb::gil_scoped_release nogil;
+                for (int i = 1; i <= faces.Extent(); ++i) {
+                    const TopoDS_Face &face = TopoDS::Face(faces.FindKey(i));
+                    TopLoc_Location loc;
+                    Handle(Poly_Triangulation) tri = BRep_Tool::Triangulation(face, loc);
+                    if (tri.IsNull())
+                        continue;  // keeps the id gap the app relies on
+
+                    const bool reversed =
+                        flip_reversed && face.Orientation() == TopAbs_REVERSED;
+                    const double sign = reversed ? -1.0 : 1.0;
+                    const gp_Trsf trsf = loc.Transformation();
+                    const bool identity = !apply_location || loc.IsIdentity();
+                    const bool normals =
+                        want_normals && tri->HasNormals();
+
+                    FaceMesh mesh;
+                    mesh.index = i - 1;
+                    mesh.nodes.resize(3 * (size_t) tri->NbNodes());
+                    if (normals)
+                        mesh.normals.resize(3 * (size_t) tri->NbNodes());
+                    mesh.tris.resize(3 * (size_t) tri->NbTriangles());
+
+                    for (int n = 1; n <= tri->NbNodes(); ++n) {
+                        gp_Pnt point = tri->Node(n);
+                        if (!identity)
+                            point.Transform(trsf);
+                        mesh.nodes[3 * (n - 1) + 0] = point.X();
+                        mesh.nodes[3 * (n - 1) + 1] = point.Y();
+                        mesh.nodes[3 * (n - 1) + 2] = point.Z();
+
+                        if (normals) {
+                            gp_Dir normal = tri->Normal(n);
+                            if (!identity)
+                                normal = normal.Transformed(trsf);
+                            mesh.normals[3 * (n - 1) + 0] = sign * normal.X();
+                            mesh.normals[3 * (n - 1) + 1] = sign * normal.Y();
+                            mesh.normals[3 * (n - 1) + 2] = sign * normal.Z();
+                        }
+                    }
+
+                    for (int t = 1; t <= tri->NbTriangles(); ++t) {
+                        Standard_Integer a = 0, b = 0, c = 0;
+                        tri->Triangle(t).Get(a, b, c);
+                        mesh.tris[3 * (t - 1) + 0] = (uint32_t) (a - 1);
+                        if (reversed) {
+                            mesh.tris[3 * (t - 1) + 1] = (uint32_t) (c - 1);
+                            mesh.tris[3 * (t - 1) + 2] = (uint32_t) (b - 1);
+                        } else {
+                            mesh.tris[3 * (t - 1) + 1] = (uint32_t) (b - 1);
+                            mesh.tris[3 * (t - 1) + 2] = (uint32_t) (c - 1);
+                        }
+                    }
+                    meshes.push_back(std::move(mesh));
+                }
+            }
+
+            nb::list out;
+            for (const FaceMesh &mesh : meshes) {
+                const size_t n = mesh.nodes.size() / 3;
+                const size_t t = mesh.tris.size() / 3;
+
+                double *nodes = new double[mesh.nodes.size()];
+                std::copy(mesh.nodes.begin(), mesh.nodes.end(), nodes);
+                uint32_t *tris = new uint32_t[mesh.tris.size()];
+                std::copy(mesh.tris.begin(), mesh.tris.end(), tris);
+
+                nb::dict entry;
+                entry["face_index"] = mesh.index;
+                entry["nodes"] = as_matrix(nodes, n, 3);
+                entry["triangles"] = as_matrix(tris, t, 3);
+                if (mesh.normals.empty()) {
+                    entry["normals"] = nb::none();
+                } else {
+                    double *normals = new double[mesh.normals.size()];
+                    std::copy(mesh.normals.begin(), mesh.normals.end(), normals);
+                    entry["normals"] = as_matrix(normals, n, 3);
+                }
+                out.append(entry);
+            }
+            return out;
+        },
+        "shape"_a, nb::kw_only(), "want_normals"_a = true,
+        "apply_location"_a = true, "flip_reversed"_a = true,
+        R"doc(Every triangulated face's mesh, as a list of dicts.
+
+Keys: ``face_index`` (0-based MapShapes(FACE) ordinal), ``nodes`` (float64
+[N,3]), ``triangles`` (uint32 [M,3], 0-based), ``normals`` (float64 [N,3] or
+None).  Faces OCCT has not meshed are omitted, so face_index has gaps — it is
+the app-wide face identity, not a position in this list.)doc");
+
+    m.def(
+        "extract_edge_polylines",
+        [](const TopoDS_Shape &shape, double angular, double linear) {
+            TopTools_IndexedMapOfShape edges;
+            TopExp::MapShapes(shape, TopAbs_EDGE, edges);
+
+            std::vector<std::pair<int, std::vector<double>>> polylines;
+            {
+                nb::gil_scoped_release nogil;
+                for (int i = 1; i <= edges.Extent(); ++i) {
+                    try {
+                        BRepAdaptor_Curve curve(TopoDS::Edge(edges.FindKey(i)));
+                        GCPnts_TangentialDeflection points(curve, angular, linear);
+                        std::vector<double> flat;
+                        flat.reserve(3 * (size_t) points.NbPoints());
+                        for (int p = 1; p <= points.NbPoints(); ++p) {
+                            const gp_Pnt point = points.Value(p);
+                            flat.push_back(point.X());
+                            flat.push_back(point.Y());
+                            flat.push_back(point.Z());
+                        }
+                        polylines.emplace_back(i - 1, std::move(flat));
+                    } catch (const Standard_Failure &) {
+                        // A curve that will not discretise is skipped, which
+                        // leaves the same id gap the Python path leaves.
+                    }
+                }
+            }
+
+            nb::list out;
+            for (const auto &entry : polylines) {
+                double *points = new double[entry.second.size()];
+                std::copy(entry.second.begin(), entry.second.end(), points);
+                nb::dict row;
+                row["edge_index"] = entry.first;
+                row["points"] = as_matrix(points, entry.second.size() / 3, 3);
+                out.append(row);
+            }
+            return out;
+        },
+        "shape"_a, "angular"_a, "linear"_a,
+        R"doc(Every edge discretised to a polyline.
+
+Returns dicts of ``edge_index`` (0-based MapShapes(EDGE) ordinal) and
+``points`` (float64 [P,3]).  An edge that will not discretise is skipped,
+leaving the same gap the per-edge Python loop leaves.)doc");
+}

tests/test_ext_tess.py

@@ -0,0 +1,197 @@
+"""n3xd_ocp.tess — the bulk mesh extractor.
+
+Every check here is against the Python loops in cad/tessellation.py, reproduced
+below.  The point of this API is to be indistinguishable from them, so the
+comparisons are exact: any difference would move what the viewport draws or
+what a face pick resolves to.
+"""
+
+from __future__ import annotations
+
+import pytest
+
+np = pytest.importorskip("numpy")
+
+import n3xd_ocp
+from OCP.BRep import BRep_Tool
+from OCP.BRepAdaptor import BRepAdaptor_Curve
+from OCP.BRepMesh import BRepMesh_IncrementalMesh
+from OCP.BRepPrimAPI import BRepPrimAPI_MakeBox, BRepPrimAPI_MakeCylinder
+from OCP.GCPnts import GCPnts_TangentialDeflection
+from OCP.gp import gp_Ax2, gp_Dir, gp_Pnt
+from OCP.TopAbs import TopAbs_EDGE, TopAbs_FACE, TopAbs_REVERSED
+from OCP.TopLoc import TopLoc_Location
+from OCP.TopoDS import TopoDS
+
+from .test_inc1_modeling import sub_shapes
+
+LINEAR, ANGULAR = 0.05, 0.15
+
+
+@pytest.fixture(scope="module")
+def meshed_cylinder():
+    """A cylinder rather than a box: it has a curved face with real normals,
+    a seam, and enough triangles for a winding mistake to show."""
+    shape = BRepPrimAPI_MakeCylinder(
+        gp_Ax2(gp_Pnt(0, 0, 0), gp_Dir(0, 0, 1)), 5.0, 12.0
+    ).Shape()
+    BRepMesh_IncrementalMesh(shape, LINEAR, False, ANGULAR, True)
+    return shape
+
+
+def _python_faces(shape):
+    """cad/tessellation.py::_extract_faces, reduced to what it produces."""
+    out = []
+    for i, raw in enumerate(sub_shapes(shape, TopAbs_FACE)):
+        face = TopoDS.Face_s(raw)
+        loc = TopLoc_Location()
+        tri = BRep_Tool.Triangulation_s(face, loc)
+        if tri is None:
+            continue  # the id gap
+        trsf = loc.Transformation()
+        identity = loc.IsIdentity()
+        reversed_ = face.Orientation() == TopAbs_REVERSED
+        sign = -1.0 if reversed_ else 1.0
+
+        nodes = []
+        for n in range(1, tri.NbNodes() + 1):
+            node = tri.Node(n)
+            if not identity:
+                node.Transform(trsf)
+            nodes.extend((node.X(), node.Y(), node.Z()))
+
+        normals = []
+        if tri.HasNormals():
+            for n in range(1, tri.NbNodes() + 1):
+                normal = tri.Normal(n)
+                if not identity:
+                    normal = normal.Transformed(trsf)
+                normals.extend((sign * normal.X(), sign * normal.Y(),
+                                sign * normal.Z()))
+
+        indices = []
+        for t in range(1, tri.NbTriangles() + 1):
+            a, b, c = tri.Triangle(t).Get()
+            if reversed_:
+                indices.extend((a - 1, c - 1, b - 1))
+            else:
+                indices.extend((a - 1, b - 1, c - 1))
+
+        out.append(
+            {"id": i, "vertices": nodes, "normals": normals, "indices": indices}
+        )
+    return out
+
+
+def test_matches_the_python_extraction_exactly(meshed_cylinder):
+    got = n3xd_ocp.tess.extract_meshes(meshed_cylinder)
+    expected = _python_faces(meshed_cylinder)
+
+    assert [f["face_index"] for f in got] == [f["id"] for f in expected]
+    for mine, theirs in zip(got, expected, strict=True):
+        assert np.array_equal(mine["nodes"].ravel(), np.array(theirs["vertices"]))
+        assert np.array_equal(
+            mine["triangles"].ravel(), np.array(theirs["indices"], dtype="uint32")
+        )
+        if theirs["normals"]:
+            assert np.array_equal(
+                mine["normals"].ravel(), np.array(theirs["normals"])
+            )
+        else:
+            assert mine["normals"] is None
+
+
+def test_face_index_is_the_map_ordinal_with_gaps_preserved():
+    """A shape meshed only in part: the untriangulated faces must be absent
+    from the list and must not shift the ids of the ones that remain."""
+    shape = BRepPrimAPI_MakeBox(10.0, 20.0, 30.0).Shape()
+    assert n3xd_ocp.tess.extract_meshes(shape) == []  # nothing meshed yet
+
+    BRepMesh_IncrementalMesh(shape, LINEAR, False, ANGULAR, True)
+    got = n3xd_ocp.tess.extract_meshes(shape)
+    assert [f["face_index"] for f in got] == list(range(6))
+
+
+def test_reversed_faces_are_rewound_and_their_normals_flipped(meshed_cylinder):
+    """Turning the flip off has to change something, or the test above would
+    pass for the wrong reason."""
+    flipped = n3xd_ocp.tess.extract_meshes(meshed_cylinder, flip_reversed=True)
+    raw = n3xd_ocp.tess.extract_meshes(meshed_cylinder, flip_reversed=False)
+
+    faces = [TopoDS.Face_s(f) for f in sub_shapes(meshed_cylinder, TopAbs_FACE)]
+    reversed_ids = {
+        i for i, f in enumerate(faces) if f.Orientation() == TopAbs_REVERSED
+    }
+    assert reversed_ids, "fixture has no reversed face, so this proves nothing"
+
+    for a, b in zip(flipped, raw, strict=True):
+        if a["face_index"] in reversed_ids:
+            assert not np.array_equal(a["triangles"], b["triangles"])
+            # BRepMesh does not compute normals unless asked, so this half of
+            # the check only fires when the mesh actually carries them.
+            if a["normals"] is not None:
+                assert np.array_equal(a["normals"], -b["normals"])
+        else:
+            assert np.array_equal(a["triangles"], b["triangles"])
+
+
+def test_normals_can_be_skipped(meshed_cylinder):
+    for face in n3xd_ocp.tess.extract_meshes(meshed_cylinder, want_normals=False):
+        assert face["normals"] is None
+
+
+def test_location_is_applied(meshed_cylinder):
+    """A located shape's nodes must arrive in world space, since that is what
+    the viewport draws."""
+    from OCP.gp import gp_Trsf, gp_Vec
+
+    trsf = gp_Trsf()
+    trsf.SetTranslation(gp_Vec(100.0, 0.0, 0.0))
+    moved = meshed_cylinder.Moved(TopLoc_Location(trsf))
+
+    placed = n3xd_ocp.tess.extract_meshes(moved, apply_location=True)
+    local = n3xd_ocp.tess.extract_meshes(moved, apply_location=False)
+    assert placed[0]["nodes"][:, 0].min() == pytest.approx(
+        local[0]["nodes"][:, 0].min() + 100.0
+    )
+
+
+def test_arrays_outlive_the_shape(meshed_cylinder):
+    """Nothing in the result points into OCCT's own storage."""
+    import gc
+
+    faces = n3xd_ocp.tess.extract_meshes(meshed_cylinder)
+    total = float(faces[0]["nodes"].sum())
+    gc.collect()
+    assert float(faces[0]["nodes"].sum()) == total
+
+
+def test_edge_polylines_match_the_python_loop(meshed_cylinder):
+    got = n3xd_ocp.tess.extract_edge_polylines(meshed_cylinder, ANGULAR, LINEAR)
+
+    expected = []
+    for i, raw in enumerate(sub_shapes(meshed_cylinder, TopAbs_EDGE)):
+        curve = BRepAdaptor_Curve(TopoDS.Edge_s(raw))
+        points = GCPnts_TangentialDeflection(curve, ANGULAR, LINEAR)
+        flat = []
+        for p in range(1, points.NbPoints() + 1):
+            point = points.Value(p)
+            flat.extend((point.X(), point.Y(), point.Z()))
+        expected.append((i, flat))
+
+    assert [e["edge_index"] for e in got] == [i for i, _ in expected]
+    for mine, (_i, flat) in zip(got, expected, strict=True):
+        assert np.array_equal(mine["points"].ravel(), np.array(flat))
+
+
+def test_results_pickle(meshed_cylinder):
+    """The cad pool sends meshes across a pipe, so the dicts have to pickle —
+    which is the reason these are plain dicts of arrays and not a bound class.
+    """
+    import pickle
+
+    faces = n3xd_ocp.tess.extract_meshes(meshed_cylinder)
+    restored = pickle.loads(pickle.dumps(faces))
+    assert len(restored) == len(faces)
+    assert np.array_equal(restored[0]["nodes"], faces[0]["nodes"])
+    assert restored[0]["face_index"] == faces[0]["face_index"]