1 // Copyright (c) 2003,2006 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/Predicate_constructions_C2.h $ 7 // $Id: Predicate_constructions_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 14 15 #ifndef CGAL_APOLLONIUS_GRAPH_2_PREDICATE_CONSTRUCTIONS_2_H 16 #define CGAL_APOLLONIUS_GRAPH_2_PREDICATE_CONSTRUCTIONS_2_H 1 17 18 #include <CGAL/license/Apollonius_graph_2.h> 19 20 21 #include <CGAL/Apollonius_graph_2/basic.h> 22 23 namespace CGAL { 24 25 namespace ApolloniusGraph_2 { 26 27 template< class K > 28 class Inverted_weighted_point_2 29 : public K::Site_2 30 { 31 public: 32 typedef typename K::Site_2 K_Site_2; 33 typedef typename K::FT FT; 34 private: 35 FT _p; 36 public: Inverted_weighted_point_2(const K_Site_2 & wp,const FT & p)37 Inverted_weighted_point_2(const K_Site_2& wp, const FT& p) 38 : K_Site_2(wp), _p(p) {} 39 p()40 inline FT p() const { return _p; } 41 }; 42 43 44 template< class K > 45 class Weighted_point_inverter_2 46 { 47 public: 48 typedef typename K::Point_2 Point_2; 49 typedef typename K::Site_2 Site_2; 50 typedef Inverted_weighted_point_2<K> Inverted_weighted_point; 51 typedef typename K::FT FT; 52 private: 53 Site_2 _pole; 54 public: Weighted_point_inverter_2(const Site_2 & pole)55 Weighted_point_inverter_2(const Site_2& pole) 56 : _pole(pole) {} 57 operator()58 Inverted_weighted_point operator()(const Site_2& wp) 59 { 60 FT xs = wp.x() - _pole.x(); 61 FT ys = wp.y() - _pole.y(); 62 FT ws = wp.weight() - _pole.weight(); 63 FT ps = CGAL::square(xs) + CGAL::square(ys) 64 - CGAL::square(ws); 65 66 return 67 Inverted_weighted_point(Site_2(Point_2(xs, ys), ws), ps); 68 } 69 pole()70 Site_2 pole() const { return _pole; } 71 }; 72 73 74 template< class K > 75 class Voronoi_radius_2 76 { 77 // this class stores the coefficients for the tritangent circle 78 // radius equation. In particular we have: 79 // a x^2 - 2 b x + c = 0; 80 // x here represents the inverse of the radius 81 public: 82 typedef typename K::FT FT; 83 typedef Inverted_weighted_point_2<K> Inverted_weighted_point; 84 85 private: 86 FT _a, _b, _c; 87 FT _c2, _delta; 88 FT _dxp, _dyp, _dwp; 89 Voronoi_radius_2(FT a,FT b,FT c,FT c2,FT delta,FT dxp,FT dyp,FT dwp)90 Voronoi_radius_2(FT a, FT b, FT c, FT c2, FT delta, 91 FT dxp, FT dyp, FT dwp) 92 : _a(a), _b(b), _c(c), _c2(c2), _delta(delta), _dxp(dxp), 93 _dyp(dyp), _dwp(dwp) {} 94 95 public: Voronoi_radius_2(const Inverted_weighted_point & u1,const Inverted_weighted_point & u2)96 Voronoi_radius_2(const Inverted_weighted_point& u1, 97 const Inverted_weighted_point& u2) 98 { 99 FT dxp = determinant(u1.x(), u1.p(), u2.x(), u2.p()); 100 FT dyp = determinant(u1.y(), u1.p(), u2.y(), u2.p()); 101 FT dwp = determinant(u1.weight(), u1.p(), u2.weight(), u2.p()); 102 FT dxy = determinant(u1.x(), u1.y(), u2.x(), u2.y()); 103 FT dxw = determinant(u1.x(), u1.weight(), u2.x(), u2.weight()); 104 FT dyw = determinant(u1.y(), u1.weight(), u2.y(), u2.weight()); 105 106 _a = CGAL::square(dxp) + CGAL::square(dyp); 107 _b = dxp * dxw + dyp * dyw; 108 _c = CGAL::square(dxw) + CGAL::square(dyw) - CGAL::square(dxy); 109 _c2 = dxy; 110 _delta = _a - CGAL::square(dwp); 111 _dxp = dxp; 112 _dyp = dyp; 113 _dwp = dwp; 114 } 115 a()116 inline FT a() const { return _a; } b()117 inline FT b() const { return _b; } c()118 inline FT c() const { return _c; } 119 c1()120 inline FT c1() const { return _b; } c2()121 inline FT c2() const { return _c2; } delta()122 inline FT delta() const { return _delta; } d()123 inline FT d() const { return _a; } 124 dxp()125 inline FT dxp() const { return _dxp; } dyp()126 inline FT dyp() const { return _dyp; } dwp()127 inline FT dwp() const { return _dwp; } 128 is_first_root()129 inline bool is_first_root() const { return CGAL::is_negative(_c2); } 130 get_symmetric()131 Voronoi_radius_2 get_symmetric() 132 { 133 return Voronoi_radius_2(_a, _b, _c, -_c2, _delta, -_dxp, -_dyp, -_dwp); 134 } 135 }; 136 137 138 template< class K > 139 class Bitangent_line_2 140 { 141 // this class computes and stores the data for the left bitangent 142 // line of the weighted points p1, p2 oriented from p1 to p2 143 // or the left bitangent line of the inverted weighted point u1 and 144 // u2, oriented from u1 to u2 145 public: 146 typedef typename K::Point_2 Point_2; 147 typedef typename K::Site_2 Site_2; 148 typedef Inverted_weighted_point_2<K> Inverted_weighted_point; 149 typedef typename K::FT FT; 150 protected: 151 FT _a1, _a2; 152 FT _b1, _b2; 153 FT _c1, _c2; 154 FT _delta; 155 FT _d; 156 FT _dw; 157 FT _dxw, _dyw; 158 Bitangent_line_2(FT a1,FT a2,FT b1,FT b2,FT c1,FT c2,FT delta,FT d,FT dw,FT dxw,FT dyw)159 Bitangent_line_2(FT a1, FT a2, FT b1, FT b2, FT c1, FT c2, 160 FT delta, FT d, FT dw, FT dxw, FT dyw) 161 : _a1(a1), _a2(a2), _b1(b1), _b2(b2), _c1(c1), _c2(c2), 162 _delta(delta), _d(d), _dw(dw),_dxw(dxw), _dyw(dyw) {} 163 164 inline void store(FT dx,FT dy,FT dw)165 store(FT dx, FT dy, FT dw) 166 { 167 _dw = dw; 168 _a1 = dx * dw; 169 _a2 = dy; 170 _b1 = dy * dw; 171 _b2 = -dx; 172 } 173 174 inline void store(FT dx,FT dy,FT dw,FT dxy,FT dxw,FT dyw)175 store(FT dx, FT dy, FT dw, FT dxy, FT dxw, FT dyw) 176 { 177 store(dx, dy, dw); 178 _c1 = dx * dxw + dy * dyw; 179 _c2 = dxy; 180 _d = CGAL::square(dx) + CGAL::square(dy); 181 _delta = _d - CGAL::square(dw); 182 _dxw = dxw; 183 _dyw = dyw; 184 } 185 186 public: Bitangent_line_2(const Site_2 & p1,const Site_2 & p2)187 Bitangent_line_2(const Site_2& p1, const Site_2& p2) 188 { 189 FT dx = p1.x() - p2.x(); 190 FT dy = p1.y() - p2.y(); 191 FT dw = p1.weight() - p2.weight(); 192 FT dxy = determinant(p1.x(), p1.y(), p2.x(), p2.y()); 193 FT dxw = determinant(p1.x(), p1.weight(), p2.x(), p2.weight()); 194 FT dyw = determinant(p1.y(), p1.weight(), p2.y(), p2.weight()); 195 196 store(dx, dy, dw, dxy, dxw, dyw); 197 } 198 199 Bitangent_line_2(const Inverted_weighted_point & u1,const Inverted_weighted_point & u2)200 Bitangent_line_2(const Inverted_weighted_point& u1, 201 const Inverted_weighted_point& u2) 202 { 203 FT dxp = determinant(u1.x(), u1.p(), u2.x(), u2.p()); 204 FT dyp = determinant(u1.y(), u1.p(), u2.y(), u2.p()); 205 FT dwp = determinant(u1.weight(), u1.p(), u2.weight(), u2.p()); 206 FT dxy = determinant(u1.x(), u1.y(), u2.x(), u2.y()); 207 FT dxw = determinant(u1.x(), u1.weight(), u2.x(), u2.weight()); 208 FT dyw = determinant(u1.y(), u1.weight(), u2.y(), u2.weight()); 209 210 store(dxp, dyp, dwp, dxy, dxw, dyw); 211 } 212 get_symmetric()213 Bitangent_line_2 get_symmetric() const 214 { 215 return 216 Bitangent_line_2(_a1, -_a2, _b1, -_b2, _c1, -_c2, _delta, _d, 217 -_dw, -_dxw, -_dyw); 218 } 219 get_rot90()220 Bitangent_line_2 get_rot90() const 221 { 222 return 223 Bitangent_line_2(-_b1, -_b2, _a1, _a2, _c1, _c2, _delta, _d, 224 _dw, -_dyw, _dxw); 225 } 226 perpendicular(const Point_2 & p)227 Bitangent_line_2 perpendicular(const Point_2& p) const 228 { 229 // THIS DOES NOT KEEP TRACK OF THE ADDITIONALLY STORED 230 // QUANTITIES; THIS IS INEVITABLE IN ANY CASE SINCE GIVEN p WE 231 // CANNOT ANY LONGER HOPE TO KEEP TRACK OF THOSE 232 Bitangent_line_2 rotated = get_rot90(); 233 rotated._c1 = _b1 * p.x() - _a1 * p.y(); 234 rotated._c2 = _b2 * p.x() - _a2 * p.y(); 235 236 return rotated; 237 } 238 perpendicular(const Inverted_weighted_point & u)239 Bitangent_line_2 perpendicular(const Inverted_weighted_point& u) const 240 { 241 // THIS DOES NOT KEEP TRACK OF THE ADDITIONALLY STORED 242 // QUANTITIES; THIS IS INEVITABLE IN ANY CASE SINCE GIVEN p WE 243 // CANNOT ANY LONGER HOPE TO KEEP TRACK OF THOSE 244 Bitangent_line_2 rotated = get_rot90(); 245 rotated._c1 = (_b1 * u.x() - _a1 * u.y()) * u.p(); 246 rotated._c2 = (_b2 * u.x() - _a2 * u.y()) * u.p(); 247 248 return rotated; 249 } 250 a1()251 inline FT a1() const { return _a1; } a2()252 inline FT a2() const { return _a2; } b1()253 inline FT b1() const { return _b1; } b2()254 inline FT b2() const { return _b2; } c1()255 inline FT c1() const { return _c1; } c2()256 inline FT c2() const { return _c2; } delta()257 inline FT delta() const { return _delta; } d()258 inline FT d() const { return _d; } 259 dx()260 inline FT dx() const { return -_b2; } dy()261 inline FT dy() const { return _a2; } dw()262 inline FT dw() const { return _dw; } dxw()263 inline FT dxw() const { return _dxw; } dyw()264 inline FT dyw() const { return _dyw; } 265 }; 266 267 268 template< class K > 269 class Voronoi_circle_2 : public Bitangent_line_2<K> 270 { 271 public: 272 typedef Inverted_weighted_point_2<K> Inverted_weighted_point; 273 typedef Bitangent_line_2<K> Bitangent_line; 274 typedef Voronoi_radius_2<K> Voronoi_radius; 275 typedef typename Bitangent_line::FT FT; 276 277 protected: 278 FT _gamma; 279 280 inline compute_gamma()281 void compute_gamma() 282 { 283 _gamma = CGAL::square(this->_dxw) + CGAL::square(this->_dyw) 284 - CGAL::square(this->_c2); 285 } 286 287 public: Voronoi_circle_2(const Voronoi_radius & vr)288 Voronoi_circle_2(const Voronoi_radius& vr) 289 : Bitangent_line(FT(0), FT(0), FT(0), FT(0), vr.b(), vr.c2(), 290 vr.delta(), vr.d(), FT(0), FT(0), FT(0)), _gamma(vr.c()) 291 { 292 this->store(vr.dxp(), vr.dyp(), vr.dwp()); 293 } 294 Voronoi_circle_2(const Bitangent_line & bl)295 Voronoi_circle_2(const Bitangent_line& bl) 296 : Bitangent_line(bl.a1(), bl.a2(), bl.b1(), bl.b2(), bl.c1(), bl.c2(), 297 bl.delta(), bl.d(), bl.dw(), bl.dxw(), bl.dyw()) 298 { 299 compute_gamma(); 300 } 301 alpha()302 inline FT alpha() const { return this->_d; } beta()303 inline FT beta() const { return this->_c1; } gamma()304 inline FT gamma() const { return _gamma; } 305 is_first_root()306 inline bool is_first_root() const { 307 return CGAL::is_negative(this->_c2); 308 } 309 compute_P4(const Inverted_weighted_point & u1,const Inverted_weighted_point & u2,const Inverted_weighted_point & u3)310 FT compute_P4(const Inverted_weighted_point& u1, 311 const Inverted_weighted_point& u2, 312 const Inverted_weighted_point& u3) const 313 { 314 FT dx1 = determinant(u2.x(), u2.p(), u1.x(), u1.p()); 315 FT dy1 = determinant(u2.y(), u2.p(), u1.y(), u1.p()); 316 FT dw1 = determinant(u2.weight(), u2.p(), u1.weight(), u1.p()); 317 318 FT dx3 = determinant(u3.x(), u3.p(), u2.x(), u2.p()); 319 FT dy3 = determinant(u3.y(), u3.p(), u2.y(), u2.p()); 320 FT dw3 = determinant(u3.weight(), u3.p(), u2.weight(), u2.p()); 321 322 FT u2Pv2 = CGAL::square(u2.x()) + CGAL::square(u2.y()); 323 FT u2Mv2 = CGAL::square(u2.x()) - CGAL::square(u2.y()); 324 FT u2v2 = FT(2) * u2.x() * u2.y(); 325 326 327 FT vvMuu = dy1 * dy3 - dx1 * dx3; 328 FT vuPuv = dy1 * dx3 + dx1 * dy3; 329 330 FT dx2Pdy2_1 = CGAL::square(dx1) + CGAL::square(dy1); 331 FT dx2Pdy2_3 = CGAL::square(dx3) + CGAL::square(dy3); 332 333 FT fr1_sq = CGAL::square(dw1) * dx2Pdy2_3; 334 FT fr3_sq = CGAL::square(dw3) * dx2Pdy2_1; 335 336 FT f1 = (fr1_sq + fr3_sq) * CGAL::square(u2Pv2); 337 338 FT f2 = FT(2) * dw1 * dw3 * u2Pv2 * (u2Mv2 * vvMuu - u2v2 * vuPuv ); 339 FT f3 = CGAL::square(u2Mv2 * vuPuv + u2v2 * vvMuu); 340 341 FT F = f1 + f2 - f3; 342 343 344 FT uuPvv = dy1 * dy3 + dx1 * dx3; 345 FT vuMuv = dy1 * dx3 - dx1 * dy3; 346 347 FT G = fr1_sq + fr3_sq - FT(2) * dw1 * dw3 * uuPvv 348 - CGAL::square(vuMuv); 349 350 return (F * G); 351 } 352 }; 353 354 } //namespace ApolloniusGraph_2 355 356 } //namespace CGAL 357 358 #endif // CGAL_APOLLONIUS_GRAPH_2_PREDICATE_CONSTRUCTIONS_2_H 359