Apply optimizations

benchmark/sweephull_benchmark.cc

@@ -21,6 +21,18 @@
 using namespace torch_delaunay;
 
 
+static void
+benchmark_shull2d_uniform_float64(benchmark::State& state)
+{
+    auto options = torch::TensorOptions().dtype(torch::kFloat64).device(torch::kCPU);
+    auto points = torch::rand({state.range(0), 2}, options);
+
+    for (auto _ : state) {
+        shull2d(points);
+    }
+}
+
+
 static void
 benchmark_shull2d_uniform_float32(benchmark::State& state)
 {
@@ -33,10 +45,17 @@ benchmark_shull2d_uniform_float32(benchmark::State& state)
 }
 
 
+BENCHMARK(benchmark_shull2d_uniform_float64)
+    ->Arg(1000)
+    ->Arg(100000)
+    ->Arg(1000000)
+    ->Unit(benchmark::kMillisecond);
+
+
 BENCHMARK(benchmark_shull2d_uniform_float32)
     ->Arg(1000)
-    ->Arg(10000)
     ->Arg(100000)
+    ->Arg(1000000)
     ->Unit(benchmark::kMillisecond);
 
 

include/torch_delaunay/predicate.h

@@ -15,8 +15,11 @@
 
 #pragma once
 
+#include <ATen/TensorAccessor.h>
 #include <torch/torch.h>
 
+#include <torch_delaunay/predicate.h>
+
 
 namespace torch_delaunay {
 
@@ -25,6 +28,22 @@ torch::Tensor
 orient2d(const torch::Tensor& p0, const torch::Tensor& p1, const torch::Tensor& p2);
 
 
+template <typename scalar_t>
+scalar_t
+orient2d_kernel(
+    const at::TensorAccessor<scalar_t, 1>& p0,
+    const at::TensorAccessor<scalar_t, 1>& p1,
+    const at::TensorAccessor<scalar_t, 1>& p2
+)
+{
+    const auto d0_x = p0[0] - p2[0];
+    const auto d0_y = p0[1] - p2[1];
+    const auto d1_x = p1[0] - p2[0];
+    const auto d1_y = p1[1] - p2[1];
+    return d0_x * d1_y - d0_y * d1_x;
+}
+
+
 torch::Tensor
 incircle2d(
     const torch::Tensor& p0,
@@ -34,8 +53,35 @@ incircle2d(
 );
 
 
-bool
+torch::Tensor
 incircle2d(const torch::Tensor& points, const torch::Tensor& q);
 
 
+template <typename scalar_t>
+scalar_t
+incircle2d_kernel(
+    const at::TensorAccessor<scalar_t, 1>& p0,
+    const at::TensorAccessor<scalar_t, 1>& p1,
+    const at::TensorAccessor<scalar_t, 1>& p2,
+    const at::TensorAccessor<scalar_t, 1>& q
+)
+{
+    const auto d0_x = p0[0] - q[0];
+    const auto d1_x = p1[0] - q[0];
+    const auto d2_x = p2[0] - q[0];
+    const auto d0_y = p0[1] - q[1];
+    const auto d1_y = p1[1] - q[1];
+    const auto d2_y = p2[1] - q[1];
+
+    const auto d0_s = d0_x * d0_x + d0_y * d0_y;
+    const auto d1_s = d1_x * d1_x + d1_y * d1_y;
+    const auto d2_s = d2_x * d2_x + d2_y * d2_y;
+
+    return (
+        d0_x * (d1_y * d2_s - d1_s * d2_y) - d0_y * (d1_x * d2_s - d1_s * d2_x)
+        + d0_s * (d1_x * d2_y - d1_y * d2_x)
+    );
+}
+
+
 } // namespace torch_delaunay

include/torch_delaunay/sweephull.h

@@ -17,14 +17,15 @@
 
 #include <optional>
 
+#include <ATen/core/Scalar.h>
 #include <torch/torch.h>
 
 
 namespace torch_delaunay {
 
 
 torch::Tensor
-shull2d(const torch::Tensor& points, std::optional<double> eps = std::nullopt);
+shull2d(const torch::Tensor& points, std::optional<const at::Scalar> eps = std::nullopt);
 
 
 } // namespace torch_delaunay

include/torch_delaunay/triangle.h

@@ -15,6 +15,8 @@
 
 #pragma once
 
+#include <cmath>
+
 #include <torch/torch.h>
 
 
@@ -32,8 +34,45 @@ circumcenter2d(const torch::Tensor& points);
 torch::Tensor
 circumradius2d(const torch::Tensor& p0, const torch::Tensor& p1, const torch::Tensor& p2);
 
+
 torch::Tensor
 circumradius2d(const torch::Tensor& points);
 
 
+template <typename scalar_t>
+std::tuple<scalar_t, scalar_t>
+circumcircle2d_kernel(
+    const at::TensorAccessor<scalar_t, 1> p0,
+    const at::TensorAccessor<scalar_t, 1> p1,
+    const at::TensorAccessor<scalar_t, 1> p2
+)
+{
+    const auto ax = p0[0], ay = p0[1];
+    const auto bx = p1[0] - p0[0], by = p1[1] - p0[1];
+    const auto cx = p2[0] - p0[0], cy = p2[1] - p0[1];
+
+    const auto d = 2 * (bx * cy - by * cx);
+
+    const auto b_norm = bx * bx + by * by;
+    const auto c_norm = cx * cx + cy * cy;
+
+    const auto ux = (cy * b_norm - by * c_norm) / d;
+    const auto uy = (bx * c_norm - cx * b_norm) / d;
+
+    return std::forward_as_tuple(ux, uy);
+}
+
+
+template <typename scalar_t>
+scalar_t
+circumradius2d_kernel(
+    const at::TensorAccessor<scalar_t, 1> p0,
+    const at::TensorAccessor<scalar_t, 1> p1,
+    const at::TensorAccessor<scalar_t, 1> p2
+)
+{
+    const auto [ux, uy] = circumcircle2d_kernel<scalar_t>(p0, p1, p2);
+    return std::sqrt(ux * ux + uy * uy);
+}
+
 } // namespace torch_delaunay

src/predicate.cc

@@ -61,14 +61,14 @@ incircle2d(
 }
 
 
-bool
+torch::Tensor
 incircle2d(const torch::Tensor& points, const torch::Tensor& q)
 {
     const auto d = points - q;
     const auto abc = d.square().sum(-1).view({-1, 1});
     const auto A = torch::hstack({d, abc});
 
-    return torch::linalg::det(A).sign().lt(0).all().item<bool>();
+    return torch::linalg::det(A).sign();
 }