92 lines
3.2 KiB
C++
92 lines
3.2 KiB
C++
#pragma once
|
|
|
|
#include <psemek/geom/sign.hpp>
|
|
#include <psemek/geom/point.hpp>
|
|
#include <psemek/geom/orientation.hpp>
|
|
#include <psemek/geom/robust.hpp>
|
|
|
|
#ifdef PSEMEK_ROBUST_PREDICATES
|
|
#include <boost/multiprecision/gmp.hpp>
|
|
#endif
|
|
|
|
#include <type_traits>
|
|
|
|
namespace psemek::geom
|
|
{
|
|
|
|
template <typename T, std::size_t N, typename ... Points>
|
|
sign_t in_circle(fast_predicate_tag, point<T, N> const & p0, Points const & ... points)
|
|
{
|
|
auto proj = [](vector<T, N> const & v)
|
|
{
|
|
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;
|
|
};
|
|
|
|
auto result = orientation(fast, proj(points - p0) ...);
|
|
if constexpr ((N % 2) == 0)
|
|
return inverse(result);
|
|
else
|
|
return result;
|
|
}
|
|
|
|
template <typename T>
|
|
std::enable_if_t<std::is_floating_point_v<T>, sign_t>
|
|
in_circle(robust_predicate_tag, point<T, 2> const & p0, point<T, 2> const & p1, point<T, 2> const & p2, point<T, 2> const & p3)
|
|
{
|
|
#ifdef PSEMEK_ROBUST_PREDICATES
|
|
constexpr T error = std::numeric_limits<T>::epsilon() * T(29) / T(2);
|
|
|
|
T const m01 = (p0[0] - p3[0]) * (p1[1] - p3[1]);
|
|
T const m02 = (p0[0] - p3[0]) * (p2[1] - p3[1]);
|
|
T const m10 = (p1[0] - p3[0]) * (p0[1] - p3[1]);
|
|
T const m12 = (p1[0] - p3[0]) * (p2[1] - p3[1]);
|
|
T const m20 = (p2[0] - p3[0]) * (p0[1] - p3[1]);
|
|
T const m21 = (p2[0] - p3[0]) * (p1[1] - p3[1]);
|
|
|
|
T const d = T(0)
|
|
+ m01 * p2[0] * p2[0] + m01 * p2[1] * p2[1] - m01 * p3[0] * p3[0] - m01 * p3[1] * p3[1]
|
|
- m02 * p1[0] * p1[0] - m02 * p1[1] * p1[1] + m02 * p3[0] * p3[0] + m02 * p3[1] * p3[1]
|
|
- m10 * p2[0] * p2[0] - m10 * p2[1] * p2[1] + m10 * p3[0] * p3[0] + m10 * p3[1] * p3[1]
|
|
+ m12 * p0[0] * p0[0] + m12 * p0[1] * p0[1] - m12 * p3[0] * p3[0] - m12 * p3[1] * p3[1]
|
|
+ m20 * p1[0] * p1[0] + m20 * p1[1] * p1[1] - m20 * p3[0] * p3[0] - m20 * p3[1] * p3[1]
|
|
- m21 * p0[0] * p0[0] - m21 * p0[1] * p0[1] + m21 * p3[0] * p3[0] + m21 * p3[1] * p3[1]
|
|
;
|
|
|
|
T const t = T(0)
|
|
+ std::abs(m01 * p2[0] * p2[0]) + std::abs(m01 * p2[1] * p2[1]) + std::abs(m01 * p3[0] * p3[0]) + std::abs(m01 * p3[1] * p3[1])
|
|
+ std::abs(m02 * p1[0] * p1[0]) + std::abs(m02 * p1[1] * p1[1]) + std::abs(m02 * p3[0] * p3[0]) + std::abs(m02 * p3[1] * p3[1])
|
|
+ std::abs(m10 * p2[0] * p2[0]) + std::abs(m10 * p2[1] * p2[1]) + std::abs(m10 * p3[0] * p3[0]) + std::abs(m10 * p3[1] * p3[1])
|
|
+ std::abs(m12 * p0[0] * p0[0]) + std::abs(m12 * p0[1] * p0[1]) + std::abs(m12 * p3[0] * p3[0]) + std::abs(m12 * p3[1] * p3[1])
|
|
+ std::abs(m20 * p1[0] * p1[0]) + std::abs(m20 * p1[1] * p1[1]) + std::abs(m20 * p3[0] * p3[0]) + std::abs(m20 * p3[1] * p3[1])
|
|
+ std::abs(m21 * p0[0] * p0[0]) + std::abs(m21 * p0[1] * p0[1]) + std::abs(m21 * p3[0] * p3[0]) + std::abs(m21 * p3[1] * p3[1])
|
|
;
|
|
|
|
if (d > t * error)
|
|
return sign_t::positive;
|
|
else if (d < - t * error)
|
|
return sign_t::negative;
|
|
else
|
|
{
|
|
using exact_type = boost::multiprecision::mpq_rational;
|
|
|
|
return in_circle(cast<exact_type>(p0), cast<exact_type>(p1), cast<exact_type>(p2), cast<exact_type>(p3));
|
|
}
|
|
#else
|
|
return in_circle(fast_predicate_tag{}, p0, p1, p2, p3);
|
|
#endif
|
|
}
|
|
|
|
template <typename ... Args>
|
|
sign_t in_circle(Args const & ... args)
|
|
{
|
|
return in_circle(default_robust_tag, args...);
|
|
}
|
|
|
|
}
|