1 // Copyright (c) 2007 INRIA Sophia-Antipolis (France). 2 // All rights reserved. 3 // 4 // This file is part of CGAL (www.cgal.org). 5 // 6 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Apollonius_graph_2/include/CGAL/Apollonius_graph_2/Orientation8_C2.h $ 7 // $Id: Orientation8_C2.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot 8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial 9 // 10 // 11 // Author(s) : Menelaos Karavelas <mkaravel@iacm.forth.gr> 12 13 #ifndef CGAL_APOLLONIUS_GRAPH_2_ORIENTATION8_C2_H 14 #define CGAL_APOLLONIUS_GRAPH_2_ORIENTATION8_C2_H 15 16 #include <CGAL/license/Apollonius_graph_2.h> 17 18 19 #include <CGAL/determinant.h> 20 #include <CGAL/Apollonius_graph_2/Orientation_2.h> 21 22 //-------------------------------------------------------------------- 23 24 namespace CGAL { 25 26 namespace ApolloniusGraph_2 { 27 28 template<class K, class MTag> 29 class Orientation8_C2 30 : public Orientation_2<K,MTag> 31 { 32 private: 33 typedef Orientation_2<K,MTag> Base; 34 35 public: 36 typedef K Kernel; 37 typedef MTag Method_tag; 38 typedef typename K::Site_2 Site_2; 39 typedef typename K::Point_2 Point_2; 40 typedef typename K::Orientation Orientation; 41 typedef typename K::FT FT; 42 43 typedef Orientation result_type; 44 typedef Site_2 argument_type; 45 46 public: 47 inline operator()48 Orientation operator()(const Site_2& s1, const Site_2& s2, 49 const Site_2& s3) const 50 { 51 return Kernel().orientation_2_object()(s1.point(), s2.point(), 52 s3.point()); 53 } 54 55 inline 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 &)56 Sign sqrt_ext_sign(const FT& A, const FT& B, 57 const FT& Exp, const FT& Eyp, const FT& Erp, 58 const FT& Exy2, const FT& Exr, const FT& Eyr, 59 const FT dx, const FT& dy, 60 const Field_with_sqrt_tag&) const 61 { 62 FT G = CGAL::square(Exp) + CGAL::square(Eyp) - CGAL::square(Erp); 63 return CGAL::sign(A + B * CGAL::sqrt(G)); 64 } 65 66 inline 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 &)67 Sign sqrt_ext_sign(const FT& A, const FT& B, 68 const FT& Exp, const FT& Eyp, const FT& Erp, 69 const FT& Exy2, const FT& Exr, const FT& Eyr, 70 const FT dx, const FT& dy, 71 const Integral_domain_without_division_tag&) const 72 { 73 Sign sA = CGAL::sign(A); 74 Sign sB = CGAL::sign(B); 75 76 if ( sA == CGAL::ZERO ) { return sB; } 77 if ( sB == CGAL::ZERO ) { return sA; } 78 if ( sA == sB ) { return sA; } 79 80 Sign s = CGAL::sign(CGAL::square(Exy2 * Exr - Erp * dy) 81 + CGAL::square(Exy2 * Eyr + Erp * dx) 82 - CGAL::square(B)); 83 return sA * s; 84 } 85 predicate(const Site_2 & s1,const Site_2 & s2,const Site_2 & s3,const Site_2 & p1,const Site_2 & p2)86 Orientation predicate(const Site_2& s1, const Site_2& s2, 87 const Site_2& s3, const Site_2& p1, 88 const Site_2& p2) const 89 { 90 // computes the orientation of the Voronoi vertex of s1, s2, s3 and 91 // the points p1 and p2 92 FT xj = s2.x() - s1.x(); 93 FT xk = s3.x() - s1.x(); 94 95 FT xl = p1.x() - s1.x(); 96 FT xm = p2.x() - s1.x(); 97 98 99 FT yj = s2.y() - s1.y(); 100 FT yk = s3.y() - s1.y(); 101 102 FT yl = p1.y() - s1.y(); 103 FT ym = p2.y() - s1.y(); 104 105 106 FT dx = xl - xm; 107 FT dy = yl - ym; 108 109 110 FT rj = s2.weight() - s1.weight(); 111 FT rk = s3.weight() - s1.weight(); 112 113 FT pj = CGAL::square(xj) + CGAL::square(yj) - CGAL::square(rj); 114 FT pk = CGAL::square(xk) + CGAL::square(yk) - CGAL::square(rk); 115 116 FT Exp = determinant(xj, pj, xk, pk); 117 FT Eyp = determinant(yj, pj, yk, pk); 118 FT Erp = determinant(rj, pj, rk, pk); 119 120 FT Exy = determinant(xj, yj, xk, yk); 121 FT Exr = determinant(xj, rj, xk, rk); 122 FT Eyr = determinant(yj, rj, yk, rk); 123 124 FT Exy2 = 2 * determinant(xl, yl, xm, ym); 125 126 FT A = (Exp * Exr + Eyp * Eyr) * Exy2 + (Eyp * dx - Exp * dy) * Erp; 127 FT B = Exy * Exy2 - Exp * dx - Eyp * dy; 128 129 return sqrt_ext_sign(A, B, Exp, Eyp, Erp, 130 Exy2, Exr, Eyr, dx, dy, Method_tag()); 131 } 132 operator()133 Orientation operator()(const Site_2& s1, const Site_2& s2, 134 const Site_2& s3, const Site_2& p1, 135 const Site_2& p2) const 136 { 137 Orientation o = predicate(s1, s2, s3, p1, p2); 138 #ifndef NDEBUG 139 Orientation o_old = Base::operator()(s1, s2, s3, p1, p2); 140 141 CGAL_assertion( o == o_old ); 142 #endif 143 return o; 144 } 145 146 }; 147 148 //-------------------------------------------------------------------- 149 150 template<class K, class MTag> 151 class Constructive_orientation8_C2 152 { 153 public: 154 typedef K Kernel; 155 typedef MTag Method_tag; 156 typedef typename K::Site_2 Site_2; 157 typedef typename K::Point_2 Point_2; 158 typedef typename K::Orientation Orientation; 159 typedef typename K::FT FT; 160 161 typedef Orientation result_type; 162 typedef Site_2 argument_type; 163 164 private: 165 FT s1x, s1y; 166 FT xj, xk, yj, yk, rj, rk, nj, nk, pj, pk, Exp, Eyp, Erp, Exy, Exr, Eyr, A1; 167 Orientation o_sym; 168 169 public: Constructive_orientation8_C2(const Site_2 & s1,const Site_2 & s2,const Site_2 & s3,bool use_xj)170 Constructive_orientation8_C2(const Site_2& s1, const Site_2& s2, 171 const Site_2& s3, bool use_xj) 172 { 173 s1x = s1.x(); 174 s1y = s1.y(); 175 176 xj = s2.x() - s1.x(); 177 xk = s3.x() - s1.x(); 178 179 yj = s2.y() - s1.y(); 180 yk = s3.y() - s1.y(); 181 182 rj = s2.weight() - s1.weight(); 183 rk = s3.weight() - s1.weight(); 184 185 nj = CGAL::square(xj) + CGAL::square(yj); 186 nk = CGAL::square(xk) + CGAL::square(yk); 187 188 pj = nj - CGAL::square(rj); 189 pk = nk - CGAL::square(rk); 190 191 Exp = determinant(xj, pj, xk, pk); 192 Eyp = determinant(yj, pj, yk, pk); 193 Erp = determinant(rj, pj, rk, pk); 194 195 Exy = determinant(xj, yj, xk, yk); 196 Exr = determinant(xj, rj, xk, rk); 197 Eyr = determinant(yj, rj, yk, rk); 198 199 A1 = Exp * Exr + Eyp * Eyr; 200 201 FT A, B, norm; 202 if ( use_xj ) { 203 A = (-Eyp * xj + Exp * yj) * Erp; 204 B = Exp * xj + Eyp * yj; 205 norm = nj; 206 } else { 207 A = (-Eyp * xk + Exp * yk) * Erp; 208 B = Exp * xk + Eyp * yk; 209 norm = nk; 210 } 211 212 o_sym = sqrt_ext_sign(A, B, norm, Method_tag()); 213 } 214 215 private: 216 inline sqrt_ext_sign(const FT & A,const FT & B,const FT & Exy2,const FT & dx,const FT & dy,const Field_with_sqrt_tag &)217 Sign sqrt_ext_sign(const FT& A, const FT& B, const FT& Exy2, 218 const FT& dx, const FT& dy, 219 const Field_with_sqrt_tag&) const 220 { 221 FT G = CGAL::square(Exp) + CGAL::square(Eyp) - CGAL::square(Erp); 222 return CGAL::sign(A + B * CGAL::sqrt(G)); 223 } 224 225 inline sqrt_ext_sign(const FT & A,const FT & B,const FT &,const Field_with_sqrt_tag &)226 Sign sqrt_ext_sign(const FT& A, const FT& B, const FT&, 227 const Field_with_sqrt_tag&) const 228 { 229 FT G = CGAL::square(Exp) + CGAL::square(Eyp) - CGAL::square(Erp); 230 return CGAL::sign(A + B * CGAL::sqrt(G)); 231 } 232 233 234 inline sqrt_ext_sign(const FT & A,const FT & B,const FT & norm,const Integral_domain_without_division_tag &)235 Sign sqrt_ext_sign(const FT& A, const FT& B, const FT& norm, 236 const Integral_domain_without_division_tag&) const 237 { 238 Sign sA = CGAL::sign(A); 239 Sign sB = CGAL::sign(B); 240 241 if ( sA == CGAL::ZERO ) { return sB; } 242 if ( sB == CGAL::ZERO ) { return sA; } 243 if ( sA == sB ) { return sA; } 244 245 Sign s = CGAL::sign(CGAL::square(Erp) * norm - CGAL::square(B)); 246 return sA * s; 247 } 248 249 250 inline sqrt_ext_sign(const FT & A,const FT & B,const FT & Exy2,const FT dx,const FT & dy,const Integral_domain_without_division_tag &)251 Sign sqrt_ext_sign(const FT& A, const FT& B, const FT& Exy2, 252 const FT dx, const FT& dy, 253 const Integral_domain_without_division_tag&) const 254 { 255 Sign sA = CGAL::sign(A); 256 Sign sB = CGAL::sign(B); 257 258 if ( sA == CGAL::ZERO ) { return sB; } 259 if ( sB == CGAL::ZERO ) { return sA; } 260 if ( sA == sB ) { return sA; } 261 262 Sign s = CGAL::sign(CGAL::square(Exy2 * Exr - Erp * dy) 263 + CGAL::square(Exy2 * Eyr + Erp * dx) 264 - CGAL::square(B)); 265 return sA * s; 266 } 267 predicate(const Site_2 & p1,const Site_2 & p2)268 Orientation predicate(const Site_2& p1, const Site_2& p2) const 269 { 270 // computes the orientation of the Voronoi vertex of s1, s2, s3 and 271 // the points p1 and p2 272 FT xl = p1.x() - s1x; 273 FT xm = p2.x() - s1x; 274 275 FT yl = p1.y() - s1y; 276 FT ym = p2.y() - s1y; 277 278 FT dx = xl - xm; 279 FT dy = yl - ym; 280 281 FT Exy2 = 2 * determinant(xl, yl, xm, ym); 282 283 FT A = A1 * Exy2 + (Eyp * dx - Exp * dy) * Erp; 284 FT B = Exy * Exy2 - Exp * dx - Eyp * dy; 285 286 return sqrt_ext_sign(A, B, Exy2, dx, dy, Method_tag()); 287 } 288 289 public: 290 inline operator()291 Orientation operator()(const Site_2& s1, const Site_2& s2, 292 const Site_2& s3) const 293 { 294 return Kernel().orientation_2_object()(s1.point(), s2.point(), 295 s3.point()); 296 } 297 298 inline operator()299 Orientation operator()(const Site_2& p1, const Site_2& p2) const 300 { 301 Orientation o = predicate(p1, p2); 302 return o; 303 } 304 305 inline operator()306 Orientation operator()() const { 307 return o_sym; 308 } 309 }; 310 311 //-------------------------------------------------------------------- 312 313 314 } //namespace ApolloniusGraph_2 315 316 } //namespace CGAL 317 318 #endif // CGAL_APOLLONIUS_GRAPH_2_ORIENTATION8_C2_H 319