// Copyright (c) 2007 INRIA Sophia-Antipolis (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). // // $URL: https://github.com/CGAL/cgal/blob/v5.2/Apollonius_graph_2/include/CGAL/Apollonius_graph_2/Orientation8_C2.h $ // $Id: Orientation8_C2.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial // // // Author(s) : Menelaos Karavelas #ifndef CGAL_APOLLONIUS_GRAPH_2_ORIENTATION8_C2_H #define CGAL_APOLLONIUS_GRAPH_2_ORIENTATION8_C2_H #include #include #include //-------------------------------------------------------------------- namespace CGAL { namespace ApolloniusGraph_2 { template class Orientation8_C2 : public Orientation_2 { private: typedef Orientation_2 Base; public: typedef K Kernel; typedef MTag Method_tag; typedef typename K::Site_2 Site_2; typedef typename K::Point_2 Point_2; typedef typename K::Orientation Orientation; typedef typename K::FT FT; typedef Orientation result_type; typedef Site_2 argument_type; public: inline Orientation operator()(const Site_2& s1, const Site_2& s2, const Site_2& s3) const { return Kernel().orientation_2_object()(s1.point(), s2.point(), s3.point()); } inline Sign sqrt_ext_sign(const FT& A, const FT& B, const FT& Exp, const FT& Eyp, const FT& Erp, const FT& Exy2, const FT& Exr, const FT& Eyr, const FT dx, const FT& dy, const Field_with_sqrt_tag&) const { FT G = CGAL::square(Exp) + CGAL::square(Eyp) - CGAL::square(Erp); return CGAL::sign(A + B * CGAL::sqrt(G)); } inline Sign sqrt_ext_sign(const FT& A, const FT& B, const FT& Exp, const FT& Eyp, const FT& Erp, const FT& Exy2, const FT& Exr, const FT& Eyr, const FT dx, const FT& dy, const Integral_domain_without_division_tag&) const { Sign sA = CGAL::sign(A); Sign sB = CGAL::sign(B); if ( sA == CGAL::ZERO ) { return sB; } if ( sB == CGAL::ZERO ) { return sA; } if ( sA == sB ) { return sA; } Sign s = CGAL::sign(CGAL::square(Exy2 * Exr - Erp * dy) + CGAL::square(Exy2 * Eyr + Erp * dx) - CGAL::square(B)); return sA * s; } Orientation predicate(const Site_2& s1, const Site_2& s2, const Site_2& s3, const Site_2& p1, const Site_2& p2) const { // computes the orientation of the Voronoi vertex of s1, s2, s3 and // the points p1 and p2 FT xj = s2.x() - s1.x(); FT xk = s3.x() - s1.x(); FT xl = p1.x() - s1.x(); FT xm = p2.x() - s1.x(); FT yj = s2.y() - s1.y(); FT yk = s3.y() - s1.y(); FT yl = p1.y() - s1.y(); FT ym = p2.y() - s1.y(); FT dx = xl - xm; FT dy = yl - ym; FT rj = s2.weight() - s1.weight(); FT rk = s3.weight() - s1.weight(); FT pj = CGAL::square(xj) + CGAL::square(yj) - CGAL::square(rj); FT pk = CGAL::square(xk) + CGAL::square(yk) - CGAL::square(rk); FT Exp = determinant(xj, pj, xk, pk); FT Eyp = determinant(yj, pj, yk, pk); FT Erp = determinant(rj, pj, rk, pk); FT Exy = determinant(xj, yj, xk, yk); FT Exr = determinant(xj, rj, xk, rk); FT Eyr = determinant(yj, rj, yk, rk); FT Exy2 = 2 * determinant(xl, yl, xm, ym); FT A = (Exp * Exr + Eyp * Eyr) * Exy2 + (Eyp * dx - Exp * dy) * Erp; FT B = Exy * Exy2 - Exp * dx - Eyp * dy; return sqrt_ext_sign(A, B, Exp, Eyp, Erp, Exy2, Exr, Eyr, dx, dy, Method_tag()); } Orientation operator()(const Site_2& s1, const Site_2& s2, const Site_2& s3, const Site_2& p1, const Site_2& p2) const { Orientation o = predicate(s1, s2, s3, p1, p2); #ifndef NDEBUG Orientation o_old = Base::operator()(s1, s2, s3, p1, p2); CGAL_assertion( o == o_old ); #endif return o; } }; //-------------------------------------------------------------------- template class Constructive_orientation8_C2 { public: typedef K Kernel; typedef MTag Method_tag; typedef typename K::Site_2 Site_2; typedef typename K::Point_2 Point_2; typedef typename K::Orientation Orientation; typedef typename K::FT FT; typedef Orientation result_type; typedef Site_2 argument_type; private: FT s1x, s1y; FT xj, xk, yj, yk, rj, rk, nj, nk, pj, pk, Exp, Eyp, Erp, Exy, Exr, Eyr, A1; Orientation o_sym; public: Constructive_orientation8_C2(const Site_2& s1, const Site_2& s2, const Site_2& s3, bool use_xj) { s1x = s1.x(); s1y = s1.y(); xj = s2.x() - s1.x(); xk = s3.x() - s1.x(); yj = s2.y() - s1.y(); yk = s3.y() - s1.y(); rj = s2.weight() - s1.weight(); rk = s3.weight() - s1.weight(); nj = CGAL::square(xj) + CGAL::square(yj); nk = CGAL::square(xk) + CGAL::square(yk); pj = nj - CGAL::square(rj); pk = nk - CGAL::square(rk); Exp = determinant(xj, pj, xk, pk); Eyp = determinant(yj, pj, yk, pk); Erp = determinant(rj, pj, rk, pk); Exy = determinant(xj, yj, xk, yk); Exr = determinant(xj, rj, xk, rk); Eyr = determinant(yj, rj, yk, rk); A1 = Exp * Exr + Eyp * Eyr; FT A, B, norm; if ( use_xj ) { A = (-Eyp * xj + Exp * yj) * Erp; B = Exp * xj + Eyp * yj; norm = nj; } else { A = (-Eyp * xk + Exp * yk) * Erp; B = Exp * xk + Eyp * yk; norm = nk; } o_sym = sqrt_ext_sign(A, B, norm, Method_tag()); } private: inline Sign sqrt_ext_sign(const FT& A, const FT& B, const FT& Exy2, const FT& dx, const FT& dy, const Field_with_sqrt_tag&) const { FT G = CGAL::square(Exp) + CGAL::square(Eyp) - CGAL::square(Erp); return CGAL::sign(A + B * CGAL::sqrt(G)); } inline Sign sqrt_ext_sign(const FT& A, const FT& B, const FT&, const Field_with_sqrt_tag&) const { FT G = CGAL::square(Exp) + CGAL::square(Eyp) - CGAL::square(Erp); return CGAL::sign(A + B * CGAL::sqrt(G)); } inline Sign sqrt_ext_sign(const FT& A, const FT& B, const FT& norm, const Integral_domain_without_division_tag&) const { Sign sA = CGAL::sign(A); Sign sB = CGAL::sign(B); if ( sA == CGAL::ZERO ) { return sB; } if ( sB == CGAL::ZERO ) { return sA; } if ( sA == sB ) { return sA; } Sign s = CGAL::sign(CGAL::square(Erp) * norm - CGAL::square(B)); return sA * s; } inline Sign sqrt_ext_sign(const FT& A, const FT& B, const FT& Exy2, const FT dx, const FT& dy, const Integral_domain_without_division_tag&) const { Sign sA = CGAL::sign(A); Sign sB = CGAL::sign(B); if ( sA == CGAL::ZERO ) { return sB; } if ( sB == CGAL::ZERO ) { return sA; } if ( sA == sB ) { return sA; } Sign s = CGAL::sign(CGAL::square(Exy2 * Exr - Erp * dy) + CGAL::square(Exy2 * Eyr + Erp * dx) - CGAL::square(B)); return sA * s; } Orientation predicate(const Site_2& p1, const Site_2& p2) const { // computes the orientation of the Voronoi vertex of s1, s2, s3 and // the points p1 and p2 FT xl = p1.x() - s1x; FT xm = p2.x() - s1x; FT yl = p1.y() - s1y; FT ym = p2.y() - s1y; FT dx = xl - xm; FT dy = yl - ym; FT Exy2 = 2 * determinant(xl, yl, xm, ym); FT A = A1 * Exy2 + (Eyp * dx - Exp * dy) * Erp; FT B = Exy * Exy2 - Exp * dx - Eyp * dy; return sqrt_ext_sign(A, B, Exy2, dx, dy, Method_tag()); } public: inline Orientation operator()(const Site_2& s1, const Site_2& s2, const Site_2& s3) const { return Kernel().orientation_2_object()(s1.point(), s2.point(), s3.point()); } inline Orientation operator()(const Site_2& p1, const Site_2& p2) const { Orientation o = predicate(p1, p2); return o; } inline Orientation operator()() const { return o_sym; } }; //-------------------------------------------------------------------- } //namespace ApolloniusGraph_2 } //namespace CGAL #endif // CGAL_APOLLONIUS_GRAPH_2_ORIENTATION8_C2_H