diff --git a/python/n3xd_ocp/__init__.py b/python/n3xd_ocp/__init__.py index da8d817..5a24113 100644 --- a/python/n3xd_ocp/__init__.py +++ b/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"] diff --git a/src/ext/ext_core.cpp b/src/ext/ext_core.cpp index 914b243..e2a98f8 100644 --- a/src/ext/ext_core.cpp +++ b/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(); } diff --git a/src/ext/ext_tess.cpp b/src/ext/ext_tess.cpp new file mode 100644 index 0000000..bf1010d --- /dev/null +++ b/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 + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace { + +template 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 +nb::object as_matrix(T *data, size_t rows, size_t cols) { + size_t shape[2] = {rows, cols}; + return nb::cast( + nb::ndarray>(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 nodes; // 3N + std::vector normals; // 3N, empty when OCCT computed none + std::vector 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 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>> 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 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"); +} diff --git a/tests/test_ext_tess.py b/tests/test_ext_tess.py new file mode 100644 index 0000000..b7586c6 --- /dev/null +++ b/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"]