From c112d9506d303ddf97213162ac486249002cb052 Mon Sep 17 00:00:00 2001 From: artem-ogre Date: Tue, 8 Sep 2026 22:13:19 -0500 Subject: [PATCH 1/8] Add triangles_as_array() and release the GIL Co-authored-by: Chris Colbert --- cdt_bindings.cpp | 77 +++++++++++++++++++++++++++++++--------- cdt_bindings_test.py | 83 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 143 insertions(+), 17 deletions(-) diff --git a/cdt_bindings.cpp b/cdt_bindings.cpp index 692dc4b..7329430 100644 --- a/cdt_bindings.cpp +++ b/cdt_bindings.cpp @@ -11,6 +11,7 @@ #include #include +#include #include #include #include @@ -204,6 +205,29 @@ PYBIND11_MODULE(PythonCDT, m) t.triangles.begin(), t.triangles.end()); }, py::keep_alive<0, 1>()) + // numpy read-back: the whole triangle vector in ONE memcpy, as a + // structured array over the dtype registered above. Reading + // t.triangles instead materialises a Python object per triangle, + // which on a large mesh costs more than the triangulation itself. + // Contents are as stored: the super-triangle's three vertices (indices + // 0-2) and the triangles touching them are present until an erase_* + // call finalizes the triangulation. A copy, so the array outlives + // any later change to the triangulation. + .def( + "triangles_as_array", + [](const Triangulation& t) { + py::array_t out( + static_cast(t.triangles.size())); + std::memcpy( + out.mutable_data(), + t.triangles.data(), + t.triangles.size() * sizeof(CDT::Triangle)); + return out; + }, + "All triangles as a structured numpy array of shape (T,) with " + "uint32 fields 'vertices' (3,) and 'neighbors' (3,); " + "arr['vertices'] is the (T, 3) vertex-index array. Includes the " + "super-triangle until erase_super_triangle()") // fixed edges .def_readonly("fixed_edges", &Triangulation::fixedEdges) .def( @@ -272,11 +296,16 @@ PYBIND11_MODULE(PythonCDT, m) }; const std::size_t n_vert = info.size / 2; const XY* const ptr = static_cast(info.ptr); - t.insertVertices( - ptr, - ptr + n_vert, - [](const XY& v) { return v.xy[0]; }, - [](const XY& v) { return v.xy[1]; }); + { + // The buffer is pinned by `info`; nothing below touches + // Python, so other threads may run meanwhile. + py::gil_scoped_release release; + t.insertVertices( + ptr, + ptr + n_vert, + [](const XY& v) { return v.xy[0]; }, + [](const XY& v) { return v.xy[1]; }); + } }, py::arg("vertex_buffer")) .def( @@ -311,11 +340,16 @@ PYBIND11_MODULE(PythonCDT, m) }; const std::size_t n_vert = info.size / 2; const EdgeData* const ptr = static_cast(info.ptr); - t.insertEdges( - ptr, - ptr + n_vert, - [](const EdgeData& e) { return e.vv[0]; }, - [](const EdgeData& e) { return e.vv[1]; }); + { + // The buffer is pinned by `info`; nothing below touches + // Python, so other threads may run meanwhile. + py::gil_scoped_release release; + t.insertEdges( + ptr, + ptr + n_vert, + [](const EdgeData& e) { return e.vv[0]; }, + [](const EdgeData& e) { return e.vv[1]; }); + } }, py::arg("edge_buffer")) .def( @@ -350,14 +384,22 @@ PYBIND11_MODULE(PythonCDT, m) }; const std::size_t n_vert = info.size / 2; const EdgeData* const ptr = static_cast(info.ptr); - t.conformToEdges( - ptr, - ptr + n_vert, - [](const EdgeData& e) { return e.vv[0]; }, - [](const EdgeData& e) { return e.vv[1]; }); + { + // The buffer is pinned by `info`; nothing below touches + // Python, so other threads may run meanwhile. + py::gil_scoped_release release; + t.conformToEdges( + ptr, + ptr + n_vert, + [](const EdgeData& e) { return e.vv[0]; }, + [](const EdgeData& e) { return e.vv[1]; }); + } }, py::arg("edge_buffer")) - .def("erase_super_triangle", &Triangulation::eraseSuperTriangle) + .def( + "erase_super_triangle", + &Triangulation::eraseSuperTriangle, + py::call_guard()) .def("erase_outer_triangles", &Triangulation::eraseOuterTriangles) .def( "erase_outer_triangles_and_holes", @@ -375,5 +417,6 @@ PYBIND11_MODULE(PythonCDT, m) m.def( "verify_topology", &CDT::verifyTopology, - py::arg("triangulation")); + py::arg("triangulation"), + py::call_guard()); } diff --git a/cdt_bindings_test.py b/cdt_bindings_test.py index c98d4e9..1835dd7 100644 --- a/cdt_bindings_test.py +++ b/cdt_bindings_test.py @@ -10,10 +10,23 @@ import pytest import tempfile import hashlib +import threading import PythonCDT as cdt +def assert_triangles_array_matches_objects(t: cdt.Triangulation) -> None: + """triangles_as_array() must agree exactly with the object API, super-triangle included.""" + arr = t.triangles_as_array() + assert arr.shape == (len(t.triangles),) + assert arr.dtype.names == ("vertices", "neighbors") + assert arr["vertices"].shape == (len(t.triangles), 3) and arr["vertices"].dtype == np.uintc + assert arr["neighbors"].shape == (len(t.triangles), 3) and arr["neighbors"].dtype == np.uintc + for i, tri in enumerate(t.triangles): + assert list(tri.vertices) == arr["vertices"][i].tolist(), "triangles_as_array vertices disagree with triangles" + assert list(tri.neighbors) == arr["neighbors"][i].tolist(), "triangles_as_array neighbors disagree with triangles" + + def test_constants() -> None: """Test that constants have proper values""" assert cdt.NO_NEIGHBOR == np.iinfo(np.uintc).max, "NO_NEIGHBOR constant has wrong value" @@ -55,6 +68,7 @@ def test_Triangulation() -> None: assert len(t.vertices) == 0, "Wrong vertex count in empty triangulation" assert len(t.triangles) == 0, "Wrong triangle count in empty triangulation" assert len(t.fixed_edges) == 0, "Wrong fixed edge count in empty triangulation" + assert_triangles_array_matches_objects(t) vv = [cdt.V2d(-1, 0), cdt.V2d(0, 0.5), cdt.V2d(1, 0), cdt.V2d(0, -0.5)] t.insert_vertices(vv) @@ -66,12 +80,14 @@ def test_Triangulation() -> None: t.insert_edges(ee) assert len(t.fixed_edges) == 1, "Wrong fixed edge count in triangulation" assert cdt.Edge(0 + 3, 2 + 3) in t.fixed_edges, "Constraint edge was not properly added" + assert_triangles_array_matches_objects(t) t.erase_super_triangle() assert cdt.Edge(0, 2) in t.fixed_edges, "Constraint edge was not properly added" assert len(t.vertices) == 4, "Wrong vertex count in triangulation" assert len(t.triangles) == 2, "Wrong triangle count in triangulation" assert len(t.fixed_edges) == 1, "Wrong fixed edge count in triangulation" + assert_triangles_array_matches_objects(t) # test retrieving triangulation data using iterators assert t.vertices_count() == len(t.vertices), "Wrong vertex count" @@ -98,6 +114,7 @@ def test_Triangulation() -> None: t.erase_super_triangle() assert len(t.vertices) == 5, "Wrong vertex count in triangulation" assert len(t.triangles) == 4, "Wrong triangle count in triangulation" + assert_triangles_array_matches_objects(t) def test_verify_topology() -> None: @@ -106,6 +123,7 @@ def test_verify_topology() -> None: t.insert_vertices([cdt.V2d(-1, 0), cdt.V2d(0, 0.5), cdt.V2d(1, 0), cdt.V2d(0, -0.5)]) t.insert_edges([cdt.Edge(0, 2), cdt.Edge(1, 3)]) assert cdt.verify_topology(t), "Verifying topology produced wrong result" + assert_triangles_array_matches_objects(t) def save_triangulation_as_off(t: cdt.Triangulation, off_file) -> None: @@ -146,6 +164,7 @@ def test_triangulate_input_file() -> None: off_file = f"{tmp_dir}/cdt.off" save_triangulation_as_off(t, off_file) assert md5_checksum(off_file) == 'db59c00d9dad866781cd96779e5262b7', "Wrong OFF file contents" + assert_triangles_array_matches_objects(t) def test_conform_to_edges() -> None: @@ -155,6 +174,7 @@ def test_conform_to_edges() -> None: t.conform_to_edges(ee) t.erase_outer_triangles_and_holes() assert triangulation_md5_checksum(t) == 'b64cae39c91a55dd4e23a146eb7df0d3', "Wrong OFF file contents" + assert_triangles_array_matches_objects(t) @pytest.mark.parametrize("vv", [[cdt.V2d(-1, 0), cdt.V2d(0, 0.5), cdt.V2d(1, 0), cdt.V2d(0, -0.5)], @@ -167,6 +187,7 @@ def test_insert_vertices(vv) -> None: assert len(t.triangles) == 9, "Wrong triangle count in triangulation" assert len(t.fixed_edges) == 0, "Wrong fixed edge count in triangulation" assert triangulation_md5_checksum(t) == 'db9176f4429942862a7a73155fb55322', "Wrong OFF file contents" + assert_triangles_array_matches_objects(t) @pytest.mark.parametrize("ee", [[cdt.Edge(0, 1), cdt.Edge(2, 3), cdt.Edge(3, 4), cdt.Edge(5, 6)], @@ -181,6 +202,7 @@ def test_insert_conform_edges(ee) -> None: assert len(t.triangles) == 15, "Wrong triangle count in triangulation" assert len(t.fixed_edges) == 4, "Wrong fixed edge count in triangulation" assert triangulation_md5_checksum(t) == '639c7a1492b2adb8f25464ec81ff6a00', "Wrong OFF file contents" + assert_triangles_array_matches_objects(t) # conform to edges t = cdt.Triangulation(cdt.VertexInsertionOrder.AS_PROVIDED, cdt.IntersectingConstraintEdges.NOT_ALLOWED, 0.0) @@ -190,3 +212,64 @@ def test_insert_conform_edges(ee) -> None: assert len(t.triangles) == 19, "Wrong triangle count in triangulation" assert len(t.fixed_edges) == 6, "Wrong fixed edge count in triangulation" assert triangulation_md5_checksum(t) == '9c87b435e247c1658ec0f04af3340dc7', "Wrong OFF file contents" + assert_triangles_array_matches_objects(t) + + +def test_triangles_as_array_before_and_after_erase() -> None: + """Before erase_super_triangle() the array holds the super-triangle's fan and input vertex i + is index i + 3; after it, the fan is gone and input vertex i is index i.""" + vv, ee = read_input_file("CDT/visualizer/data/ditch.txt") + t = cdt.Triangulation(cdt.VertexInsertionOrder.AS_PROVIDED, cdt.IntersectingConstraintEdges.TRY_RESOLVE, 0.0) + t.insert_vertices(vv) + t.insert_edges(ee) + + before = t.triangles_as_array()["vertices"] + assert (before < 3).any(axis=1).sum() > 0, "unfinalized array must include triangles touching the super-triangle" + kept = before[(before >= 3).all(axis=1)] - 3 + + t.erase_super_triangle() + after = t.triangles_as_array()["vertices"] + assert np.array_equal(np.sort(np.sort(kept, axis=1), axis=0), np.sort(np.sort(after, axis=1), axis=0)), \ + "dropping the super-triangle's fan and shifting by 3 must give the erased mesh" + assert after.max() < len(vv), "after erase every index addresses an input vertex" + + +def test_triangles_as_array_empty() -> None: + t = cdt.Triangulation(cdt.VertexInsertionOrder.AS_PROVIDED, cdt.IntersectingConstraintEdges.NOT_ALLOWED, 0.0) + arr = t.triangles_as_array() + assert arr.shape == (0,) and arr["vertices"].shape == (0, 3) + + +def test_triangles_as_array_is_a_copy() -> None: + """The array is a copy: writing to it must not alter the triangulation.""" + t = cdt.Triangulation(cdt.VertexInsertionOrder.AS_PROVIDED, cdt.IntersectingConstraintEdges.NOT_ALLOWED, 0.0) + t.insert_vertices(np.array([[-1, 0], [0, 0.5], [1, 0], [0, -0.5]], dtype=np.float64)) + arr = t.triangles_as_array() + arr["vertices"][:] = 0 + assert all(list(tri.vertices) != [0, 0, 0] for tri in t.triangles), "triangles_as_array must return a copy" + + +def test_insert_releases_the_gil() -> None: + """Concurrent builds on threads must each equal the single-threaded build.""" + vv, ee = read_input_file("CDT/visualizer/data/Constrained Sweden.txt") + verts = np.array([[v.x, v.y] for v in vv], dtype=np.float64) + edges = np.array([[e.v1, e.v2] for e in ee], dtype=np.uintc) + + def build() -> str: + t = cdt.Triangulation(cdt.VertexInsertionOrder.AS_PROVIDED, cdt.IntersectingConstraintEdges.TRY_RESOLVE, 0.0) + t.insert_vertices(verts) + t.insert_edges(edges) + t.erase_super_triangle() + assert cdt.verify_topology(t) + return triangulation_md5_checksum(t) + + expected = build() + results = [None] * 4 + def worker(i): + results[i] = build() + threads = [threading.Thread(target=worker, args=(i,)) for i in range(4)] + for th in threads: + th.start() + for th in threads: + th.join() + assert results == [expected] * 4, "Threaded builds must match the single-threaded build" From fb884d9f37845b4a40fabfcb836bd76d39755c9c Mon Sep 17 00:00:00 2001 From: artem-ogre Date: Wed, 9 Sep 2026 12:55:24 -0500 Subject: [PATCH 2/8] Clean up formatting Co-authored-by: Chris Colbert --- cdt_bindings.cpp | 23 ++++++----------------- 1 file changed, 6 insertions(+), 17 deletions(-) diff --git a/cdt_bindings.cpp b/cdt_bindings.cpp index 7329430..718ee57 100644 --- a/cdt_bindings.cpp +++ b/cdt_bindings.cpp @@ -205,14 +205,6 @@ PYBIND11_MODULE(PythonCDT, m) t.triangles.begin(), t.triangles.end()); }, py::keep_alive<0, 1>()) - // numpy read-back: the whole triangle vector in ONE memcpy, as a - // structured array over the dtype registered above. Reading - // t.triangles instead materialises a Python object per triangle, - // which on a large mesh costs more than the triangulation itself. - // Contents are as stored: the super-triangle's three vertices (indices - // 0-2) and the triangles touching them are present until an erase_* - // call finalizes the triangulation. A copy, so the array outlives - // any later change to the triangulation. .def( "triangles_as_array", [](const Triangulation& t) { @@ -297,10 +289,9 @@ PYBIND11_MODULE(PythonCDT, m) const std::size_t n_vert = info.size / 2; const XY* const ptr = static_cast(info.ptr); { - // The buffer is pinned by `info`; nothing below touches - // Python, so other threads may run meanwhile. + // The buffer is pinned and nothing below touches Python. py::gil_scoped_release release; - t.insertVertices( + t.insertVertices( ptr, ptr + n_vert, [](const XY& v) { return v.xy[0]; }, @@ -341,10 +332,9 @@ PYBIND11_MODULE(PythonCDT, m) const std::size_t n_vert = info.size / 2; const EdgeData* const ptr = static_cast(info.ptr); { - // The buffer is pinned by `info`; nothing below touches - // Python, so other threads may run meanwhile. + // The buffer is pinned and nothing below touches Python. py::gil_scoped_release release; - t.insertEdges( + t.insertEdges( ptr, ptr + n_vert, [](const EdgeData& e) { return e.vv[0]; }, @@ -385,10 +375,9 @@ PYBIND11_MODULE(PythonCDT, m) const std::size_t n_vert = info.size / 2; const EdgeData* const ptr = static_cast(info.ptr); { - // The buffer is pinned by `info`; nothing below touches - // Python, so other threads may run meanwhile. + // The buffer is pinned and nothing below touches Python. py::gil_scoped_release release; - t.conformToEdges( + t.conformToEdges( ptr, ptr + n_vert, [](const EdgeData& e) { return e.vv[0]; }, From a009408c4603b1378a9d0ded8d1bf83c5266d574 Mon Sep 17 00:00:00 2001 From: artem-ogre Date: Wed, 9 Sep 2026 13:04:03 -0500 Subject: [PATCH 3/8] Release the GIL in more places Co-authored-by: Chris Colbert --- cdt_bindings.cpp | 26 ++++++++++++++++---------- 1 file changed, 16 insertions(+), 10 deletions(-) diff --git a/cdt_bindings.cpp b/cdt_bindings.cpp index 718ee57..78dc112 100644 --- a/cdt_bindings.cpp +++ b/cdt_bindings.cpp @@ -261,7 +261,8 @@ PYBIND11_MODULE(PythonCDT, m) "insert_vertices", static_cast&)>( &Triangulation::insertVertices), - py::arg("vertices")) + py::arg("vertices"), + py::call_guard()) .def( "insert_vertices", [](Triangulation& t, py::buffer b) { @@ -289,7 +290,6 @@ PYBIND11_MODULE(PythonCDT, m) const std::size_t n_vert = info.size / 2; const XY* const ptr = static_cast(info.ptr); { - // The buffer is pinned and nothing below touches Python. py::gil_scoped_release release; t.insertVertices( ptr, @@ -303,7 +303,8 @@ PYBIND11_MODULE(PythonCDT, m) "insert_edges", static_cast&)>( &Triangulation::insertEdges), - py::arg("edges")) + py::arg("edges"), + py::call_guard()) .def( "insert_edges", [](Triangulation& t, py::buffer b) { @@ -332,7 +333,6 @@ PYBIND11_MODULE(PythonCDT, m) const std::size_t n_vert = info.size / 2; const EdgeData* const ptr = static_cast(info.ptr); { - // The buffer is pinned and nothing below touches Python. py::gil_scoped_release release; t.insertEdges( ptr, @@ -346,7 +346,8 @@ PYBIND11_MODULE(PythonCDT, m) "conform_to_edges", static_cast&)>( &Triangulation::conformToEdges), - py::arg("edges")) + py::arg("edges"), + py::call_guard()) .def( "conform_to_edges", [](Triangulation& t, py::buffer b) { @@ -375,7 +376,6 @@ PYBIND11_MODULE(PythonCDT, m) const std::size_t n_vert = info.size / 2; const EdgeData* const ptr = static_cast(info.ptr); { - // The buffer is pinned and nothing below touches Python. py::gil_scoped_release release; t.conformToEdges( ptr, @@ -389,19 +389,25 @@ PYBIND11_MODULE(PythonCDT, m) "erase_super_triangle", &Triangulation::eraseSuperTriangle, py::call_guard()) - .def("erase_outer_triangles", &Triangulation::eraseOuterTriangles) + .def( + "erase_outer_triangles", + &Triangulation::eraseOuterTriangles, + py::call_guard()) .def( "erase_outer_triangles_and_holes", - &Triangulation::eraseOuterTrianglesAndHoles) + &Triangulation::eraseOuterTrianglesAndHoles, + py::call_guard()) .def("is_finalized", &Triangulation::isFinalized) .def( "calculate_triangle_depths", - &Triangulation::calculateTriangleDepths) + &Triangulation::calculateTriangleDepths, + py::call_guard()) .def( "remove_triangles", static_cast( &Triangulation::removeTriangles), - py::arg("triangle_indices")); + py::arg("triangle_indices"), + py::call_guard()); m.def( "verify_topology", From df1f43c220f04ba0857271f990c2237a78a833d9 Mon Sep 17 00:00:00 2001 From: artem-ogre Date: Mon, 14 Sep 2026 23:58:07 +0200 Subject: [PATCH 4/8] Validate shape and contiguity of input buffers --- cdt_bindings.cpp | 132 ++++++++++++++++++------------------------- cdt_bindings_test.py | 25 ++++++++ 2 files changed, 81 insertions(+), 76 deletions(-) diff --git a/cdt_bindings.cpp b/cdt_bindings.cpp index 78dc112..86c5766 100644 --- a/cdt_bindings.cpp +++ b/cdt_bindings.cpp @@ -29,6 +29,48 @@ using Triangulation = CDT::Triangulation; namespace { +template +struct BufferPair +{ + T v[2]; +}; + +// Reads a C-contiguous buffer of shape (2N,) or (N, 2) as N pairs +template +std::pair*, std::size_t> buffer_pairs( + const py::buffer_info& info, + const std::string& type_name, + const std::string& item_name) +{ + if (info.format != py::format_descriptor::format()) + { + throw std::runtime_error( + "Incompatible format: expected a " + type_name + " array!"); + } + if (info.ndim != 1 && info.ndim != 2) + { + throw std::runtime_error("Incompatible buffer dimension!"); + } + if (info.ndim == 2 ? info.shape[1] != 2 : info.shape[0] % 2 != 0) + { + throw std::runtime_error( + "Buffer must hold " + type_name + " pairs (2 per " + item_name + + "): shape (N, 2) or (2N,)!"); + } + // strides of dimensions with a single element are irrelevant + py::ssize_t expected_stride = info.itemsize; + for (py::ssize_t d = info.ndim - 1; d >= 0; --d) + { + if (info.shape[d] > 1 && info.strides[d] != expected_stride) + { + throw std::runtime_error( + "Buffer must be C-contiguous: use numpy.ascontiguousarray!"); + } + expected_stride *= info.shape[d]; + } + return {static_cast*>(info.ptr), info.size / 2}; +} + std::string TriInd2str(CDT::TriInd it) { return it != CDT::noNeighbor ? std::to_string(it) : "-"; @@ -267,35 +309,15 @@ PYBIND11_MODULE(PythonCDT, m) "insert_vertices", [](Triangulation& t, py::buffer b) { const py::buffer_info info = b.request(); - // sanity checks - if (info.format != py::format_descriptor::format()) - { - throw std::runtime_error( - "Incompatible format: expected a double array!"); - } - if (info.ndim != 1 && info.ndim != 2) - { - throw std::runtime_error("Incompatible buffer dimension!"); - } - if (info.size % 2 != 0) - { - throw std::runtime_error("Buffer must hold even number of " - "coordinates (2 per vertex)!"); - } - // create from buffer - struct XY - { - coord_t xy[2]; - }; - const std::size_t n_vert = info.size / 2; - const XY* const ptr = static_cast(info.ptr); + const auto [ptr, n_vert] = + buffer_pairs(info, "double", "vertex"); { py::gil_scoped_release release; t.insertVertices( ptr, ptr + n_vert, - [](const XY& v) { return v.xy[0]; }, - [](const XY& v) { return v.xy[1]; }); + [](const BufferPair& v) { return v.v[0]; }, + [](const BufferPair& v) { return v.v[1]; }); } }, py::arg("vertex_buffer")) @@ -309,36 +331,15 @@ PYBIND11_MODULE(PythonCDT, m) "insert_edges", [](Triangulation& t, py::buffer b) { const py::buffer_info info = b.request(); - // sanity checks - if (info.format != - py::format_descriptor::format()) - { - throw std::runtime_error( - "Incompatible format: expected a CDT::VertInd array!"); - } - if (info.ndim != 1 && info.ndim != 2) - { - throw std::runtime_error("Incompatible buffer dimension!"); - } - if (info.size % 2 != 0) - { - throw std::runtime_error("Buffer must hold even number of " - "CDT::VertInd (2 per edge)!"); - } - // create from buffer - struct EdgeData - { - CDT::VertInd vv[2]; - }; - const std::size_t n_vert = info.size / 2; - const EdgeData* const ptr = static_cast(info.ptr); + const auto [ptr, n_edges] = + buffer_pairs(info, "CDT::VertInd", "edge"); { py::gil_scoped_release release; t.insertEdges( ptr, - ptr + n_vert, - [](const EdgeData& e) { return e.vv[0]; }, - [](const EdgeData& e) { return e.vv[1]; }); + ptr + n_edges, + [](const BufferPair& e) { return e.v[0]; }, + [](const BufferPair& e) { return e.v[1]; }); } }, py::arg("edge_buffer")) @@ -352,36 +353,15 @@ PYBIND11_MODULE(PythonCDT, m) "conform_to_edges", [](Triangulation& t, py::buffer b) { const py::buffer_info info = b.request(); - // sanity checks - if (info.format != - py::format_descriptor::format()) - { - throw std::runtime_error( - "Incompatible format: expected a CDT::VertInd array!"); - } - if (info.ndim != 1 && info.ndim != 2) - { - throw std::runtime_error("Incompatible buffer dimension!"); - } - if (info.size % 2 != 0) - { - throw std::runtime_error("Buffer must hold even number of " - "CDT::VertInd (2 per edge)!"); - } - // create from buffer - struct EdgeData - { - CDT::VertInd vv[2]; - }; - const std::size_t n_vert = info.size / 2; - const EdgeData* const ptr = static_cast(info.ptr); + const auto [ptr, n_edges] = + buffer_pairs(info, "CDT::VertInd", "edge"); { py::gil_scoped_release release; t.conformToEdges( ptr, - ptr + n_vert, - [](const EdgeData& e) { return e.vv[0]; }, - [](const EdgeData& e) { return e.vv[1]; }); + ptr + n_edges, + [](const BufferPair& e) { return e.v[0]; }, + [](const BufferPair& e) { return e.v[1]; }); } }, py::arg("edge_buffer")) diff --git a/cdt_bindings_test.py b/cdt_bindings_test.py index 1835dd7..d235e27 100644 --- a/cdt_bindings_test.py +++ b/cdt_bindings_test.py @@ -273,3 +273,28 @@ def worker(i): for th in threads: th.join() assert results == [expected] * 4, "Threaded builds must match the single-threaded build" + + +def test_insert_buffers_must_be_contiguous_pairs() -> None: + """Buffer overloads reject buffers that aren't a C-contiguous run of pairs""" + vertices = np.array([[0, 0], [9, 9], [1, 0], [9, 9], [0, 1], [9, 9], [1, 1], [9, 9]], dtype=np.float64) + edges = np.array([[0, 1], [9, 9], [1, 3], [9, 9]], dtype=np.uintc) + + def triangulation(): + return cdt.Triangulation(cdt.VertexInsertionOrder.AS_PROVIDED, cdt.IntersectingConstraintEdges.NOT_ALLOWED, 0.0) + + for bad in [vertices[::2], vertices[:, ::-1], np.asfortranarray(vertices[::2]), np.zeros((2, 3))]: + with pytest.raises(RuntimeError): + triangulation().insert_vertices(bad) + for insert in [cdt.Triangulation.insert_edges, cdt.Triangulation.conform_to_edges]: + t = triangulation() + t.insert_vertices(vertices[::2].copy()) + with pytest.raises(RuntimeError): + insert(t, edges[::2]) + + for good in [vertices[::2].copy(), vertices[::2].ravel()]: + t = triangulation() + t.insert_vertices(good) + assert [(v.x, v.y) for v in t.vertices][3:] == [(0, 0), (1, 0), (0, 1), (1, 1)], "Wrong vertices inserted" + t.insert_edges(edges[::2].copy()) + assert t.fixed_edges == {cdt.Edge(3, 4), cdt.Edge(4, 6)}, "Wrong edges inserted" From 2e8b758916eccdb0d0e5adeb58d480b7064d84b5 Mon Sep 17 00:00:00 2001 From: artem-ogre Date: Mon, 14 Sep 2026 23:58:49 +0200 Subject: [PATCH 5/8] Replace triangles_as_array() with vertices_array() and triangles_array() --- README.md | 19 +++++++++ benchmark_list_vs_array_access.py | 17 ++++++++ cdt_bindings.cpp | 47 +++++++++++++-------- cdt_bindings_test.py | 68 ++++++++----------------------- 4 files changed, 82 insertions(+), 69 deletions(-) create mode 100644 benchmark_list_vs_array_access.py diff --git a/README.md b/README.md index e874810..34a586f 100644 --- a/README.md +++ b/README.md @@ -16,6 +16,25 @@ pip3 install . pytest ./cdt_bindings_test.py ``` +## Usage + +```python +import numpy as np +import PythonCDT as cdt + +t = cdt.Triangulation(cdt.VertexInsertionOrder.AUTO, cdt.IntersectingConstraintEdges.TRY_RESOLVE, 0.0) +t.insert_vertices(np.array([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]])) +t.insert_edges(np.array([[0, 3]], dtype=np.uintc)) +t.erase_super_triangle() + +vertices = t.vertices_array() # numpy array with fields 'x' and 'y' +triangles = t.triangles_array() # numpy array with fields 'vertices' and 'neighbors' +triangles["vertices"] # (T, 3) vertex indices into vertices +``` + +`vertices_array()` and `triangles_array()` return copies. With `copy=False` they return read-only views of the +triangulation's memory without copying; a view is invalidated by any call that modifies the triangulation. + ## License [Mozilla Public License, v. 2.0](https://www.mozilla.org/en-US/MPL/2.0/FAQ/) diff --git a/benchmark_list_vs_array_access.py b/benchmark_list_vs_array_access.py new file mode 100644 index 0000000..aaad84a --- /dev/null +++ b/benchmark_list_vs_array_access.py @@ -0,0 +1,17 @@ +"""Compares reading a triangulation as Python lists vs numpy arrays""" + +import timeit + +import numpy as np + +import PythonCDT as cdt + +t = cdt.Triangulation(cdt.VertexInsertionOrder.AUTO, cdt.IntersectingConstraintEdges.NOT_ALLOWED, 0.0) +t.insert_vertices(np.random.default_rng(0).random((200_000, 2))) + +print(f"{'':<10}{'count':>9}{'list':>12}{'copy':>12}{'view':>12}") +for name in ["vertices", "triangles"]: + array = getattr(t, f"{name}_array") + ways = [lambda: getattr(t, name), array, lambda: array(copy=False)] + ms = [1e3 * min(timeit.repeat(way, number=1, repeat=5)) for way in ways] + print(f"{name:<10}{len(array()):>9}" + "".join(f"{m:>9.3f} ms" for m in ms)) diff --git a/cdt_bindings.cpp b/cdt_bindings.cpp index 86c5766..6e8b10b 100644 --- a/cdt_bindings.cpp +++ b/cdt_bindings.cpp @@ -11,7 +11,6 @@ #include #include -#include #include #include #include @@ -29,6 +28,32 @@ using Triangulation = CDT::Triangulation; namespace { +// Binds a vector member as `name_array(*, copy=True)` returning a numpy array +template +void def_array( + py::class_& cls, + const std::string& name, + std::vector Triangulation::*member) +{ + cls.def( + (name + "_array").c_str(), + [member](const Triangulation& t, bool copy) { + const std::vector& items = t.*member; + const auto size = static_cast(items.size()); + if (copy) + return py::array_t(size, items.data()); + // the view's base is the Python wrapper of t, keeping t alive + py::array_t view(size, items.data(), py::cast(&t)); + view.attr("setflags")(py::arg("write") = false); + return view; + }, + py::kw_only(), + py::arg("copy") = true, + "Copy as a numpy structured array. copy=False returns a read-only view " + "of the triangulation's memory instead; it is invalidated by any call " + "that modifies the triangulation."); +} + template struct BufferPair { @@ -215,7 +240,8 @@ PYBIND11_MODULE(PythonCDT, m) return oss.str(); }); - py::class_(m, "Triangulation") + py::class_ triangulation(m, "Triangulation"); + triangulation .def( py::init< CDT::VertexInsertionOrder::Enum, @@ -247,21 +273,6 @@ PYBIND11_MODULE(PythonCDT, m) t.triangles.begin(), t.triangles.end()); }, py::keep_alive<0, 1>()) - .def( - "triangles_as_array", - [](const Triangulation& t) { - py::array_t out( - static_cast(t.triangles.size())); - std::memcpy( - out.mutable_data(), - t.triangles.data(), - t.triangles.size() * sizeof(CDT::Triangle)); - return out; - }, - "All triangles as a structured numpy array of shape (T,) with " - "uint32 fields 'vertices' (3,) and 'neighbors' (3,); " - "arr['vertices'] is the (T, 3) vertex-index array. Includes the " - "super-triangle until erase_super_triangle()") // fixed edges .def_readonly("fixed_edges", &Triangulation::fixedEdges) .def( @@ -388,6 +399,8 @@ PYBIND11_MODULE(PythonCDT, m) &Triangulation::removeTriangles), py::arg("triangle_indices"), py::call_guard()); + def_array(triangulation, "vertices", &Triangulation::vertices); + def_array(triangulation, "triangles", &Triangulation::triangles); m.def( "verify_topology", diff --git a/cdt_bindings_test.py b/cdt_bindings_test.py index d235e27..9f335d3 100644 --- a/cdt_bindings_test.py +++ b/cdt_bindings_test.py @@ -15,18 +15,6 @@ import PythonCDT as cdt -def assert_triangles_array_matches_objects(t: cdt.Triangulation) -> None: - """triangles_as_array() must agree exactly with the object API, super-triangle included.""" - arr = t.triangles_as_array() - assert arr.shape == (len(t.triangles),) - assert arr.dtype.names == ("vertices", "neighbors") - assert arr["vertices"].shape == (len(t.triangles), 3) and arr["vertices"].dtype == np.uintc - assert arr["neighbors"].shape == (len(t.triangles), 3) and arr["neighbors"].dtype == np.uintc - for i, tri in enumerate(t.triangles): - assert list(tri.vertices) == arr["vertices"][i].tolist(), "triangles_as_array vertices disagree with triangles" - assert list(tri.neighbors) == arr["neighbors"][i].tolist(), "triangles_as_array neighbors disagree with triangles" - - def test_constants() -> None: """Test that constants have proper values""" assert cdt.NO_NEIGHBOR == np.iinfo(np.uintc).max, "NO_NEIGHBOR constant has wrong value" @@ -68,7 +56,6 @@ def test_Triangulation() -> None: assert len(t.vertices) == 0, "Wrong vertex count in empty triangulation" assert len(t.triangles) == 0, "Wrong triangle count in empty triangulation" assert len(t.fixed_edges) == 0, "Wrong fixed edge count in empty triangulation" - assert_triangles_array_matches_objects(t) vv = [cdt.V2d(-1, 0), cdt.V2d(0, 0.5), cdt.V2d(1, 0), cdt.V2d(0, -0.5)] t.insert_vertices(vv) @@ -80,14 +67,12 @@ def test_Triangulation() -> None: t.insert_edges(ee) assert len(t.fixed_edges) == 1, "Wrong fixed edge count in triangulation" assert cdt.Edge(0 + 3, 2 + 3) in t.fixed_edges, "Constraint edge was not properly added" - assert_triangles_array_matches_objects(t) t.erase_super_triangle() assert cdt.Edge(0, 2) in t.fixed_edges, "Constraint edge was not properly added" assert len(t.vertices) == 4, "Wrong vertex count in triangulation" assert len(t.triangles) == 2, "Wrong triangle count in triangulation" assert len(t.fixed_edges) == 1, "Wrong fixed edge count in triangulation" - assert_triangles_array_matches_objects(t) # test retrieving triangulation data using iterators assert t.vertices_count() == len(t.vertices), "Wrong vertex count" @@ -114,7 +99,6 @@ def test_Triangulation() -> None: t.erase_super_triangle() assert len(t.vertices) == 5, "Wrong vertex count in triangulation" assert len(t.triangles) == 4, "Wrong triangle count in triangulation" - assert_triangles_array_matches_objects(t) def test_verify_topology() -> None: @@ -123,7 +107,6 @@ def test_verify_topology() -> None: t.insert_vertices([cdt.V2d(-1, 0), cdt.V2d(0, 0.5), cdt.V2d(1, 0), cdt.V2d(0, -0.5)]) t.insert_edges([cdt.Edge(0, 2), cdt.Edge(1, 3)]) assert cdt.verify_topology(t), "Verifying topology produced wrong result" - assert_triangles_array_matches_objects(t) def save_triangulation_as_off(t: cdt.Triangulation, off_file) -> None: @@ -164,7 +147,6 @@ def test_triangulate_input_file() -> None: off_file = f"{tmp_dir}/cdt.off" save_triangulation_as_off(t, off_file) assert md5_checksum(off_file) == 'db59c00d9dad866781cd96779e5262b7', "Wrong OFF file contents" - assert_triangles_array_matches_objects(t) def test_conform_to_edges() -> None: @@ -174,7 +156,6 @@ def test_conform_to_edges() -> None: t.conform_to_edges(ee) t.erase_outer_triangles_and_holes() assert triangulation_md5_checksum(t) == 'b64cae39c91a55dd4e23a146eb7df0d3', "Wrong OFF file contents" - assert_triangles_array_matches_objects(t) @pytest.mark.parametrize("vv", [[cdt.V2d(-1, 0), cdt.V2d(0, 0.5), cdt.V2d(1, 0), cdt.V2d(0, -0.5)], @@ -187,7 +168,6 @@ def test_insert_vertices(vv) -> None: assert len(t.triangles) == 9, "Wrong triangle count in triangulation" assert len(t.fixed_edges) == 0, "Wrong fixed edge count in triangulation" assert triangulation_md5_checksum(t) == 'db9176f4429942862a7a73155fb55322', "Wrong OFF file contents" - assert_triangles_array_matches_objects(t) @pytest.mark.parametrize("ee", [[cdt.Edge(0, 1), cdt.Edge(2, 3), cdt.Edge(3, 4), cdt.Edge(5, 6)], @@ -202,7 +182,6 @@ def test_insert_conform_edges(ee) -> None: assert len(t.triangles) == 15, "Wrong triangle count in triangulation" assert len(t.fixed_edges) == 4, "Wrong fixed edge count in triangulation" assert triangulation_md5_checksum(t) == '639c7a1492b2adb8f25464ec81ff6a00', "Wrong OFF file contents" - assert_triangles_array_matches_objects(t) # conform to edges t = cdt.Triangulation(cdt.VertexInsertionOrder.AS_PROVIDED, cdt.IntersectingConstraintEdges.NOT_ALLOWED, 0.0) @@ -212,41 +191,26 @@ def test_insert_conform_edges(ee) -> None: assert len(t.triangles) == 19, "Wrong triangle count in triangulation" assert len(t.fixed_edges) == 6, "Wrong fixed edge count in triangulation" assert triangulation_md5_checksum(t) == '9c87b435e247c1658ec0f04af3340dc7', "Wrong OFF file contents" - assert_triangles_array_matches_objects(t) - - -def test_triangles_as_array_before_and_after_erase() -> None: - """Before erase_super_triangle() the array holds the super-triangle's fan and input vertex i - is index i + 3; after it, the fan is gone and input vertex i is index i.""" - vv, ee = read_input_file("CDT/visualizer/data/ditch.txt") - t = cdt.Triangulation(cdt.VertexInsertionOrder.AS_PROVIDED, cdt.IntersectingConstraintEdges.TRY_RESOLVE, 0.0) - t.insert_vertices(vv) - t.insert_edges(ee) - - before = t.triangles_as_array()["vertices"] - assert (before < 3).any(axis=1).sum() > 0, "unfinalized array must include triangles touching the super-triangle" - kept = before[(before >= 3).all(axis=1)] - 3 - - t.erase_super_triangle() - after = t.triangles_as_array()["vertices"] - assert np.array_equal(np.sort(np.sort(kept, axis=1), axis=0), np.sort(np.sort(after, axis=1), axis=0)), \ - "dropping the super-triangle's fan and shifting by 3 must give the erased mesh" - assert after.max() < len(vv), "after erase every index addresses an input vertex" -def test_triangles_as_array_empty() -> None: +@pytest.mark.parametrize("copy", [True, False]) +def test_arrays(copy) -> None: t = cdt.Triangulation(cdt.VertexInsertionOrder.AS_PROVIDED, cdt.IntersectingConstraintEdges.NOT_ALLOWED, 0.0) - arr = t.triangles_as_array() - assert arr.shape == (0,) and arr["vertices"].shape == (0, 3) + assert t.vertices_array(copy=copy).shape == t.triangles_array(copy=copy).shape == (0,), \ + "Empty triangulation must give empty arrays" - -def test_triangles_as_array_is_a_copy() -> None: - """The array is a copy: writing to it must not alter the triangulation.""" - t = cdt.Triangulation(cdt.VertexInsertionOrder.AS_PROVIDED, cdt.IntersectingConstraintEdges.NOT_ALLOWED, 0.0) - t.insert_vertices(np.array([[-1, 0], [0, 0.5], [1, 0], [0, -0.5]], dtype=np.float64)) - arr = t.triangles_as_array() - arr["vertices"][:] = 0 - assert all(list(tri.vertices) != [0, 0, 0] for tri in t.triangles), "triangles_as_array must return a copy" + t.insert_vertices([cdt.V2d(-1, 0), cdt.V2d(0, 0.5), cdt.V2d(1, 0), cdt.V2d(0, -0.5)]) + vertices = t.vertices_array(copy=copy) + triangles = t.triangles_array(copy=copy) + assert vertices.tolist() == [(v.x, v.y) for v in t.vertices] + assert triangles["vertices"].tolist() == [list(tri.vertices) for tri in t.triangles] + assert triangles["neighbors"].tolist() == [list(tri.neighbors) for tri in t.triangles] + + expected = [vertices.copy(), triangles.copy()] + del t + for arr, before in zip([vertices, triangles], expected): + assert arr.flags.owndata == arr.flags.writeable == copy, "Copy must be owned and writeable, view neither" + assert np.array_equal(arr, before), "Array must stay valid after the triangulation is deleted" def test_insert_releases_the_gil() -> None: From 8bb42b1c1e7c4085fccd918689e35a9ef378b9da Mon Sep 17 00:00:00 2001 From: artem-ogre Date: Mon, 14 Sep 2026 23:59:18 +0200 Subject: [PATCH 6/8] Make Triangulation safe to share between threads --- README.md | 6 + cdt_bindings.cpp | 285 ++++++++++++++++++++++++------------------- cdt_bindings_test.py | 45 +++---- 3 files changed, 186 insertions(+), 150 deletions(-) diff --git a/README.md b/README.md index 34a586f..b741c5e 100644 --- a/README.md +++ b/README.md @@ -35,6 +35,12 @@ triangles["vertices"] # (T, 3) vertex indices into vertices `vertices_array()` and `triangles_array()` return copies. With `copy=False` they return read-only views of the triangulation's memory without copying; a view is invalidated by any call that modifies the triangulation. +### Threads + +Triangulating releases the GIL, so separate triangulations can be built in parallel on Python threads. +One triangulation can be shared between threads: calls on it wait for each other. +Iterators (`*_iter()`) and `copy=False` views are not protected: don't use them while another thread modifies the triangulation. + ## License [Mozilla Public License, v. 2.0](https://www.mozilla.org/en-US/MPL/2.0/FAQ/) diff --git a/cdt_bindings.cpp b/cdt_bindings.cpp index 6e8b10b..53933ef 100644 --- a/cdt_bindings.cpp +++ b/cdt_bindings.cpp @@ -11,6 +11,7 @@ #include #include +#include #include #include #include @@ -23,22 +24,81 @@ namespace py = pybind11; using coord_t = double; using NearPointLocator_t = CDT::LocatorKDTree; using V2d = CDT::V2d; -using Triangulation = CDT::Triangulation; +using CdtTriangulation = CDT::Triangulation; namespace { +// CDT triangulation shared by Python threads, reachable only through lock() +class Triangulation +{ +public: + Triangulation( + CDT::VertexInsertionOrder::Enum vertexInsertionOrder, + CDT::IntersectingConstraintEdges::Enum intersectingEdgesStrategy, + coord_t minDistToConstraintEdge) + : m_cdt( + vertexInsertionOrder, + intersectingEdgesStrategy, + minDistToConstraintEdge) + {} + + // Waits without the GIL, since the thread holding the mutex may need it + std::pair> lock() + { + py::gil_scoped_release release; + return {m_cdt, std::unique_lock(m_mutex)}; + } + +private: + CdtTriangulation m_cdt; + std::mutex m_mutex; +}; + +// Binds a container member as property `name`, `name_count()` and `name_iter()` +template +void def_container( + py::class_& cls, + const std::string& name, + Container CdtTriangulation::*member) +{ + cls.def_property_readonly( + name.c_str(), + [member](Triangulation& t) { + auto [cdt, lock] = t.lock(); + return py::cast( + cdt.*member, + py::return_value_policy::reference_internal, + py::cast(&t)); + }) + .def( + (name + "_count").c_str(), + [member](Triangulation& t) { + auto [cdt, lock] = t.lock(); + return (cdt.*member).size(); + }) + .def( + (name + "_iter").c_str(), + [member](Triangulation& t) { + auto [cdt, lock] = t.lock(); + const Container& container = cdt.*member; + return py::make_iterator(container.begin(), container.end()); + }, + py::keep_alive<0, 1>()); +} + // Binds a vector member as `name_array(*, copy=True)` returning a numpy array template void def_array( py::class_& cls, const std::string& name, - std::vector Triangulation::*member) + std::vector CdtTriangulation::*member) { cls.def( (name + "_array").c_str(), - [member](const Triangulation& t, bool copy) { - const std::vector& items = t.*member; + [member](Triangulation& t, bool copy) { + auto [cdt, lock] = t.lock(); + const std::vector& items = cdt.*member; const auto size = static_cast(items.size()); if (copy) return py::array_t(size, items.data()); @@ -241,170 +301,145 @@ PYBIND11_MODULE(PythonCDT, m) }); py::class_ triangulation(m, "Triangulation"); + triangulation.def( + py::init< + CDT::VertexInsertionOrder::Enum, + CDT::IntersectingConstraintEdges::Enum, + coord_t>(), + py::arg("vertex_insertion_order"), + py::arg("intersecting_edges_strategy"), + py::arg("min_dist_to_constraint_edge")); + def_container(triangulation, "vertices", &CdtTriangulation::vertices); + def_container(triangulation, "triangles", &CdtTriangulation::triangles); + def_container(triangulation, "fixed_edges", &CdtTriangulation::fixedEdges); + def_container( + triangulation, "overlap_count", &CdtTriangulation::overlapCount); + def_container( + triangulation, + "piece_to_originals", + &CdtTriangulation::pieceToOriginals); + def_array(triangulation, "vertices", &CdtTriangulation::vertices); + def_array(triangulation, "triangles", &CdtTriangulation::triangles); triangulation - .def( - py::init< - CDT::VertexInsertionOrder::Enum, - CDT::IntersectingConstraintEdges::Enum, - coord_t>(), - py::arg("vertex_insertion_order"), - py::arg("intersecting_edges_strategy"), - py::arg("min_dist_to_constraint_edge")) - // vertices - .def_readonly("vertices", &Triangulation::vertices) - .def( - "vertices_count", - [](const Triangulation& t) { return t.vertices.size(); }) - .def( - "vertices_iter", - [](const Triangulation& t) -> py::iterator { - return py::make_iterator(t.vertices.begin(), t.vertices.end()); - }, - py::keep_alive<0, 1>()) - // triangles - .def_readonly("triangles", &Triangulation::triangles) - .def( - "triangles_count", - [](const Triangulation& t) { return t.triangles.size(); }) - .def( - "triangles_iter", - [](const Triangulation& t) -> py::iterator { - return py::make_iterator( - t.triangles.begin(), t.triangles.end()); - }, - py::keep_alive<0, 1>()) - // fixed edges - .def_readonly("fixed_edges", &Triangulation::fixedEdges) - .def( - "fixed_edges_count", - [](const Triangulation& t) { return t.fixedEdges.size(); }) - .def( - "fixed_edges_iter", - [](const Triangulation& t) -> py::iterator { - return py::make_iterator( - t.fixedEdges.begin(), t.fixedEdges.end()); - }, - py::keep_alive<0, 1>()) - // overlaps - .def_readonly("overlap_count", &Triangulation::overlapCount) - .def( - "overlap_count_count", - [](const Triangulation& t) { return t.overlapCount.size(); }) - .def( - "overlap_count_iter", - [](const Triangulation& t) -> py::iterator { - return py::make_iterator( - t.overlapCount.begin(), t.overlapCount.end()); - }, - py::keep_alive<0, 1>()) - // piece to originals mapping for edges - .def_readonly("piece_to_originals", &Triangulation::pieceToOriginals) - .def( - "piece_to_originals_count", - [](const Triangulation& t) { return t.pieceToOriginals.size(); }) - .def( - "piece_to_originals_iter", - [](const Triangulation& t) -> py::iterator { - return py::make_iterator( - t.pieceToOriginals.begin(), t.pieceToOriginals.end()); - }, - py::keep_alive<0, 1>()) // methods .def( "insert_vertices", - static_cast&)>( - &Triangulation::insertVertices), - py::arg("vertices"), - py::call_guard()) + [](Triangulation& t, const std::vector& vertices) { + auto [cdt, lock] = t.lock(); + py::gil_scoped_release release; + cdt.insertVertices(vertices); + }, + py::arg("vertices")) .def( "insert_vertices", [](Triangulation& t, py::buffer b) { const py::buffer_info info = b.request(); const auto [ptr, n_vert] = buffer_pairs(info, "double", "vertex"); - { - py::gil_scoped_release release; - t.insertVertices( - ptr, - ptr + n_vert, - [](const BufferPair& v) { return v.v[0]; }, - [](const BufferPair& v) { return v.v[1]; }); - } + auto [cdt, lock] = t.lock(); + py::gil_scoped_release release; + cdt.insertVertices( + ptr, + ptr + n_vert, + [](const BufferPair& v) { return v.v[0]; }, + [](const BufferPair& v) { return v.v[1]; }); }, py::arg("vertex_buffer")) .def( "insert_edges", - static_cast&)>( - &Triangulation::insertEdges), - py::arg("edges"), - py::call_guard()) + [](Triangulation& t, const std::vector& edges) { + auto [cdt, lock] = t.lock(); + py::gil_scoped_release release; + cdt.insertEdges(edges); + }, + py::arg("edges")) .def( "insert_edges", [](Triangulation& t, py::buffer b) { const py::buffer_info info = b.request(); const auto [ptr, n_edges] = buffer_pairs(info, "CDT::VertInd", "edge"); - { - py::gil_scoped_release release; - t.insertEdges( - ptr, - ptr + n_edges, - [](const BufferPair& e) { return e.v[0]; }, - [](const BufferPair& e) { return e.v[1]; }); - } + auto [cdt, lock] = t.lock(); + py::gil_scoped_release release; + cdt.insertEdges( + ptr, + ptr + n_edges, + [](const BufferPair& e) { return e.v[0]; }, + [](const BufferPair& e) { return e.v[1]; }); }, py::arg("edge_buffer")) .def( "conform_to_edges", - static_cast&)>( - &Triangulation::conformToEdges), - py::arg("edges"), - py::call_guard()) + [](Triangulation& t, const std::vector& edges) { + auto [cdt, lock] = t.lock(); + py::gil_scoped_release release; + cdt.conformToEdges(edges); + }, + py::arg("edges")) .def( "conform_to_edges", [](Triangulation& t, py::buffer b) { const py::buffer_info info = b.request(); const auto [ptr, n_edges] = buffer_pairs(info, "CDT::VertInd", "edge"); - { - py::gil_scoped_release release; - t.conformToEdges( - ptr, - ptr + n_edges, - [](const BufferPair& e) { return e.v[0]; }, - [](const BufferPair& e) { return e.v[1]; }); - } + auto [cdt, lock] = t.lock(); + py::gil_scoped_release release; + cdt.conformToEdges( + ptr, + ptr + n_edges, + [](const BufferPair& e) { return e.v[0]; }, + [](const BufferPair& e) { return e.v[1]; }); }, py::arg("edge_buffer")) .def( "erase_super_triangle", - &Triangulation::eraseSuperTriangle, - py::call_guard()) + [](Triangulation& t) { + auto [cdt, lock] = t.lock(); + py::gil_scoped_release release; + cdt.eraseSuperTriangle(); + }) .def( "erase_outer_triangles", - &Triangulation::eraseOuterTriangles, - py::call_guard()) + [](Triangulation& t) { + auto [cdt, lock] = t.lock(); + py::gil_scoped_release release; + cdt.eraseOuterTriangles(); + }) .def( "erase_outer_triangles_and_holes", - &Triangulation::eraseOuterTrianglesAndHoles, - py::call_guard()) - .def("is_finalized", &Triangulation::isFinalized) + [](Triangulation& t) { + auto [cdt, lock] = t.lock(); + py::gil_scoped_release release; + cdt.eraseOuterTrianglesAndHoles(); + }) + .def( + "is_finalized", + [](Triangulation& t) { + auto [cdt, lock] = t.lock(); + return cdt.isFinalized(); + }) .def( "calculate_triangle_depths", - &Triangulation::calculateTriangleDepths, - py::call_guard()) + [](Triangulation& t) { + auto [cdt, lock] = t.lock(); + py::gil_scoped_release release; + return cdt.calculateTriangleDepths(); + }) .def( "remove_triangles", - static_cast( - &Triangulation::removeTriangles), - py::arg("triangle_indices"), - py::call_guard()); - def_array(triangulation, "vertices", &Triangulation::vertices); - def_array(triangulation, "triangles", &Triangulation::triangles); + [](Triangulation& t, const CDT::TriIndUSet& triangle_indices) { + auto [cdt, lock] = t.lock(); + py::gil_scoped_release release; + cdt.removeTriangles(triangle_indices); + }, + py::arg("triangle_indices")); m.def( "verify_topology", - &CDT::verifyTopology, - py::arg("triangulation"), - py::call_guard()); + [](Triangulation& t) { + auto [cdt, lock] = t.lock(); + py::gil_scoped_release release; + return CDT::verifyTopology(cdt); + }, + py::arg("triangulation")); } diff --git a/cdt_bindings_test.py b/cdt_bindings_test.py index 9f335d3..0a0954a 100644 --- a/cdt_bindings_test.py +++ b/cdt_bindings_test.py @@ -10,7 +10,7 @@ import pytest import tempfile import hashlib -import threading +from concurrent.futures import ThreadPoolExecutor import PythonCDT as cdt @@ -74,6 +74,12 @@ def test_Triangulation() -> None: assert len(t.triangles) == 2, "Wrong triangle count in triangulation" assert len(t.fixed_edges) == 1, "Wrong fixed edge count in triangulation" + # properties give references into the triangulation + vertices = t.vertices + vertices[3].x = 42 + assert t.vertices[3] == cdt.V2d(42, -0.5), "Vertex property must reference triangulation's vertices" + vertices[3].x = 0 + # test retrieving triangulation data using iterators assert t.vertices_count() == len(t.vertices), "Wrong vertex count" assert t.triangles_count() == len(t.triangles), "Wrong triangle count" @@ -213,30 +219,19 @@ def test_arrays(copy) -> None: assert np.array_equal(arr, before), "Array must stay valid after the triangulation is deleted" -def test_insert_releases_the_gil() -> None: - """Concurrent builds on threads must each equal the single-threaded build.""" - vv, ee = read_input_file("CDT/visualizer/data/Constrained Sweden.txt") - verts = np.array([[v.x, v.y] for v in vv], dtype=np.float64) - edges = np.array([[e.v1, e.v2] for e in ee], dtype=np.uintc) - - def build() -> str: - t = cdt.Triangulation(cdt.VertexInsertionOrder.AS_PROVIDED, cdt.IntersectingConstraintEdges.TRY_RESOLVE, 0.0) - t.insert_vertices(verts) - t.insert_edges(edges) - t.erase_super_triangle() - assert cdt.verify_topology(t) - return triangulation_md5_checksum(t) - - expected = build() - results = [None] * 4 - def worker(i): - results[i] = build() - threads = [threading.Thread(target=worker, args=(i,)) for i in range(4)] - for th in threads: - th.start() - for th in threads: - th.join() - assert results == [expected] * 4, "Threaded builds must match the single-threaded build" +def test_shared_triangulation() -> None: + """Threads can insert into and read from one triangulation concurrently""" + vertices = np.random.default_rng(0).random((160_000, 2)) + t = cdt.Triangulation(cdt.VertexInsertionOrder.AUTO, cdt.IntersectingConstraintEdges.NOT_ALLOWED, 0.0) + + def insert_and_read(batch) -> None: + t.insert_vertices(batch) + assert t.triangles_array()["vertices"].max() < t.vertices_count() + + with ThreadPoolExecutor(max_workers=8) as pool: + list(pool.map(insert_and_read, np.array_split(vertices, 8))) + assert t.vertices_count() == len(vertices) + 3, "All vertices plus the super-triangle must be inserted" + assert cdt.verify_topology(t) def test_insert_buffers_must_be_contiguous_pairs() -> None: From 33316696d85177ae8a73240ecb2e34451144452d Mon Sep 17 00:00:00 2001 From: artem-ogre Date: Mon, 14 Sep 2026 23:59:18 +0200 Subject: [PATCH 7/8] Add GIL release benchmark --- benchmark_gil_release.py | 47 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) create mode 100644 benchmark_gil_release.py diff --git a/benchmark_gil_release.py b/benchmark_gil_release.py new file mode 100644 index 0000000..8aa9398 --- /dev/null +++ b/benchmark_gil_release.py @@ -0,0 +1,47 @@ +"""Shows the effect of releasing the GIL during triangulation""" + +import os +import threading +import time +from concurrent.futures import ThreadPoolExecutor + +import numpy as np + +import PythonCDT as cdt + +N_THREADS = min(8, os.cpu_count()) +N_VERTICES = 200_000 + + +def triangulate(vertices) -> None: + t = cdt.Triangulation(cdt.VertexInsertionOrder.AUTO, cdt.IntersectingConstraintEdges.NOT_ALLOWED, 0.0) + t.insert_vertices(vertices) + + +def seconds(f) -> float: + start = time.perf_counter() + f() + return time.perf_counter() - start + + +rng = np.random.default_rng(0) +batches = [rng.random((N_VERTICES, 2)) for _ in range(N_THREADS)] + +# separate triangulations built in parallel threads +sequential = seconds(lambda: [triangulate(b) for b in batches]) +with ThreadPoolExecutor(N_THREADS) as pool: + threaded = seconds(lambda: list(pool.map(triangulate, batches))) +print(f"{N_THREADS} triangulations of {N_VERTICES} vertices each") +print(f" sequential {sequential:8.3f} s") +print(f" {N_THREADS} threads {threaded:8.3f} s ({sequential / threaded:.1f}x speedup)") + +# a Python thread keeps running while another thread triangulates +worker = threading.Thread(target=triangulate, args=(batches[0],)) +ticks = [time.perf_counter()] +worker.start() +while worker.is_alive(): + time.sleep(0.001) + ticks.append(time.perf_counter()) +longest_stall = max(np.diff(ticks)) +print(f"main thread ticking every 1 ms during {ticks[-1] - ticks[0]:.3f} s triangulation") +print(f" ticks {len(ticks) - 1:>6}, longest stall {1e3 * longest_stall:8.1f} ms") From d42f7b2d488e2294ec4d53cafad2f62363c06275 Mon Sep 17 00:00:00 2001 From: artem-ogre Date: Tue, 15 Sep 2026 22:07:30 +0200 Subject: [PATCH 8/8] Make locking more elegant --- cdt_bindings.cpp | 129 +++++++++++++++++++++++------------------------ 1 file changed, 62 insertions(+), 67 deletions(-) diff --git a/cdt_bindings.cpp b/cdt_bindings.cpp index 53933ef..02e598c 100644 --- a/cdt_bindings.cpp +++ b/cdt_bindings.cpp @@ -29,30 +29,33 @@ using CdtTriangulation = CDT::Triangulation; namespace { -// CDT triangulation shared by Python threads, reachable only through lock() -class Triangulation +struct Triangulation : CdtTriangulation { -public: - Triangulation( - CDT::VertexInsertionOrder::Enum vertexInsertionOrder, - CDT::IntersectingConstraintEdges::Enum intersectingEdgesStrategy, - coord_t minDistToConstraintEdge) - : m_cdt( - vertexInsertionOrder, - intersectingEdgesStrategy, - minDistToConstraintEdge) + using CdtTriangulation::CdtTriangulation; + std::mutex mutex; +}; + +struct LockWithoutGil +{ + explicit LockWithoutGil(Triangulation& t) + : lock(t.mutex) {} + py::gil_scoped_release release; // constructed first, destroyed last + std::lock_guard lock; +}; - // Waits without the GIL, since the thread holding the mutex may need it - std::pair> lock() +struct LockWithGil +{ + explicit LockWithGil(Triangulation& t) + : lock(t.mutex, std::try_to_lock) { - py::gil_scoped_release release; - return {m_cdt, std::unique_lock(m_mutex)}; + if(!lock.owns_lock()) + { + py::gil_scoped_release release; + lock.lock(); + } } - -private: - CdtTriangulation m_cdt; - std::mutex m_mutex; + std::unique_lock lock; }; // Binds a container member as property `name`, `name_count()` and `name_iter()` @@ -65,26 +68,30 @@ void def_container( cls.def_property_readonly( name.c_str(), [member](Triangulation& t) { - auto [cdt, lock] = t.lock(); + LockWithGil lock(t); return py::cast( - cdt.*member, + t.*member, py::return_value_policy::reference_internal, py::cast(&t)); - }) + }, + "Items reference the triangulation's memory; they are invalidated " + "by any call that modifies the triangulation.") .def( (name + "_count").c_str(), [member](Triangulation& t) { - auto [cdt, lock] = t.lock(); - return (cdt.*member).size(); + LockWithGil lock(t); + return (t.*member).size(); }) .def( (name + "_iter").c_str(), [member](Triangulation& t) { - auto [cdt, lock] = t.lock(); - const Container& container = cdt.*member; + LockWithGil lock(t); + const Container& container = t.*member; return py::make_iterator(container.begin(), container.end()); }, - py::keep_alive<0, 1>()); + py::keep_alive<0, 1>(), + "Iterates the triangulation's memory; the iterator is invalidated " + "by any call that modifies the triangulation."); } // Binds a vector member as `name_array(*, copy=True)` returning a numpy array @@ -97,8 +104,8 @@ void def_array( cls.def( (name + "_array").c_str(), [member](Triangulation& t, bool copy) { - auto [cdt, lock] = t.lock(); - const std::vector& items = cdt.*member; + LockWithGil lock(t); + const std::vector& items = t.*member; const auto size = static_cast(items.size()); if (copy) return py::array_t(size, items.data()); @@ -325,9 +332,8 @@ PYBIND11_MODULE(PythonCDT, m) .def( "insert_vertices", [](Triangulation& t, const std::vector& vertices) { - auto [cdt, lock] = t.lock(); - py::gil_scoped_release release; - cdt.insertVertices(vertices); + LockWithoutGil lock(t); + t.insertVertices(vertices); }, py::arg("vertices")) .def( @@ -336,9 +342,8 @@ PYBIND11_MODULE(PythonCDT, m) const py::buffer_info info = b.request(); const auto [ptr, n_vert] = buffer_pairs(info, "double", "vertex"); - auto [cdt, lock] = t.lock(); - py::gil_scoped_release release; - cdt.insertVertices( + LockWithoutGil lock(t); + t.insertVertices( ptr, ptr + n_vert, [](const BufferPair& v) { return v.v[0]; }, @@ -348,9 +353,8 @@ PYBIND11_MODULE(PythonCDT, m) .def( "insert_edges", [](Triangulation& t, const std::vector& edges) { - auto [cdt, lock] = t.lock(); - py::gil_scoped_release release; - cdt.insertEdges(edges); + LockWithoutGil lock(t); + t.insertEdges(edges); }, py::arg("edges")) .def( @@ -359,9 +363,8 @@ PYBIND11_MODULE(PythonCDT, m) const py::buffer_info info = b.request(); const auto [ptr, n_edges] = buffer_pairs(info, "CDT::VertInd", "edge"); - auto [cdt, lock] = t.lock(); - py::gil_scoped_release release; - cdt.insertEdges( + LockWithoutGil lock(t); + t.insertEdges( ptr, ptr + n_edges, [](const BufferPair& e) { return e.v[0]; }, @@ -371,9 +374,8 @@ PYBIND11_MODULE(PythonCDT, m) .def( "conform_to_edges", [](Triangulation& t, const std::vector& edges) { - auto [cdt, lock] = t.lock(); - py::gil_scoped_release release; - cdt.conformToEdges(edges); + LockWithoutGil lock(t); + t.conformToEdges(edges); }, py::arg("edges")) .def( @@ -382,9 +384,8 @@ PYBIND11_MODULE(PythonCDT, m) const py::buffer_info info = b.request(); const auto [ptr, n_edges] = buffer_pairs(info, "CDT::VertInd", "edge"); - auto [cdt, lock] = t.lock(); - py::gil_scoped_release release; - cdt.conformToEdges( + LockWithoutGil lock(t); + t.conformToEdges( ptr, ptr + n_edges, [](const BufferPair& e) { return e.v[0]; }, @@ -394,52 +395,46 @@ PYBIND11_MODULE(PythonCDT, m) .def( "erase_super_triangle", [](Triangulation& t) { - auto [cdt, lock] = t.lock(); - py::gil_scoped_release release; - cdt.eraseSuperTriangle(); + LockWithoutGil lock(t); + t.eraseSuperTriangle(); }) .def( "erase_outer_triangles", [](Triangulation& t) { - auto [cdt, lock] = t.lock(); - py::gil_scoped_release release; - cdt.eraseOuterTriangles(); + LockWithoutGil lock(t); + t.eraseOuterTriangles(); }) .def( "erase_outer_triangles_and_holes", [](Triangulation& t) { - auto [cdt, lock] = t.lock(); - py::gil_scoped_release release; - cdt.eraseOuterTrianglesAndHoles(); + LockWithoutGil lock(t); + t.eraseOuterTrianglesAndHoles(); }) .def( "is_finalized", [](Triangulation& t) { - auto [cdt, lock] = t.lock(); - return cdt.isFinalized(); + LockWithGil lock(t); + return t.isFinalized(); }) .def( "calculate_triangle_depths", [](Triangulation& t) { - auto [cdt, lock] = t.lock(); - py::gil_scoped_release release; - return cdt.calculateTriangleDepths(); + LockWithoutGil lock(t); + return t.calculateTriangleDepths(); }) .def( "remove_triangles", [](Triangulation& t, const CDT::TriIndUSet& triangle_indices) { - auto [cdt, lock] = t.lock(); - py::gil_scoped_release release; - cdt.removeTriangles(triangle_indices); + LockWithoutGil lock(t); + t.removeTriangles(triangle_indices); }, py::arg("triangle_indices")); m.def( "verify_topology", [](Triangulation& t) { - auto [cdt, lock] = t.lock(); - py::gil_scoped_release release; - return CDT::verifyTopology(cdt); + LockWithoutGil lock(t); + return CDT::verifyTopology(t); }, py::arg("triangulation")); }