From ae7a08424b5ddcf4970bf6fdaa2ac4d3af0b3119 Mon Sep 17 00:00:00 2001 From: Menelaos Karavelas Date: Wed, 18 Jul 2007 21:46:47 +0000 Subject: [PATCH] code polishing; added constructive version of the predicate; added versions for different algebraic structures --- .../CGAL/Apollonius_graph_2/Orientation8_C2.h | 230 ++++++++++++++++-- 1 file changed, 216 insertions(+), 14 deletions(-) diff --git a/Apollonius_graph_2/include/CGAL/Apollonius_graph_2/Orientation8_C2.h b/Apollonius_graph_2/include/CGAL/Apollonius_graph_2/Orientation8_C2.h index 46813b78950..1cc66298105 100644 --- a/Apollonius_graph_2/include/CGAL/Apollonius_graph_2/Orientation8_C2.h +++ b/Apollonius_graph_2/include/CGAL/Apollonius_graph_2/Orientation8_C2.h @@ -57,6 +57,37 @@ public: 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 @@ -97,17 +128,147 @@ public: FT Exy2 = 2 * det2x2_by_formula(xl, yl, xm, ym); -#if 0 FT A = (Exp * Exr + Eyp * Eyr) * Exy2 + (Eyp * dx - Exp * dy) * Erp; FT B = Exy * Exy2 - Exp * dx - Eyp * dy; - FT C = CGAL::square(Exp) + CGAL::square(Eyp) - CGAL::square(Erp); - return sign_a_plus_b_x_sqrt_c(A, B, C); -#else - Sign sA = CGAL::sign((Exp * Exr + Eyp * Eyr) * Exy2 - + (Eyp * dx - Exp * dy) * Erp); + return sqrt_ext_sign(A, B, Exp, Eyp, Erp, + Exy2, Exr, Eyr, dx, dy, Method_tag()); + } - FT B = Exy * Exy2 - Exp * dx - Eyp * dy; + 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 Arity_tag<3> Arity; + typedef Site_2 argument_type; + +private: + FT s1x, s1y; + FT xj, xk, yj, yk, rj, rk, pj, pk, Exp, Eyp, Erp, Exy, Exr, Eyr, A1; + + FT xl_, yl_, ca1, ca2, ca3, cb1, cb2, cb3; + +public: + Constructive_orientation8_C2(const Site_2& s1, const Site_2& s2, + const Site_2& s3) + { + 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(); + + pj = CGAL::square(xj) + CGAL::square(yj) - CGAL::square(rj); + pk = CGAL::square(xk) + CGAL::square(yk) - CGAL::square(rk); + + Exp = det2x2_by_formula(xj, pj, xk, pk); + Eyp = det2x2_by_formula(yj, pj, yk, pk); + Erp = det2x2_by_formula(rj, pj, rk, pk); + + Exy = det2x2_by_formula(xj, yj, xk, yk); + Exr = det2x2_by_formula(xj, rj, xk, rk); + Eyr = det2x2_by_formula(yj, rj, yk, rk); + + A1 = Exp * Exr + Eyp * Eyr; + } + + Constructive_orientation8_C2(const Site_2& s1, const Site_2& s2, + const Site_2& s3, const Site_2& p1) + { + 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(); + + pj = CGAL::square(xj) + CGAL::square(yj) - CGAL::square(rj); + pk = CGAL::square(xk) + CGAL::square(yk) - CGAL::square(rk); + + Exp = det2x2_by_formula(xj, pj, xk, pk); + Eyp = det2x2_by_formula(yj, pj, yk, pk); + Erp = det2x2_by_formula(rj, pj, rk, pk); + + Exy = det2x2_by_formula(xj, yj, xk, yk); + Exr = det2x2_by_formula(xj, rj, xk, rk); + Eyr = det2x2_by_formula(yj, rj, yk, rk); + + A1 = Exp * Exr + Eyp * Eyr; + + xl_ = p1.x() - s1.x(); + yl_ = p1.y() - s1.y(); + + ca1 = 2 * A1 * xl_ - Eyp * Erp; + ca2 = -2 * A1 * yl_ + Exp * Erp; + ca3 = (Eyp * xl_ - Exp * yl_) * Erp; + + cb1 = Exp + Exy * xl_; + cb2 = Eyp - Exy * yl_; + cb3 = -Exp * xl_ - Eyp * yl_; + } + + + +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& 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& 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; } @@ -118,17 +279,57 @@ public: + CGAL::square(Exy2 * Eyr + Erp * dx) - CGAL::square(B)); return sA * s; -#endif } - Orientation operator()(const Site_2& s1, const Site_2& s2, - const Site_2& s3, const Site_2& p1, - const Site_2& p2) const + Orientation predicate(const Site_2& p1, const Site_2& p2) const { - Orientation o = predicate(s1, s2, s3, p1, p2); - Orientation o_old = Base::operator()(s1, s2, s3, p1, p2); + // 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; - CGAL_assertion( o == o_old ); + FT yl = p1.y() - s1y; + FT ym = p2.y() - s1y; + + FT dx = xl - xm; + FT dy = yl - ym; + + FT Exy2 = 2 * det2x2_by_formula(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()); + } + + Orientation predicate(const Site_2& p2) const + { + // computes the orientation of the Voronoi vertex of s1, s2, s3 and + // the points p1 and p2 + FT xm = p2.x() - s1x; + FT ym = p2.y() - s1y; + + FT dx = xl_ - xm; + FT dy = yl_ - ym; + + FT Exy2 = 2 * det2x2_by_formula(xl_, yl_, xm, ym); + + FT A = ca1 * xm + ca2 * ym + ca3; + FT B = cb1 * xm + cb2 * ym + cb3; + + return sqrt_ext_sign(A, B, Exy2, dx, dy, Method_tag()); + } + + + Orientation operator()(const Site_2& p1, const Site_2& p2) const + { + Orientation o = predicate(p1, p2); + return o; + } + + Orientation operator()(const Site_2& p2) const + { + Orientation o = predicate(p2); return o; } @@ -136,6 +337,7 @@ public: //-------------------------------------------------------------------- + CGAL_APOLLONIUS_GRAPH_2_END_NAMESPACE CGAL_END_NAMESPACE