From 46458abd71d57de4a1adf1ba33517430f74d84d2 Mon Sep 17 00:00:00 2001
From: lisyarus <lisyarus@gmail.com>
Date: Fri, 16 Sep 2022 17:56:11 +0300
Subject: [PATCH] Improve the non-robust geom::in_circle predicate and
 generalize to any directions

---
 libs/geom/include/psemek/geom/incircle.hpp | 23 +++++++++++++++-------
 1 file changed, 16 insertions(+), 7 deletions(-)

diff --git a/libs/geom/include/psemek/geom/incircle.hpp b/libs/geom/include/psemek/geom/incircle.hpp
index 92e5ed24..ed152fda 100644
--- a/libs/geom/include/psemek/geom/incircle.hpp
+++ b/libs/geom/include/psemek/geom/incircle.hpp
@@ -13,22 +13,31 @@
 namespace psemek::geom
 {
 
-	template <typename T>
+	template <typename T, std::size_t N, typename ... Points>
 #ifdef PSEMEK_GEOM_ROBUST_PREDICATES
 	std::enable_if_t<!std::is_floating_point_v<T>, sign_t>
 #else
 	sign_t
 #endif
-		in_circle(point<T, 2> const & p0, point<T, 2> const & p1, point<T, 2> const & p2, point<T, 2> const & p3)
+		in_circle(point<T, N> const & p0, Points const & ... points)
 	{
-		auto proj = [](point<T, 2> const & p)
+		auto proj = [](vector<T, N> const & v)
 		{
-			auto const x = p[0];
-			auto const y = p[1];
-			return point<T, 3>{ x, y, x*x + y*y };
+			vector<T, N + 1> u;
+			u[N] = T{0};
+			for (std::size_t i = 0; i < N; ++i)
+			{
+				u[i] = v[i];
+				u[N] += v[i] * v[i];
+			}
+			return u;
 		};
 
-		return orientation(proj(p0), proj(p1), proj(p2), proj(p3));
+		auto result = orientation(proj(points - p0) ...);
+		if constexpr ((N % 2) == 0)
+			return inverse(result);
+		else
+			return result;
 	}
 
 #ifdef PSEMEK_GEOM_ROBUST_PREDICATES