ocp/tests/test_ext_tess.py

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