diff --git a/pyproject.toml b/pyproject.toml index 90dc41d..31c4102 100644 --- a/pyproject.toml +++ b/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" ] diff --git a/src/triangle.cc b/src/triangle.cc index 07422cc..776fd7d 100644 --- a/src/triangle.cc +++ b/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() > 0)) { - std::cout << "flip " << i << std::endl; - out.index_put_({i, 0}, p0[i].item()); - out.index_put_({i, 1}, pr[i].item()); - out.index_put_({i, 2}, p1[i].item()); - } - } -} - - -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(); - auto v1 = triangles.index({i, 1}).item(); - auto v2 = triangles.index({i, 2}).item(); - - 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 hash; std::vector 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()) { + if (incircle2d(p0, pr, pl, p1).lt(0).all().item()) { 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(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}); } diff --git a/torch_delaunay/__bind__/python_module.cc b/torch_delaunay/__bind__/python_module.cc index ebd0b46..b1e7e02 100644 --- a/torch_delaunay/__bind__/python_module.cc +++ b/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 diff --git a/torch_delaunay/functional.py b/torch_delaunay/functional.py index 0e51f84..f4fb3bc 100644 --- a/torch_delaunay/functional.py +++ b/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) diff --git a/torch_delaunay/functional_test.py b/torch_delaunay/functional_test.py index caa1456..740a4a6 100644 --- a/torch_delaunay/functional_test.py +++ b/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