Update bindings to CDT 2.0.0

CMakeLists.txt

@@ -4,7 +4,7 @@ set(CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake" ${CMAKE_MODULE_PATH})
 
 project(
     PythonCDT
-    VERSION 1.4.5
+    VERSION 2.0.0
     DESCRIPTION "Python bindings for CDT: constrained Delaunay triangulation"
     LANGUAGES CXX
 )

README.md

@@ -16,24 +16,65 @@ pip3 install .
 pytest ./cdt_bindings_test.py
 ```
 
-## Usage
+## Usage examples
+
+### Constrained triangulation
 
 ```python
 import numpy as np
 import PythonCDT as cdt
 
+vertices = np.array([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0], [0.4, 0.4]])
+edges = np.array([[0, 1], [1, 2], [2, 3], [3, 0]], dtype=np.uintc)
+
+t = cdt.Triangulation(cdt.VertexInsertionOrder.AUTO, cdt.IntersectingConstraintEdges.TRY_RESOLVE, 0.0)
+t.insert_vertices(vertices)
+t.insert_edges(edges)
+t.erase_outer_triangles_and_holes()
+
+vv = t.vertices_array()   # numpy array with fields 'x' and 'y'
+tt = t.triangles_array()  # numpy array with fields 'vertices' and 'neighbors'
+tt["vertices"]            # (T, 3) vertex indices into vv
+```
+
+#### Notes
+
+- `vertices_array()` and `triangles_array()` return copies. With `copy=False` they return read-only views of the
+  triangulation's memory instead; a view is invalidated by any call that modifies the triangulation.
+
+### Conforming triangulation
+
+```python
 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()
+t.insert_vertices(vertices)
+t.conform_to_edges(edges)
+t.erase_outer_triangles_and_holes()
+```
+
+#### Notes
+
+- `conform_to_edges()` splits the constraint edges as needed, instead of keeping them as they are.
+
+### Refined triangulation
 
-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
+```python
+t = cdt.Triangulation(cdt.VertexInsertionOrder.AUTO, cdt.IntersectingConstraintEdges.TRY_RESOLVE, 0.0)
+t.insert_vertices(vertices)
+t.insert_edges(edges)
+
+to_erase = t.collect_outer_triangles_and_holes()
+unrefined = t.refine_triangles(1000, cdt.RefinementCriterion.SMALLEST_ANGLE, cdt.deg_to_rad(20.0), to_erase)
+t.finalize_triangulation(to_erase)
 ```
 
-`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.
+#### Notes
+
+- `refine_triangles()` improves the shape of the triangles by inserting new points (Steiner points). It must be called
+  before the triangulation is finalized: collect the triangles to remove first so that they are not refined, then
+  remove them with `finalize_triangulation()` (the set is updated in place).
+- Some places can not be refined: e.g., a sharp angle between two constraint edges comes from the input and can not be
+  made any larger. `refine_triangles()` returns the counts of such refinements; `find_unrefined_triangles()` and
+  `find_encroached_fixed_edges()` locate them in the resulting triangulation.
 
 ### Threads
 
@@ -46,3 +87,4 @@ Iterators (`*_iter()`) and `copy=False` views are not protected: don't use them
 
 ## Contributors
 - [SioulisChris](https://github.com/SioulisChris): fixing the tests on Windows
+- [sccolbert](https://github.com/sccolbert): reading the triangulation back as numpy arrays, releasing the GIL

cdt_bindings.cpp

@@ -202,18 +202,41 @@ PYBIND11_MODULE(PythonCDT, m)
         .value("TRY_RESOLVE", CDT::IntersectingConstraintEdges::TryResolve)
         .value("DONT_CHECK", CDT::IntersectingConstraintEdges::DontCheck);
 
+    py::enum_<CDT::RefinementCriterion::Enum>(m, "RefinementCriterion")
+        .value("SMALLEST_ANGLE", CDT::RefinementCriterion::SmallestAngle)
+        .value("LARGEST_AREA", CDT::RefinementCriterion::LargestArea);
+
+    py::class_<CDT::Unrefined>(m, "Unrefined")
+        .def_readonly(
+            "short_edge_triangles", &CDT::Unrefined::shortEdgeTriangles)
+        .def_readonly(
+            "circumcenter_outside", &CDT::Unrefined::circumcenterOutside)
+        .def_readonly(
+            "circumcenter_on_vertex", &CDT::Unrefined::circumcenterOnVertex)
+        .def_readonly("sharp_fixed_corner", &CDT::Unrefined::sharpFixedCorner)
+        .def_readonly("short_edges", &CDT::Unrefined::shortEdges)
+        .def_readonly(
+            "split_vertex_invalid", &CDT::Unrefined::splitVertexInvalid)
+        .def("__repr__", [](const CDT::Unrefined& u) {
+            std::ostringstream oss;
+            oss << "Unrefined(short_edge_triangles=" << u.shortEdgeTriangles
+                << ", circumcenter_outside=" << u.circumcenterOutside
+                << ", circumcenter_on_vertex=" << u.circumcenterOnVertex
+                << ", sharp_fixed_corner=" << u.sharpFixedCorner
+                << ", short_edges=" << u.shortEdges
+                << ", split_vertex_invalid=" << u.splitVertexInvalid << ")";
+            return oss.str();
+        });
+
     py::class_<V2d>(m, "V2d", py::buffer_protocol())
         .def(py::init<coord_t, coord_t>(), py::arg("x"), py::arg("y"))
         .def(py::init([](py::buffer b) {
-            // Request a buffer descriptor from Python
             py::buffer_info info = b.request();
-            // Some sanity checks ...
             if (info.format != py::format_descriptor<coord_t>::format())
                 throw std::runtime_error(
                     "Incompatible format: expected a double array!");
             if (info.ndim != 1)
                 throw std::runtime_error("Incompatible buffer dimension!");
-            // create from buffer
             const coord_t* const ptr = static_cast<coord_t*>(info.ptr);
             return V2d{ptr[0], ptr[1]};
         }))
@@ -300,6 +323,7 @@ PYBIND11_MODULE(PythonCDT, m)
         .def_property_readonly("v2", &CDT::Edge::v2)
         .def(py::self == py::self)
         .def(py::self != py::self)
+        .def(py::self < py::self)
         .def(py::hash(py::self))
         .def("__repr__", [](const CDT::Edge& e) {
             std::ostringstream oss;
@@ -428,7 +452,83 @@ PYBIND11_MODULE(PythonCDT, m)
                 LockWithoutGil lock(t);
                 t.removeTriangles(triangle_indices);
             },
-            py::arg("triangle_indices"));
+            py::arg("triangle_indices"))
+        .def(
+            "refine_triangles",
+            [](Triangulation& t,
+               CDT::VertInd max_vertices_to_insert,
+               CDT::RefinementCriterion::Enum refinement_criterion,
+               coord_t refinement_threshold,
+               py::object to_erase,
+               coord_t min_edge_length) {
+                const bool has_to_erase = !to_erase.is_none();
+                if (has_to_erase && !py::isinstance<py::set>(to_erase))
+                    throw py::type_error("to_erase must be a set or None!");
+                CDT::TriIndUSet triangles_to_erase;
+                if (has_to_erase)
+                    triangles_to_erase = to_erase.cast<CDT::TriIndUSet>();
+                CDT::Unrefined unrefined;
+                {
+                    LockWithoutGil lock(t);
+                    unrefined = t.refineTriangles(
+                        max_vertices_to_insert,
+                        refinement_criterion,
+                        refinement_threshold,
+                        has_to_erase ? &triangles_to_erase : NULL,
+                        min_edge_length);
+                }
+                if (has_to_erase)
+                {
+                    py::set marked = py::reinterpret_borrow<py::set>(to_erase);
+                    marked.clear();
+                    for (const CDT::TriInd it : triangles_to_erase)
+                        marked.add(py::cast(it));
+                }
+                return unrefined;
+            },
+            py::arg("max_vertices_to_insert"),
+            py::arg("refinement_criterion") =
+                CDT::RefinementCriterion::SmallestAngle,
+            py::arg("refinement_threshold") = CDT::degToRad(coord_t(20)),
+            py::arg("to_erase") = py::none(),
+            py::arg("min_edge_length") = coord_t(1e-6),
+            "Triangles from `to_erase` are not refined; the triangles that "
+            "replace them are added to it in place.")
+        .def("find_encroached_fixed_edges", [](Triangulation& t) {
+            LockWithoutGil lock(t);
+            return t.findEncroachedFixedEdges();
+        })
+        .def(
+            "find_unrefined_triangles",
+            [](Triangulation& t,
+               CDT::RefinementCriterion::Enum refinement_criterion,
+               coord_t refinement_threshold) {
+                LockWithoutGil lock(t);
+                return t.findUnrefinedTriangles(
+                    refinement_criterion, refinement_threshold);
+            },
+            py::arg("refinement_criterion") =
+                CDT::RefinementCriterion::SmallestAngle,
+            py::arg("refinement_threshold") = CDT::degToRad(coord_t(20)))
+        .def("collect_super_triangle", [](Triangulation& t) {
+            LockWithoutGil lock(t);
+            return t.collectSuperTriangle();
+        })
+        .def("collect_outer_triangles", [](Triangulation& t) {
+            LockWithoutGil lock(t);
+            return t.collectOuterTriangles();
+        })
+        .def("collect_outer_triangles_and_holes", [](Triangulation& t) {
+            LockWithoutGil lock(t);
+            return t.collectOuterTrianglesAndHoles();
+        })
+        .def(
+            "finalize_triangulation",
+            [](Triangulation& t, const CDT::TriIndUSet& removed_triangles) {
+                LockWithoutGil lock(t);
+                t.finalizeTriangulation(removed_triangles);
+            },
+            py::arg("removed_triangles"));
 
     m.def(
         "verify_topology",
@@ -437,4 +537,35 @@ PYBIND11_MODULE(PythonCDT, m)
             return CDT::verifyTopology(t);
         },
         py::arg("triangulation"));
+
+    m.def(
+        "verify_winding",
+        [](Triangulation& t) {
+            LockWithoutGil lock(t);
+            return CDT::verifyWinding(t);
+        },
+        py::arg("triangulation"));
+
+    m.def("deg_to_rad", &CDT::degToRad<coord_t>, py::arg("degrees"));
+
+    m.def(
+        "area",
+        &CDT::area<coord_t>,
+        py::arg("a"),
+        py::arg("b"),
+        py::arg("c"));
+
+    m.def(
+        "circumcenter",
+        &CDT::circumcenter<coord_t>,
+        py::arg("a"),
+        py::arg("b"),
+        py::arg("c"));
+
+    m.def(
+        "smallest_angle",
+        &CDT::smallestAngle<coord_t>,
+        py::arg("a"),
+        py::arg("b"),
+        py::arg("c"));
 }

cdt_bindings_test.py

@@ -49,6 +49,9 @@ def test_Edge() -> None:
 
     assert cdt.Edge(1, 2).__repr__() == "Edge(1, 2)", "Wrong __repr__ output for Edge"
 
+    ee = [cdt.Edge(2, 3), cdt.Edge(0, 5), cdt.Edge(0, 1)]
+    assert sorted(ee) == [cdt.Edge(0, 1), cdt.Edge(0, 5), cdt.Edge(2, 3)], "Edges are ordered wrong"
+
 
 def test_Triangulation() -> None:
     """Test Triangulation class"""
@@ -115,6 +118,24 @@ def test_verify_topology() -> None:
     assert cdt.verify_topology(t), "Verifying topology produced wrong result"
 
 
+def test_verify_winding() -> None:
+    t = cdt.Triangulation(cdt.VertexInsertionOrder.AS_PROVIDED, cdt.IntersectingConstraintEdges.TRY_RESOLVE, 0.0)
+    t.insert_vertices([cdt.V2d(-1, 0), cdt.V2d(0, 0.5), cdt.V2d(1, 0), cdt.V2d(0, -0.5)])
+    assert cdt.verify_winding(t), "Verifying winding produced wrong result"
+
+    t.triangles[0].vertices = t.triangles[0].vertices[::-1]
+    assert not cdt.verify_winding(t), "Inverted triangle must not pass winding verification"
+
+
+def test_triangle_geometry() -> None:
+    assert cdt.deg_to_rad(180.0) == np.pi, "Wrong degrees to radians conversion"
+
+    a, b, c = cdt.V2d(0, 0), cdt.V2d(1, 0), cdt.V2d(0, 1)
+    assert cdt.area(a, b, c) == 0.5, "Wrong triangle area"
+    assert cdt.smallest_angle(a, b, c) == pytest.approx(np.pi / 4), "Wrong smallest triangle angle"
+    assert cdt.circumcenter(a, b, c) == cdt.V2d(0.5, 0.5), "Wrong triangle circumcenter"
+
+
 def save_triangulation_as_off(t: cdt.Triangulation, off_file) -> None:
     with open(off_file, "w") as f:
         f.write(f"OFF\n")
@@ -164,6 +185,97 @@ def test_conform_to_edges() -> None:
     assert triangulation_md5_checksum(t) == 'b64cae39c91a55dd4e23a146eb7df0d3', "Wrong OFF file contents"
 
 
+def triangle_smallest_angles(t: cdt.Triangulation):
+    return [cdt.smallest_angle(*(t.vertices[int(i)] for i in tri.vertices)) for tri in t.triangles_iter()]
+
+
+def triangle_areas(t: cdt.Triangulation):
+    return [cdt.area(*(t.vertices[int(i)] for i in tri.vertices)) for tri in t.triangles_iter()]
+
+
+def triangulation_with_bad_triangles() -> cdt.Triangulation:
+    t = cdt.Triangulation(cdt.VertexInsertionOrder.AS_PROVIDED, cdt.IntersectingConstraintEdges.NOT_ALLOWED, 0.0)
+    t.insert_vertices(np.array([[0.0, 0.0], [10.0, 0.0], [10.0, 1.0], [0.0, 1.0]]))
+    return t
+
+
+def triangulation_with_sharp_input_corner() -> cdt.Triangulation:
+    t = cdt.Triangulation(cdt.VertexInsertionOrder.AS_PROVIDED, cdt.IntersectingConstraintEdges.NOT_ALLOWED, 0.0)
+    t.insert_vertices(np.array([[0.0, 0.0], [10.0, 0.0], [10.0, 1.0]]))
+    t.insert_edges(np.array([[0, 1], [1, 2], [2, 0]], dtype=np.uintc))
+    return t
+
+
+@pytest.mark.parametrize("criterion, threshold, is_criterion_fulfilled", [
+    (cdt.RefinementCriterion.SMALLEST_ANGLE, cdt.deg_to_rad(20.0),
+     lambda t, threshold: min(triangle_smallest_angles(t)) >= threshold),
+    (cdt.RefinementCriterion.LARGEST_AREA, 2.0,
+     lambda t, threshold: max(triangle_areas(t)) <= threshold),
+])
+def test_refine_triangles_fulfills_criterion(criterion, threshold, is_criterion_fulfilled) -> None:
+    vv, ee = read_input_file("CDT/visualizer/data/Capital A.txt")
+    t = cdt.Triangulation(cdt.VertexInsertionOrder.AS_PROVIDED, cdt.IntersectingConstraintEdges.NOT_ALLOWED, 0.0)
+    t.insert_vertices(vv)
+    t.insert_edges(ee)
+
+    to_erase = t.collect_outer_triangles_and_holes()
+    assert to_erase == t.collect_outer_triangles_and_holes(), "Collecting triangles must not change the triangulation"
+    unrefined = t.refine_triangles(1000, criterion, threshold, to_erase)
+    gave_up = [unrefined.short_edge_triangles, unrefined.circumcenter_outside, unrefined.circumcenter_on_vertex,
+               unrefined.sharp_fixed_corner, unrefined.short_edges, unrefined.split_vertex_invalid]
+    assert not any(gave_up), f"Refinement gave up: {unrefined}"
+    assert t.find_encroached_fixed_edges() == [], "Triangulation has encroached fixed edges"
+
+    t.finalize_triangulation(to_erase)
+    assert t.find_unrefined_triangles(criterion, threshold) == [], "Kept triangles must all be refined"
+    assert cdt.verify_topology(t) and cdt.verify_winding(t), "Refinement broke the triangulation"
+    assert is_criterion_fulfilled(t, threshold), "Kept triangles must fulfill the refinement criterion"
+
+
+def test_refine_triangles_with_zero_threshold_inserts_nothing() -> None:
+    t = triangulation_with_bad_triangles()
+    n_vertices_before = t.vertices_count()
+    t.refine_triangles(1000, cdt.RefinementCriterion.SMALLEST_ANGLE, 0.0)
+    assert t.vertices_count() == n_vertices_before, "Any triangle fulfills a zero threshold"
+
+
+def test_refine_triangles_stops_at_vertex_budget() -> None:
+    t = triangulation_with_bad_triangles()
+    n_vertices_before = t.vertices_count()
+    t.refine_triangles(5, cdt.RefinementCriterion.SMALLEST_ANGLE, cdt.deg_to_rad(30.0))
+    assert t.vertices_count() == n_vertices_before + 5, "Refinement must insert exactly the budgeted vertices"
+
+
+def test_refine_triangles_reports_sharp_input_corner() -> None:
+    t = triangulation_with_sharp_input_corner()
+    threshold = cdt.deg_to_rad(20.0)
+    to_erase = t.collect_outer_triangles()
+    unrefined = t.refine_triangles(1000, cdt.RefinementCriterion.SMALLEST_ANGLE, threshold, to_erase)
+    assert unrefined.sharp_fixed_corner > 0, f"Sharp input corner was not reported: {unrefined}"
+    assert len(t.find_unrefined_triangles(cdt.RefinementCriterion.SMALLEST_ANGLE, threshold)) == 1, \
+        "Triangle left sharp by the input must be found"
+
+    t.finalize_triangulation(to_erase)
+    assert t.triangles_count() == 1, "Only the sharp input triangle must be left"
+
+
+def test_refine_triangles_rejects_to_erase_that_is_not_a_set() -> None:
+    t = triangulation_with_bad_triangles()
+    with pytest.raises(TypeError):
+        t.refine_triangles(10, to_erase=[0, 1])
+
+
+def test_refining_and_collecting_fail_on_finalized_triangulation() -> None:
+    t = triangulation_with_bad_triangles()
+    t.erase_super_triangle()
+    assert t.is_finalized(), "Triangulation must be finalized"
+    for call in [lambda: t.refine_triangles(10), t.find_encroached_fixed_edges, t.collect_super_triangle,
+                 t.collect_outer_triangles, t.collect_outer_triangles_and_holes,
+                 lambda: t.finalize_triangulation(set())]:
+        with pytest.raises(RuntimeError):
+            call()
+
+
 @pytest.mark.parametrize("vv", [[cdt.V2d(-1, 0), cdt.V2d(0, 0.5), cdt.V2d(1, 0), cdt.V2d(0, -0.5)],
                                 np.array([[-1, 0], [0, 0.5], [1, 0], [0, -0.5]], dtype=np.float64),
                                 np.array([-1, 0, 0, 0.5, 1, 0, 0, -0.5], dtype=np.float64)])