Add a Delaunay-property test

pyproject.toml

@@ -42,9 +42,10 @@ Source = "https://github.com/ybubnov/torch_delaunay"
 
 [project.optional-dependencies]
 test = [
+  "geopandas >= 1.0.0",
+  "numpy >= 1.4.0, <2.0.0",
   "pytest >= 8.0.0",
   "shapely >= 2.0.0",
-  "numpy >= 1.4.0, <2.0.0"
 ]
 
 

src/triangle.cc

@@ -66,99 +66,6 @@ dist(const torch::Tensor& p, const torch::Tensor& q)
 }
 
 
-void
-_lawson_flip_out(
-    const torch::Tensor& triangles,
-    const torch::Tensor& halfedges,
-    const torch::Tensor& points,
-    const torch::Tensor& pr,
-    const torch::Tensor& pl,
-    const torch::Tensor& p0,
-    torch::Tensor& out
-)
-{
-    auto index_twin = halfedges.index_select(0, pl).index_select(1, pr);
-    std::cout << "OPPOSITE" << std::endl << index_twin.to_dense() << std::endl;
-
-    // TODO: implement diag for sparse coordinate tensor.
-    auto twins = index_twin.to_dense().diag();
-    twins = twins - 1;
-    std::cout << "TWINS" << std::endl << twins << std::endl;
-
-    twins = torch::where(twins.ne(-1), twins, p0);
-    std::cout << "TWINS (defaulted)" << std::endl << twins << std::endl;
-
-    // torch::Tensor p1 = triangles.index({twins - 1, dim});
-    const auto p1 = twins;
-    std::cout << "P1: " << p1 << std::endl;
-
-    const auto incircle = incircle2d(
-        points.index({p0}), points.index({pr}), points.index({pl}), points.index({p1})
-    );
-    // std::cout << "IN CIRCLE" << std::endl << incircle << std::endl;
-    std::cout << "PP: " << torch::column_stack({p0, pr, pl, p1, incircle}) << std::endl;
-
-    for (int64_t i = 0; i < out.sizes()[0]; i++) {
-        if ((incircle[i].item<int64_t>() > 0)) {
-            std::cout << "flip " << i << std::endl;
-            out.index_put_({i, 0}, p0[i].item<int64_t>());
-            out.index_put_({i, 1}, pr[i].item<int64_t>());
-            out.index_put_({i, 2}, p1[i].item<int64_t>());
-        }
-    }
-}
-
-
-torch::Tensor
-lawson_flip(const torch::Tensor& triangles, const torch::Tensor& points)
-{
-    torch::Tensor out = triangles.clone();
-
-    const auto n = triangles.size(0);
-
-    torch::Tensor he_indices = torch::full({2, n * 3}, -1, triangles.options());
-    torch::Tensor he_values = torch::full({n * 3}, -1, triangles.options());
-
-    for (const auto i : c10::irange(n)) {
-        auto v0 = triangles.index({i, 0}).item<int64_t>();
-        auto v1 = triangles.index({i, 1}).item<int64_t>();
-        auto v2 = triangles.index({i, 2}).item<int64_t>();
-
-        std::cout << i;
-        he_indices.index_put_({0, i * 3 + 0}, v0);
-        he_indices.index_put_({1, i * 3 + 0}, v1);
-
-        he_indices.index_put_({0, i * 3 + 1}, v1);
-        he_indices.index_put_({1, i * 3 + 1}, v2);
-
-        he_indices.index_put_({0, i * 3 + 2}, v2);
-        he_indices.index_put_({1, i * 3 + 2}, v0);
-        std::cout << " ind... ";
-
-        he_values.index_put_({i * 3 + 0}, v2 + 1);
-        he_values.index_put_({i * 3 + 1}, v0 + 1);
-        he_values.index_put_({i * 3 + 2}, v1 + 1);
-        std::cout << "vals..." << std::endl;
-    }
-
-    const auto halfedges = torch::sparse_coo_tensor(he_indices, he_values);
-    std::cout << "REV INDEX" << std::endl << halfedges.to_dense() << std::endl;
-
-    torch::Tensor pr = triangles.index({Slice(), 0}).contiguous();
-    torch::Tensor pl = triangles.index({Slice(), 1}).contiguous();
-    torch::Tensor p0 = triangles.index({Slice(), 2}).contiguous();
-
-    _lawson_flip_out(triangles, halfedges, points, pr, pl, p0, out);
-    std::cout << "OUT" << std::endl << out << std::endl;
-    _lawson_flip_out(triangles, halfedges, points, p0, pr, pl, out);
-    std::cout << "OUT" << std::endl << out << std::endl;
-    _lawson_flip_out(triangles, halfedges, points, pl, p0, pr, out);
-    std::cout << "OUT" << std::endl << out << std::endl;
-
-    return out;
-}
-
-
 struct _SHull {
     std::vector<int64_t> hash;
     std::vector<int64_t> triangles;
@@ -297,8 +204,7 @@ struct _SHull {
             auto pl = points[triangles[al]].unsqueeze(0);
             auto p1 = points[triangles[bl]].unsqueeze(0);
 
-            // TODO: is it >0 or <0?
-            if (incircle2d(p0, pr, pl, p1).gt(0).all().item<bool>()) {
+            if (incircle2d(p0, pr, pl, p1).lt(0).all().item<bool>()) {
                 triangles[a] = triangles[bl];
                 triangles[b] = triangles[ar];
 
@@ -536,9 +442,9 @@ shull2d(const torch::Tensor& points)
         hull.set(hull.key(points[ie].unsqueeze(0)), ie);
     }
 
-    const auto tn = static_cast<int64_t>(hull.triangles.size() / 3);
-    const auto opts = torch::dtype(torch::kInt64);
-    return at::from_blob(hull.triangles.data(), {tn, 3}, opts);
+    int64_t tn = hull.triangles.size() / 3;
+    auto answer = torch::tensor(hull.triangles, torch::TensorOptions().dtype(torch::kInt64));
+    return answer.reshape({tn, 3});
 }
 
 

torch_delaunay/__bind__/python_module.cc

@@ -6,7 +6,12 @@
 namespace torch_delaunay {
 
 
-PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) { m.def("shull2d", &shull2d); };
+PYBIND11_MODULE(TORCH_EXTENSION_NAME, m)
+{
+    m.def("shull2d", &shull2d);
+    m.def("circumcenter2d", &circumcenter2d);
+    m.def("circumradius2d", &circumradius2d);
+};
 
 
 } // namespace torch_delaunay

torch_delaunay/functional.py

@@ -1,6 +1,14 @@
-import torch
 import torch_delaunay._C as _C
+from torch import Tensor
 
 
-def shull2d(points: torch.Tensor) -> torch.Tensor:
+def shull2d(points: Tensor) -> Tensor:
     return _C.shull2d(points)
+
+
+def circumcenter2d(p0: Tensor, p1: Tensor, p2: Tensor) -> Tensor:
+    return _C.circumcenter2d(p0, p1, p2)
+
+
+def circumradius2d(p0: Tensor, p1: Tensor, p2: Tensor) -> Tensor:
+    return _C.circumradius2d(p0, p1, p2)

torch_delaunay/functional_test.py

@@ -1,11 +1,29 @@
 import torch
+from shapely import Point
+from geopandas import GeoDataFrame
+from geopandas import points_from_xy
 
-from torch_delaunay.functional import shull2d
+from torch_delaunay.functional import shull2d, circumcenter2d, circumradius2d
 
 
 def test_shull2d() -> None:
     points = torch.randn((1000, 2), dtype=torch.float64)
     simplices = shull2d(points)
 
+    assert (simplices >= 1000).sum() == 0, "simplices are outside of points array limit"
+
     assert simplices.shape[1] == 3
     assert simplices.shape[0] > 0
+
+    vertices = GeoDataFrame(geometry=points_from_xy(points[:, 0], points[:, 1]))
+
+    # Ensure that triangles comply with Delaunay definition.
+    tri = points[simplices]
+
+    centers = circumcenter2d(tri[:, 0], tri[:, 1], tri[:, 2]).detach().numpy()
+    radii = circumradius2d(tri[:, 0], tri[:, 1], tri[:, 2]).detach().numpy()
+
+    for center, radius in zip(centers, radii):
+        circle = Point(*center).buffer(radius)
+        incircle = vertices.intersection(circle)
+        assert (~incircle.is_empty).sum() <= 3