1 // Boost.Polygon library detail/voronoi_predicates.hpp header file 2 3 // Copyright Andrii Sydorchuk 2010-2012. 4 // Distributed under the Boost Software License, Version 1.0. 5 // (See accompanying file LICENSE_1_0.txt or copy at 6 // http://www.boost.org/LICENSE_1_0.txt) 7 8 // See http://www.boost.org for updates, documentation, and revision history. 9 10 #ifndef BOOST_POLYGON_DETAIL_VORONOI_PREDICATES 11 #define BOOST_POLYGON_DETAIL_VORONOI_PREDICATES 12 13 #include <utility> 14 15 #include "voronoi_robust_fpt.hpp" 16 17 namespace boost { 18 namespace polygon { 19 namespace detail { 20 21 // Predicate utilities. Operates with the coordinate types that could 22 // be converted to the 32-bit signed integer without precision loss. 23 template <typename CTYPE_TRAITS> 24 class voronoi_predicates { 25 public: 26 typedef typename CTYPE_TRAITS::int_type int_type; 27 typedef typename CTYPE_TRAITS::int_x2_type int_x2_type; 28 typedef typename CTYPE_TRAITS::uint_x2_type uint_x2_type; 29 typedef typename CTYPE_TRAITS::big_int_type big_int_type; 30 typedef typename CTYPE_TRAITS::fpt_type fpt_type; 31 typedef typename CTYPE_TRAITS::efpt_type efpt_type; 32 typedef typename CTYPE_TRAITS::ulp_cmp_type ulp_cmp_type; 33 typedef typename CTYPE_TRAITS::to_fpt_converter_type to_fpt_converter; 34 typedef typename CTYPE_TRAITS::to_efpt_converter_type to_efpt_converter; 35 36 enum { 37 ULPS = 64, 38 ULPSx2 = 128 39 }; 40 41 template <typename Point> is_vertical(const Point & point1,const Point & point2)42 static bool is_vertical(const Point& point1, const Point& point2) { 43 return point1.x() == point2.x(); 44 } 45 46 template <typename Site> is_vertical(const Site & site)47 static bool is_vertical(const Site& site) { 48 return is_vertical(site.point0(), site.point1()); 49 } 50 51 // Compute robust cross_product: a1 * b2 - b1 * a2. 52 // It was mathematically proven that the result is correct 53 // with epsilon relative error equal to 1EPS. robust_cross_product(int_x2_type a1_,int_x2_type b1_,int_x2_type a2_,int_x2_type b2_)54 static fpt_type robust_cross_product(int_x2_type a1_, 55 int_x2_type b1_, 56 int_x2_type a2_, 57 int_x2_type b2_) { 58 static to_fpt_converter to_fpt; 59 uint_x2_type a1 = static_cast<uint_x2_type>(is_neg(a1_) ? -a1_ : a1_); 60 uint_x2_type b1 = static_cast<uint_x2_type>(is_neg(b1_) ? -b1_ : b1_); 61 uint_x2_type a2 = static_cast<uint_x2_type>(is_neg(a2_) ? -a2_ : a2_); 62 uint_x2_type b2 = static_cast<uint_x2_type>(is_neg(b2_) ? -b2_ : b2_); 63 64 uint_x2_type l = a1 * b2; 65 uint_x2_type r = b1 * a2; 66 67 if (is_neg(a1_) ^ is_neg(b2_)) { 68 if (is_neg(a2_) ^ is_neg(b1_)) 69 return (l > r) ? -to_fpt(l - r) : to_fpt(r - l); 70 else 71 return -to_fpt(l + r); 72 } else { 73 if (is_neg(a2_) ^ is_neg(b1_)) 74 return to_fpt(l + r); 75 else 76 return (l < r) ? -to_fpt(r - l) : to_fpt(l - r); 77 } 78 } 79 80 typedef struct orientation_test { 81 public: 82 // Represents orientation test result. 83 enum Orientation { 84 RIGHT = -1, 85 COLLINEAR = 0, 86 LEFT = 1 87 }; 88 89 // Value is a determinant of two vectors (e.g. x1 * y2 - x2 * y1). 90 // Return orientation based on the sign of the determinant. 91 template <typename T> evalboost::polygon::detail::voronoi_predicates::orientation_test92 static Orientation eval(T value) { 93 if (is_zero(value)) return COLLINEAR; 94 return (is_neg(value)) ? RIGHT : LEFT; 95 } 96 evalboost::polygon::detail::voronoi_predicates::orientation_test97 static Orientation eval(int_x2_type dif_x1_, 98 int_x2_type dif_y1_, 99 int_x2_type dif_x2_, 100 int_x2_type dif_y2_) { 101 return eval(robust_cross_product(dif_x1_, dif_y1_, dif_x2_, dif_y2_)); 102 } 103 104 template <typename Point> evalboost::polygon::detail::voronoi_predicates::orientation_test105 static Orientation eval(const Point& point1, 106 const Point& point2, 107 const Point& point3) { 108 int_x2_type dx1 = static_cast<int_x2_type>(point1.x()) - 109 static_cast<int_x2_type>(point2.x()); 110 int_x2_type dx2 = static_cast<int_x2_type>(point2.x()) - 111 static_cast<int_x2_type>(point3.x()); 112 int_x2_type dy1 = static_cast<int_x2_type>(point1.y()) - 113 static_cast<int_x2_type>(point2.y()); 114 int_x2_type dy2 = static_cast<int_x2_type>(point2.y()) - 115 static_cast<int_x2_type>(point3.y()); 116 return eval(robust_cross_product(dx1, dy1, dx2, dy2)); 117 } 118 } ot; 119 120 template <typename Point> 121 class point_comparison_predicate { 122 public: 123 typedef Point point_type; 124 operator ()(const point_type & lhs,const point_type & rhs) const125 bool operator()(const point_type& lhs, const point_type& rhs) const { 126 if (lhs.x() == rhs.x()) 127 return lhs.y() < rhs.y(); 128 return lhs.x() < rhs.x(); 129 } 130 }; 131 132 template <typename Site, typename Circle> 133 class event_comparison_predicate { 134 public: 135 typedef Site site_type; 136 typedef Circle circle_type; 137 operator ()(const site_type & lhs,const site_type & rhs) const138 bool operator()(const site_type& lhs, const site_type& rhs) const { 139 if (lhs.x0() != rhs.x0()) 140 return lhs.x0() < rhs.x0(); 141 if (!lhs.is_segment()) { 142 if (!rhs.is_segment()) 143 return lhs.y0() < rhs.y0(); 144 if (is_vertical(rhs)) 145 return lhs.y0() <= rhs.y0(); 146 return true; 147 } else { 148 if (is_vertical(rhs)) { 149 if (is_vertical(lhs)) 150 return lhs.y0() < rhs.y0(); 151 return false; 152 } 153 if (is_vertical(lhs)) 154 return true; 155 if (lhs.y0() != rhs.y0()) 156 return lhs.y0() < rhs.y0(); 157 return ot::eval(lhs.point1(), lhs.point0(), rhs.point1()) == ot::LEFT; 158 } 159 } 160 operator ()(const site_type & lhs,const circle_type & rhs) const161 bool operator()(const site_type& lhs, const circle_type& rhs) const { 162 typename ulp_cmp_type::Result xCmp = 163 ulp_cmp(to_fpt(lhs.x0()), to_fpt(rhs.lower_x()), ULPS); 164 return xCmp == ulp_cmp_type::LESS; 165 } 166 operator ()(const circle_type & lhs,const site_type & rhs) const167 bool operator()(const circle_type& lhs, const site_type& rhs) const { 168 typename ulp_cmp_type::Result xCmp = 169 ulp_cmp(to_fpt(lhs.lower_x()), to_fpt(rhs.x0()), ULPS); 170 return xCmp == ulp_cmp_type::LESS; 171 } 172 operator ()(const circle_type & lhs,const circle_type & rhs) const173 bool operator()(const circle_type& lhs, const circle_type& rhs) const { 174 if (lhs.lower_x() != rhs.lower_x()) { 175 return lhs.lower_x() < rhs.lower_x(); 176 } 177 return lhs.y() < rhs.y(); 178 } 179 180 private: 181 ulp_cmp_type ulp_cmp; 182 to_fpt_converter to_fpt; 183 }; 184 185 template <typename Site> 186 class distance_predicate { 187 public: 188 typedef Site site_type; 189 typedef typename site_type::point_type point_type; 190 191 // Returns true if a horizontal line going through a new site intersects 192 // right arc at first, else returns false. If horizontal line goes 193 // through intersection point of the given two arcs returns false also. operator ()(const site_type & left_site,const site_type & right_site,const point_type & new_point) const194 bool operator()(const site_type& left_site, 195 const site_type& right_site, 196 const point_type& new_point) const { 197 if (!left_site.is_segment()) { 198 if (!right_site.is_segment()) { 199 return pp(left_site, right_site, new_point); 200 } else { 201 return ps(left_site, right_site, new_point, false); 202 } 203 } else { 204 if (!right_site.is_segment()) { 205 return ps(right_site, left_site, new_point, true); 206 } else { 207 return ss(left_site, right_site, new_point); 208 } 209 } 210 } 211 212 private: 213 // Represents the result of the epsilon robust predicate. If the 214 // result is undefined some further processing is usually required. 215 enum kPredicateResult { 216 LESS = -1, 217 UNDEFINED = 0, 218 MORE = 1 219 }; 220 221 // Robust predicate, avoids using high-precision libraries. 222 // Returns true if a horizontal line going through the new point site 223 // intersects right arc at first, else returns false. If horizontal line 224 // goes through intersection point of the given two arcs returns false. pp(const site_type & left_site,const site_type & right_site,const point_type & new_point) const225 bool pp(const site_type& left_site, 226 const site_type& right_site, 227 const point_type& new_point) const { 228 const point_type& left_point = left_site.point0(); 229 const point_type& right_point = right_site.point0(); 230 if (left_point.x() > right_point.x()) { 231 if (new_point.y() <= left_point.y()) 232 return false; 233 } else if (left_point.x() < right_point.x()) { 234 if (new_point.y() >= right_point.y()) 235 return true; 236 } else { 237 return static_cast<int_x2_type>(left_point.y()) + 238 static_cast<int_x2_type>(right_point.y()) < 239 static_cast<int_x2_type>(new_point.y()) * 2; 240 } 241 242 fpt_type dist1 = find_distance_to_point_arc(left_site, new_point); 243 fpt_type dist2 = find_distance_to_point_arc(right_site, new_point); 244 245 // The undefined ulp range is equal to 3EPS + 3EPS <= 6ULP. 246 return dist1 < dist2; 247 } 248 ps(const site_type & left_site,const site_type & right_site,const point_type & new_point,bool reverse_order) const249 bool ps(const site_type& left_site, const site_type& right_site, 250 const point_type& new_point, bool reverse_order) const { 251 kPredicateResult fast_res = fast_ps( 252 left_site, right_site, new_point, reverse_order); 253 if (fast_res != UNDEFINED) { 254 return fast_res == LESS; 255 } 256 257 fpt_type dist1 = find_distance_to_point_arc(left_site, new_point); 258 fpt_type dist2 = find_distance_to_segment_arc(right_site, new_point); 259 260 // The undefined ulp range is equal to 3EPS + 7EPS <= 10ULP. 261 return reverse_order ^ (dist1 < dist2); 262 } 263 ss(const site_type & left_site,const site_type & right_site,const point_type & new_point) const264 bool ss(const site_type& left_site, 265 const site_type& right_site, 266 const point_type& new_point) const { 267 // Handle temporary segment sites. 268 if (left_site.sorted_index() == right_site.sorted_index()) { 269 return ot::eval( 270 left_site.point0(), left_site.point1(), new_point) == ot::LEFT; 271 } 272 273 fpt_type dist1 = find_distance_to_segment_arc(left_site, new_point); 274 fpt_type dist2 = find_distance_to_segment_arc(right_site, new_point); 275 276 // The undefined ulp range is equal to 7EPS + 7EPS <= 14ULP. 277 return dist1 < dist2; 278 } 279 find_distance_to_point_arc(const site_type & site,const point_type & point) const280 fpt_type find_distance_to_point_arc( 281 const site_type& site, const point_type& point) const { 282 fpt_type dx = to_fpt(site.x()) - to_fpt(point.x()); 283 fpt_type dy = to_fpt(site.y()) - to_fpt(point.y()); 284 // The relative error is at most 3EPS. 285 return (dx * dx + dy * dy) / (to_fpt(2.0) * dx); 286 } 287 find_distance_to_segment_arc(const site_type & site,const point_type & point) const288 fpt_type find_distance_to_segment_arc( 289 const site_type& site, const point_type& point) const { 290 if (is_vertical(site)) { 291 return (to_fpt(site.x()) - to_fpt(point.x())) * to_fpt(0.5); 292 } else { 293 const point_type& segment0 = site.point0(); 294 const point_type& segment1 = site.point1(); 295 fpt_type a1 = to_fpt(segment1.x()) - to_fpt(segment0.x()); 296 fpt_type b1 = to_fpt(segment1.y()) - to_fpt(segment0.y()); 297 fpt_type k = get_sqrt(a1 * a1 + b1 * b1); 298 // Avoid subtraction while computing k. 299 if (!is_neg(b1)) { 300 k = to_fpt(1.0) / (b1 + k); 301 } else { 302 k = (k - b1) / (a1 * a1); 303 } 304 // The relative error is at most 7EPS. 305 return k * robust_cross_product( 306 static_cast<int_x2_type>(segment1.x()) - 307 static_cast<int_x2_type>(segment0.x()), 308 static_cast<int_x2_type>(segment1.y()) - 309 static_cast<int_x2_type>(segment0.y()), 310 static_cast<int_x2_type>(point.x()) - 311 static_cast<int_x2_type>(segment0.x()), 312 static_cast<int_x2_type>(point.y()) - 313 static_cast<int_x2_type>(segment0.y())); 314 } 315 } 316 fast_ps(const site_type & left_site,const site_type & right_site,const point_type & new_point,bool reverse_order) const317 kPredicateResult fast_ps( 318 const site_type& left_site, const site_type& right_site, 319 const point_type& new_point, bool reverse_order) const { 320 const point_type& site_point = left_site.point0(); 321 const point_type& segment_start = right_site.point0(); 322 const point_type& segment_end = right_site.point1(); 323 324 if (ot::eval(segment_start, segment_end, new_point) != ot::RIGHT) 325 return (!right_site.is_inverse()) ? LESS : MORE; 326 327 fpt_type dif_x = to_fpt(new_point.x()) - to_fpt(site_point.x()); 328 fpt_type dif_y = to_fpt(new_point.y()) - to_fpt(site_point.y()); 329 fpt_type a = to_fpt(segment_end.x()) - to_fpt(segment_start.x()); 330 fpt_type b = to_fpt(segment_end.y()) - to_fpt(segment_start.y()); 331 332 if (is_vertical(right_site)) { 333 if (new_point.y() < site_point.y() && !reverse_order) 334 return MORE; 335 else if (new_point.y() > site_point.y() && reverse_order) 336 return LESS; 337 return UNDEFINED; 338 } else { 339 typename ot::Orientation orientation = ot::eval( 340 static_cast<int_x2_type>(segment_end.x()) - 341 static_cast<int_x2_type>(segment_start.x()), 342 static_cast<int_x2_type>(segment_end.y()) - 343 static_cast<int_x2_type>(segment_start.y()), 344 static_cast<int_x2_type>(new_point.x()) - 345 static_cast<int_x2_type>(site_point.x()), 346 static_cast<int_x2_type>(new_point.y()) - 347 static_cast<int_x2_type>(site_point.y())); 348 if (orientation == ot::LEFT) { 349 if (!right_site.is_inverse()) 350 return reverse_order ? LESS : UNDEFINED; 351 return reverse_order ? UNDEFINED : MORE; 352 } 353 } 354 355 fpt_type fast_left_expr = a * (dif_y + dif_x) * (dif_y - dif_x); 356 fpt_type fast_right_expr = (to_fpt(2.0) * b) * dif_x * dif_y; 357 typename ulp_cmp_type::Result expr_cmp = 358 ulp_cmp(fast_left_expr, fast_right_expr, 4); 359 if (expr_cmp != ulp_cmp_type::EQUAL) { 360 if ((expr_cmp == ulp_cmp_type::MORE) ^ reverse_order) 361 return reverse_order ? LESS : MORE; 362 return UNDEFINED; 363 } 364 return UNDEFINED; 365 } 366 367 private: 368 ulp_cmp_type ulp_cmp; 369 to_fpt_converter to_fpt; 370 }; 371 372 template <typename Node> 373 class node_comparison_predicate { 374 public: 375 typedef Node node_type; 376 typedef typename Node::site_type site_type; 377 typedef typename site_type::point_type point_type; 378 typedef typename point_type::coordinate_type coordinate_type; 379 typedef point_comparison_predicate<point_type> point_comparison_type; 380 typedef distance_predicate<site_type> distance_predicate_type; 381 382 // Compares nodes in the balanced binary search tree. Nodes are 383 // compared based on the y coordinates of the arcs intersection points. 384 // Nodes with less y coordinate of the intersection point go first. 385 // Comparison is only called during the new site events processing. 386 // That's why one of the nodes will always lie on the sweepline and may 387 // be represented as a straight horizontal line. operator ()(const node_type & node1,const node_type & node2) const388 bool operator() (const node_type& node1, 389 const node_type& node2) const { 390 // Get x coordinate of the rightmost site from both nodes. 391 const site_type& site1 = get_comparison_site(node1); 392 const site_type& site2 = get_comparison_site(node2); 393 const point_type& point1 = get_comparison_point(site1); 394 const point_type& point2 = get_comparison_point(site2); 395 396 if (point1.x() < point2.x()) { 397 // The second node contains a new site. 398 return distance_predicate_( 399 node1.left_site(), node1.right_site(), point2); 400 } else if (point1.x() > point2.x()) { 401 // The first node contains a new site. 402 return !distance_predicate_( 403 node2.left_site(), node2.right_site(), point1); 404 } else { 405 // This checks were evaluated experimentally. 406 if (site1.sorted_index() == site2.sorted_index()) { 407 // Both nodes are new (inserted during same site event processing). 408 return get_comparison_y(node1) < get_comparison_y(node2); 409 } else if (site1.sorted_index() < site2.sorted_index()) { 410 std::pair<coordinate_type, int> y1 = get_comparison_y(node1, false); 411 std::pair<coordinate_type, int> y2 = get_comparison_y(node2, true); 412 if (y1.first != y2.first) return y1.first < y2.first; 413 return (!site1.is_segment()) ? (y1.second < 0) : false; 414 } else { 415 std::pair<coordinate_type, int> y1 = get_comparison_y(node1, true); 416 std::pair<coordinate_type, int> y2 = get_comparison_y(node2, false); 417 if (y1.first != y2.first) return y1.first < y2.first; 418 return (!site2.is_segment()) ? (y2.second > 0) : true; 419 } 420 } 421 } 422 423 private: 424 // Get the newer site. get_comparison_site(const node_type & node) const425 const site_type& get_comparison_site(const node_type& node) const { 426 if (node.left_site().sorted_index() > node.right_site().sorted_index()) { 427 return node.left_site(); 428 } 429 return node.right_site(); 430 } 431 get_comparison_point(const site_type & site) const432 const point_type& get_comparison_point(const site_type& site) const { 433 return point_comparison_(site.point0(), site.point1()) ? 434 site.point0() : site.point1(); 435 } 436 437 // Get comparison pair: y coordinate and direction of the newer site. get_comparison_y(const node_type & node,bool is_new_node=true) const438 std::pair<coordinate_type, int> get_comparison_y( 439 const node_type& node, bool is_new_node = true) const { 440 if (node.left_site().sorted_index() == 441 node.right_site().sorted_index()) { 442 return std::make_pair(node.left_site().y0(), 0); 443 } 444 if (node.left_site().sorted_index() > node.right_site().sorted_index()) { 445 if (!is_new_node && 446 node.left_site().is_segment() && 447 is_vertical(node.left_site())) { 448 return std::make_pair(node.left_site().y0(), 1); 449 } 450 return std::make_pair(node.left_site().y1(), 1); 451 } 452 return std::make_pair(node.right_site().y0(), -1); 453 } 454 455 point_comparison_type point_comparison_; 456 distance_predicate_type distance_predicate_; 457 }; 458 459 template <typename Site> 460 class circle_existence_predicate { 461 public: 462 typedef typename Site::point_type point_type; 463 typedef Site site_type; 464 ppp(const site_type & site1,const site_type & site2,const site_type & site3) const465 bool ppp(const site_type& site1, 466 const site_type& site2, 467 const site_type& site3) const { 468 return ot::eval(site1.point0(), 469 site2.point0(), 470 site3.point0()) == ot::RIGHT; 471 } 472 pps(const site_type & site1,const site_type & site2,const site_type & site3,int segment_index) const473 bool pps(const site_type& site1, 474 const site_type& site2, 475 const site_type& site3, 476 int segment_index) const { 477 if (segment_index != 2) { 478 typename ot::Orientation orient1 = ot::eval( 479 site1.point0(), site2.point0(), site3.point0()); 480 typename ot::Orientation orient2 = ot::eval( 481 site1.point0(), site2.point0(), site3.point1()); 482 if (segment_index == 1 && site1.x0() >= site2.x0()) { 483 if (orient1 != ot::RIGHT) 484 return false; 485 } else if (segment_index == 3 && site2.x0() >= site1.x0()) { 486 if (orient2 != ot::RIGHT) 487 return false; 488 } else if (orient1 != ot::RIGHT && orient2 != ot::RIGHT) { 489 return false; 490 } 491 } else { 492 return (site3.point0() != site1.point0()) || 493 (site3.point1() != site2.point0()); 494 } 495 return true; 496 } 497 pss(const site_type & site1,const site_type & site2,const site_type & site3,int point_index) const498 bool pss(const site_type& site1, 499 const site_type& site2, 500 const site_type& site3, 501 int point_index) const { 502 if (site2.sorted_index() == site3.sorted_index()) { 503 return false; 504 } 505 if (point_index == 2) { 506 if (!site2.is_inverse() && site3.is_inverse()) 507 return false; 508 if (site2.is_inverse() == site3.is_inverse() && 509 ot::eval(site2.point0(), 510 site1.point0(), 511 site3.point1()) != ot::RIGHT) 512 return false; 513 } 514 return true; 515 } 516 sss(const site_type & site1,const site_type & site2,const site_type & site3) const517 bool sss(const site_type& site1, 518 const site_type& site2, 519 const site_type& site3) const { 520 return (site1.sorted_index() != site2.sorted_index()) && 521 (site2.sorted_index() != site3.sorted_index()); 522 } 523 }; 524 525 template <typename Site, typename Circle> 526 class mp_circle_formation_functor { 527 public: 528 typedef typename Site::point_type point_type; 529 typedef Site site_type; 530 typedef Circle circle_type; 531 typedef robust_sqrt_expr<big_int_type, efpt_type, to_efpt_converter> 532 robust_sqrt_expr_type; 533 ppp(const site_type & site1,const site_type & site2,const site_type & site3,circle_type & circle,bool recompute_c_x=true,bool recompute_c_y=true,bool recompute_lower_x=true)534 void ppp(const site_type& site1, 535 const site_type& site2, 536 const site_type& site3, 537 circle_type& circle, 538 bool recompute_c_x = true, 539 bool recompute_c_y = true, 540 bool recompute_lower_x = true) { 541 big_int_type dif_x[3], dif_y[3], sum_x[2], sum_y[2]; 542 dif_x[0] = static_cast<int_x2_type>(site1.x()) - 543 static_cast<int_x2_type>(site2.x()); 544 dif_x[1] = static_cast<int_x2_type>(site2.x()) - 545 static_cast<int_x2_type>(site3.x()); 546 dif_x[2] = static_cast<int_x2_type>(site1.x()) - 547 static_cast<int_x2_type>(site3.x()); 548 dif_y[0] = static_cast<int_x2_type>(site1.y()) - 549 static_cast<int_x2_type>(site2.y()); 550 dif_y[1] = static_cast<int_x2_type>(site2.y()) - 551 static_cast<int_x2_type>(site3.y()); 552 dif_y[2] = static_cast<int_x2_type>(site1.y()) - 553 static_cast<int_x2_type>(site3.y()); 554 sum_x[0] = static_cast<int_x2_type>(site1.x()) + 555 static_cast<int_x2_type>(site2.x()); 556 sum_x[1] = static_cast<int_x2_type>(site2.x()) + 557 static_cast<int_x2_type>(site3.x()); 558 sum_y[0] = static_cast<int_x2_type>(site1.y()) + 559 static_cast<int_x2_type>(site2.y()); 560 sum_y[1] = static_cast<int_x2_type>(site2.y()) + 561 static_cast<int_x2_type>(site3.y()); 562 fpt_type inv_denom = to_fpt(0.5) / to_fpt(static_cast<big_int_type>( 563 dif_x[0] * dif_y[1] - dif_x[1] * dif_y[0])); 564 big_int_type numer1 = dif_x[0] * sum_x[0] + dif_y[0] * sum_y[0]; 565 big_int_type numer2 = dif_x[1] * sum_x[1] + dif_y[1] * sum_y[1]; 566 567 if (recompute_c_x || recompute_lower_x) { 568 big_int_type c_x = numer1 * dif_y[1] - numer2 * dif_y[0]; 569 circle.x(to_fpt(c_x) * inv_denom); 570 571 if (recompute_lower_x) { 572 // Evaluate radius of the circle. 573 big_int_type sqr_r = (dif_x[0] * dif_x[0] + dif_y[0] * dif_y[0]) * 574 (dif_x[1] * dif_x[1] + dif_y[1] * dif_y[1]) * 575 (dif_x[2] * dif_x[2] + dif_y[2] * dif_y[2]); 576 fpt_type r = get_sqrt(to_fpt(sqr_r)); 577 578 // If c_x >= 0 then lower_x = c_x + r, 579 // else lower_x = (c_x * c_x - r * r) / (c_x - r). 580 // To guarantee epsilon relative error. 581 if (!is_neg(circle.x())) { 582 if (!is_neg(inv_denom)) { 583 circle.lower_x(circle.x() + r * inv_denom); 584 } else { 585 circle.lower_x(circle.x() - r * inv_denom); 586 } 587 } else { 588 big_int_type numer = c_x * c_x - sqr_r; 589 fpt_type lower_x = to_fpt(numer) * inv_denom / (to_fpt(c_x) + r); 590 circle.lower_x(lower_x); 591 } 592 } 593 } 594 595 if (recompute_c_y) { 596 big_int_type c_y = numer2 * dif_x[0] - numer1 * dif_x[1]; 597 circle.y(to_fpt(c_y) * inv_denom); 598 } 599 } 600 601 // Recompute parameters of the circle event using high-precision library. pps(const site_type & site1,const site_type & site2,const site_type & site3,int segment_index,circle_type & c_event,bool recompute_c_x=true,bool recompute_c_y=true,bool recompute_lower_x=true)602 void pps(const site_type& site1, 603 const site_type& site2, 604 const site_type& site3, 605 int segment_index, 606 circle_type& c_event, 607 bool recompute_c_x = true, 608 bool recompute_c_y = true, 609 bool recompute_lower_x = true) { 610 big_int_type cA[4], cB[4]; 611 big_int_type line_a = static_cast<int_x2_type>(site3.y1()) - 612 static_cast<int_x2_type>(site3.y0()); 613 big_int_type line_b = static_cast<int_x2_type>(site3.x0()) - 614 static_cast<int_x2_type>(site3.x1()); 615 big_int_type segm_len = line_a * line_a + line_b * line_b; 616 big_int_type vec_x = static_cast<int_x2_type>(site2.y()) - 617 static_cast<int_x2_type>(site1.y()); 618 big_int_type vec_y = static_cast<int_x2_type>(site1.x()) - 619 static_cast<int_x2_type>(site2.x()); 620 big_int_type sum_x = static_cast<int_x2_type>(site1.x()) + 621 static_cast<int_x2_type>(site2.x()); 622 big_int_type sum_y = static_cast<int_x2_type>(site1.y()) + 623 static_cast<int_x2_type>(site2.y()); 624 big_int_type teta = line_a * vec_x + line_b * vec_y; 625 big_int_type denom = vec_x * line_b - vec_y * line_a; 626 627 big_int_type dif0 = static_cast<int_x2_type>(site3.y1()) - 628 static_cast<int_x2_type>(site1.y()); 629 big_int_type dif1 = static_cast<int_x2_type>(site1.x()) - 630 static_cast<int_x2_type>(site3.x1()); 631 big_int_type A = line_a * dif1 - line_b * dif0; 632 dif0 = static_cast<int_x2_type>(site3.y1()) - 633 static_cast<int_x2_type>(site2.y()); 634 dif1 = static_cast<int_x2_type>(site2.x()) - 635 static_cast<int_x2_type>(site3.x1()); 636 big_int_type B = line_a * dif1 - line_b * dif0; 637 big_int_type sum_AB = A + B; 638 639 if (is_zero(denom)) { 640 big_int_type numer = teta * teta - sum_AB * sum_AB; 641 denom = teta * sum_AB; 642 cA[0] = denom * sum_x * 2 + numer * vec_x; 643 cB[0] = segm_len; 644 cA[1] = denom * sum_AB * 2 + numer * teta; 645 cB[1] = 1; 646 cA[2] = denom * sum_y * 2 + numer * vec_y; 647 fpt_type inv_denom = to_fpt(1.0) / to_fpt(denom); 648 if (recompute_c_x) 649 c_event.x(to_fpt(0.25) * to_fpt(cA[0]) * inv_denom); 650 if (recompute_c_y) 651 c_event.y(to_fpt(0.25) * to_fpt(cA[2]) * inv_denom); 652 if (recompute_lower_x) { 653 c_event.lower_x(to_fpt(0.25) * to_fpt(sqrt_expr_.eval2(cA, cB)) * 654 inv_denom / get_sqrt(to_fpt(segm_len))); 655 } 656 return; 657 } 658 659 big_int_type det = (teta * teta + denom * denom) * A * B * 4; 660 fpt_type inv_denom_sqr = to_fpt(1.0) / to_fpt(denom); 661 inv_denom_sqr *= inv_denom_sqr; 662 663 if (recompute_c_x || recompute_lower_x) { 664 cA[0] = sum_x * denom * denom + teta * sum_AB * vec_x; 665 cB[0] = 1; 666 cA[1] = (segment_index == 2) ? -vec_x : vec_x; 667 cB[1] = det; 668 if (recompute_c_x) { 669 c_event.x(to_fpt(0.5) * to_fpt(sqrt_expr_.eval2(cA, cB)) * 670 inv_denom_sqr); 671 } 672 } 673 674 if (recompute_c_y || recompute_lower_x) { 675 cA[2] = sum_y * denom * denom + teta * sum_AB * vec_y; 676 cB[2] = 1; 677 cA[3] = (segment_index == 2) ? -vec_y : vec_y; 678 cB[3] = det; 679 if (recompute_c_y) { 680 c_event.y(to_fpt(0.5) * to_fpt(sqrt_expr_.eval2(&cA[2], &cB[2])) * 681 inv_denom_sqr); 682 } 683 } 684 685 if (recompute_lower_x) { 686 cB[0] = cB[0] * segm_len; 687 cB[1] = cB[1] * segm_len; 688 cA[2] = sum_AB * (denom * denom + teta * teta); 689 cB[2] = 1; 690 cA[3] = (segment_index == 2) ? -teta : teta; 691 cB[3] = det; 692 c_event.lower_x(to_fpt(0.5) * to_fpt(sqrt_expr_.eval4(cA, cB)) * 693 inv_denom_sqr / get_sqrt(to_fpt(segm_len))); 694 } 695 } 696 697 // Recompute parameters of the circle event using high-precision library. pss(const site_type & site1,const site_type & site2,const site_type & site3,int point_index,circle_type & c_event,bool recompute_c_x=true,bool recompute_c_y=true,bool recompute_lower_x=true)698 void pss(const site_type& site1, 699 const site_type& site2, 700 const site_type& site3, 701 int point_index, 702 circle_type& c_event, 703 bool recompute_c_x = true, 704 bool recompute_c_y = true, 705 bool recompute_lower_x = true) { 706 big_int_type a[2], b[2], c[2], cA[4], cB[4]; 707 const point_type& segm_start1 = site2.point1(); 708 const point_type& segm_end1 = site2.point0(); 709 const point_type& segm_start2 = site3.point0(); 710 const point_type& segm_end2 = site3.point1(); 711 a[0] = static_cast<int_x2_type>(segm_end1.x()) - 712 static_cast<int_x2_type>(segm_start1.x()); 713 b[0] = static_cast<int_x2_type>(segm_end1.y()) - 714 static_cast<int_x2_type>(segm_start1.y()); 715 a[1] = static_cast<int_x2_type>(segm_end2.x()) - 716 static_cast<int_x2_type>(segm_start2.x()); 717 b[1] = static_cast<int_x2_type>(segm_end2.y()) - 718 static_cast<int_x2_type>(segm_start2.y()); 719 big_int_type orientation = a[1] * b[0] - a[0] * b[1]; 720 if (is_zero(orientation)) { 721 fpt_type denom = to_fpt(2.0) * to_fpt( 722 static_cast<big_int_type>(a[0] * a[0] + b[0] * b[0])); 723 c[0] = b[0] * (static_cast<int_x2_type>(segm_start2.x()) - 724 static_cast<int_x2_type>(segm_start1.x())) - 725 a[0] * (static_cast<int_x2_type>(segm_start2.y()) - 726 static_cast<int_x2_type>(segm_start1.y())); 727 big_int_type dx = a[0] * (static_cast<int_x2_type>(site1.y()) - 728 static_cast<int_x2_type>(segm_start1.y())) - 729 b[0] * (static_cast<int_x2_type>(site1.x()) - 730 static_cast<int_x2_type>(segm_start1.x())); 731 big_int_type dy = b[0] * (static_cast<int_x2_type>(site1.x()) - 732 static_cast<int_x2_type>(segm_start2.x())) - 733 a[0] * (static_cast<int_x2_type>(site1.y()) - 734 static_cast<int_x2_type>(segm_start2.y())); 735 cB[0] = dx * dy; 736 cB[1] = 1; 737 738 if (recompute_c_y) { 739 cA[0] = b[0] * ((point_index == 2) ? 2 : -2); 740 cA[1] = a[0] * a[0] * (static_cast<int_x2_type>(segm_start1.y()) + 741 static_cast<int_x2_type>(segm_start2.y())) - 742 a[0] * b[0] * (static_cast<int_x2_type>(segm_start1.x()) + 743 static_cast<int_x2_type>(segm_start2.x()) - 744 static_cast<int_x2_type>(site1.x()) * 2) + 745 b[0] * b[0] * (static_cast<int_x2_type>(site1.y()) * 2); 746 fpt_type c_y = to_fpt(sqrt_expr_.eval2(cA, cB)); 747 c_event.y(c_y / denom); 748 } 749 750 if (recompute_c_x || recompute_lower_x) { 751 cA[0] = a[0] * ((point_index == 2) ? 2 : -2); 752 cA[1] = b[0] * b[0] * (static_cast<int_x2_type>(segm_start1.x()) + 753 static_cast<int_x2_type>(segm_start2.x())) - 754 a[0] * b[0] * (static_cast<int_x2_type>(segm_start1.y()) + 755 static_cast<int_x2_type>(segm_start2.y()) - 756 static_cast<int_x2_type>(site1.y()) * 2) + 757 a[0] * a[0] * (static_cast<int_x2_type>(site1.x()) * 2); 758 759 if (recompute_c_x) { 760 fpt_type c_x = to_fpt(sqrt_expr_.eval2(cA, cB)); 761 c_event.x(c_x / denom); 762 } 763 764 if (recompute_lower_x) { 765 cA[2] = is_neg(c[0]) ? -c[0] : c[0]; 766 cB[2] = a[0] * a[0] + b[0] * b[0]; 767 fpt_type lower_x = to_fpt(sqrt_expr_.eval3(cA, cB)); 768 c_event.lower_x(lower_x / denom); 769 } 770 } 771 return; 772 } 773 c[0] = b[0] * segm_end1.x() - a[0] * segm_end1.y(); 774 c[1] = a[1] * segm_end2.y() - b[1] * segm_end2.x(); 775 big_int_type ix = a[0] * c[1] + a[1] * c[0]; 776 big_int_type iy = b[0] * c[1] + b[1] * c[0]; 777 big_int_type dx = ix - orientation * site1.x(); 778 big_int_type dy = iy - orientation * site1.y(); 779 if (is_zero(dx) && is_zero(dy)) { 780 fpt_type denom = to_fpt(orientation); 781 fpt_type c_x = to_fpt(ix) / denom; 782 fpt_type c_y = to_fpt(iy) / denom; 783 c_event = circle_type(c_x, c_y, c_x); 784 return; 785 } 786 787 big_int_type sign = ((point_index == 2) ? 1 : -1) * 788 (is_neg(orientation) ? 1 : -1); 789 cA[0] = a[1] * -dx + b[1] * -dy; 790 cA[1] = a[0] * -dx + b[0] * -dy; 791 cA[2] = sign; 792 cA[3] = 0; 793 cB[0] = a[0] * a[0] + b[0] * b[0]; 794 cB[1] = a[1] * a[1] + b[1] * b[1]; 795 cB[2] = a[0] * a[1] + b[0] * b[1]; 796 cB[3] = (a[0] * dy - b[0] * dx) * (a[1] * dy - b[1] * dx) * -2; 797 fpt_type temp = to_fpt( 798 sqrt_expr_evaluator_pss4<big_int_type, efpt_type>(cA, cB)); 799 fpt_type denom = temp * to_fpt(orientation); 800 801 if (recompute_c_y) { 802 cA[0] = b[1] * (dx * dx + dy * dy) - iy * (dx * a[1] + dy * b[1]); 803 cA[1] = b[0] * (dx * dx + dy * dy) - iy * (dx * a[0] + dy * b[0]); 804 cA[2] = iy * sign; 805 fpt_type cy = to_fpt( 806 sqrt_expr_evaluator_pss4<big_int_type, efpt_type>(cA, cB)); 807 c_event.y(cy / denom); 808 } 809 810 if (recompute_c_x || recompute_lower_x) { 811 cA[0] = a[1] * (dx * dx + dy * dy) - ix * (dx * a[1] + dy * b[1]); 812 cA[1] = a[0] * (dx * dx + dy * dy) - ix * (dx * a[0] + dy * b[0]); 813 cA[2] = ix * sign; 814 815 if (recompute_c_x) { 816 fpt_type cx = to_fpt( 817 sqrt_expr_evaluator_pss4<big_int_type, efpt_type>(cA, cB)); 818 c_event.x(cx / denom); 819 } 820 821 if (recompute_lower_x) { 822 cA[3] = orientation * (dx * dx + dy * dy) * (is_neg(temp) ? -1 : 1); 823 fpt_type lower_x = to_fpt( 824 sqrt_expr_evaluator_pss4<big_int_type, efpt_type>(cA, cB)); 825 c_event.lower_x(lower_x / denom); 826 } 827 } 828 } 829 830 // Recompute parameters of the circle event using high-precision library. sss(const site_type & site1,const site_type & site2,const site_type & site3,circle_type & c_event,bool recompute_c_x=true,bool recompute_c_y=true,bool recompute_lower_x=true)831 void sss(const site_type& site1, 832 const site_type& site2, 833 const site_type& site3, 834 circle_type& c_event, 835 bool recompute_c_x = true, 836 bool recompute_c_y = true, 837 bool recompute_lower_x = true) { 838 big_int_type a[3], b[3], c[3], cA[4], cB[4]; 839 // cA - corresponds to the cross product. 840 // cB - corresponds to the squared length. 841 a[0] = static_cast<int_x2_type>(site1.x1()) - 842 static_cast<int_x2_type>(site1.x0()); 843 a[1] = static_cast<int_x2_type>(site2.x1()) - 844 static_cast<int_x2_type>(site2.x0()); 845 a[2] = static_cast<int_x2_type>(site3.x1()) - 846 static_cast<int_x2_type>(site3.x0()); 847 848 b[0] = static_cast<int_x2_type>(site1.y1()) - 849 static_cast<int_x2_type>(site1.y0()); 850 b[1] = static_cast<int_x2_type>(site2.y1()) - 851 static_cast<int_x2_type>(site2.y0()); 852 b[2] = static_cast<int_x2_type>(site3.y1()) - 853 static_cast<int_x2_type>(site3.y0()); 854 855 c[0] = static_cast<int_x2_type>(site1.x0()) * 856 static_cast<int_x2_type>(site1.y1()) - 857 static_cast<int_x2_type>(site1.y0()) * 858 static_cast<int_x2_type>(site1.x1()); 859 c[1] = static_cast<int_x2_type>(site2.x0()) * 860 static_cast<int_x2_type>(site2.y1()) - 861 static_cast<int_x2_type>(site2.y0()) * 862 static_cast<int_x2_type>(site2.x1()); 863 c[2] = static_cast<int_x2_type>(site3.x0()) * 864 static_cast<int_x2_type>(site3.y1()) - 865 static_cast<int_x2_type>(site3.y0()) * 866 static_cast<int_x2_type>(site3.x1()); 867 868 for (int i = 0; i < 3; ++i) 869 cB[i] = a[i] * a[i] + b[i] * b[i]; 870 871 for (int i = 0; i < 3; ++i) { 872 int j = (i+1) % 3; 873 int k = (i+2) % 3; 874 cA[i] = a[j] * b[k] - a[k] * b[j]; 875 } 876 fpt_type denom = to_fpt(sqrt_expr_.eval3(cA, cB)); 877 878 if (recompute_c_y) { 879 for (int i = 0; i < 3; ++i) { 880 int j = (i+1) % 3; 881 int k = (i+2) % 3; 882 cA[i] = b[j] * c[k] - b[k] * c[j]; 883 } 884 fpt_type c_y = to_fpt(sqrt_expr_.eval3(cA, cB)); 885 c_event.y(c_y / denom); 886 } 887 888 if (recompute_c_x || recompute_lower_x) { 889 cA[3] = 0; 890 for (int i = 0; i < 3; ++i) { 891 int j = (i+1) % 3; 892 int k = (i+2) % 3; 893 cA[i] = a[j] * c[k] - a[k] * c[j]; 894 if (recompute_lower_x) { 895 cA[3] = cA[3] + cA[i] * b[i]; 896 } 897 } 898 899 if (recompute_c_x) { 900 fpt_type c_x = to_fpt(sqrt_expr_.eval3(cA, cB)); 901 c_event.x(c_x / denom); 902 } 903 904 if (recompute_lower_x) { 905 cB[3] = 1; 906 fpt_type lower_x = to_fpt(sqrt_expr_.eval4(cA, cB)); 907 c_event.lower_x(lower_x / denom); 908 } 909 } 910 } 911 912 private: 913 // Evaluates A[3] + A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) + 914 // A[2] * sqrt(B[3] * (sqrt(B[0] * B[1]) + B[2])). 915 template <typename _int, typename _fpt> sqrt_expr_evaluator_pss4(_int * A,_int * B)916 _fpt sqrt_expr_evaluator_pss4(_int *A, _int *B) { 917 _int cA[4], cB[4]; 918 if (is_zero(A[3])) { 919 _fpt lh = sqrt_expr_.eval2(A, B); 920 cA[0] = 1; 921 cB[0] = B[0] * B[1]; 922 cA[1] = B[2]; 923 cB[1] = 1; 924 _fpt rh = sqrt_expr_.eval1(A+2, B+3) * 925 get_sqrt(sqrt_expr_.eval2(cA, cB)); 926 if ((!is_neg(lh) && !is_neg(rh)) || (!is_pos(lh) && !is_pos(rh))) 927 return lh + rh; 928 cA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1] - 929 A[2] * A[2] * B[3] * B[2]; 930 cB[0] = 1; 931 cA[1] = A[0] * A[1] * 2 - A[2] * A[2] * B[3]; 932 cB[1] = B[0] * B[1]; 933 _fpt numer = sqrt_expr_.eval2(cA, cB); 934 return numer / (lh - rh); 935 } 936 cA[0] = 1; 937 cB[0] = B[0] * B[1]; 938 cA[1] = B[2]; 939 cB[1] = 1; 940 _fpt rh = sqrt_expr_.eval1(A+2, B+3) * get_sqrt(sqrt_expr_.eval2(cA, cB)); 941 cA[0] = A[0]; 942 cB[0] = B[0]; 943 cA[1] = A[1]; 944 cB[1] = B[1]; 945 cA[2] = A[3]; 946 cB[2] = 1; 947 _fpt lh = sqrt_expr_.eval3(cA, cB); 948 if ((!is_neg(lh) && !is_neg(rh)) || (!is_pos(lh) && !is_pos(rh))) 949 return lh + rh; 950 cA[0] = A[3] * A[0] * 2; 951 cA[1] = A[3] * A[1] * 2; 952 cA[2] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1] + 953 A[3] * A[3] - A[2] * A[2] * B[2] * B[3]; 954 cA[3] = A[0] * A[1] * 2 - A[2] * A[2] * B[3]; 955 cB[3] = B[0] * B[1]; 956 _fpt numer = sqrt_expr_evaluator_pss3<_int, _fpt>(cA, cB); 957 return numer / (lh - rh); 958 } 959 960 template <typename _int, typename _fpt> 961 // Evaluates A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) + 962 // A[2] + A[3] * sqrt(B[0] * B[1]). 963 // B[3] = B[0] * B[1]. sqrt_expr_evaluator_pss3(_int * A,_int * B)964 _fpt sqrt_expr_evaluator_pss3(_int *A, _int *B) { 965 _int cA[2], cB[2]; 966 _fpt lh = sqrt_expr_.eval2(A, B); 967 _fpt rh = sqrt_expr_.eval2(A+2, B+2); 968 if ((!is_neg(lh) && !is_neg(rh)) || (!is_pos(lh) && !is_pos(rh))) 969 return lh + rh; 970 cA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1] - 971 A[2] * A[2] - A[3] * A[3] * B[0] * B[1]; 972 cB[0] = 1; 973 cA[1] = (A[0] * A[1] - A[2] * A[3]) * 2; 974 cB[1] = B[3]; 975 _fpt numer = sqrt_expr_.eval2(cA, cB); 976 return numer / (lh - rh); 977 } 978 979 robust_sqrt_expr_type sqrt_expr_; 980 to_fpt_converter to_fpt; 981 }; 982 983 template <typename Site, typename Circle> 984 class lazy_circle_formation_functor { 985 public: 986 typedef robust_fpt<fpt_type> robust_fpt_type; 987 typedef robust_dif<robust_fpt_type> robust_dif_type; 988 typedef typename Site::point_type point_type; 989 typedef Site site_type; 990 typedef Circle circle_type; 991 typedef mp_circle_formation_functor<site_type, circle_type> 992 exact_circle_formation_functor_type; 993 ppp(const site_type & site1,const site_type & site2,const site_type & site3,circle_type & c_event)994 void ppp(const site_type& site1, 995 const site_type& site2, 996 const site_type& site3, 997 circle_type& c_event) { 998 fpt_type dif_x1 = to_fpt(site1.x()) - to_fpt(site2.x()); 999 fpt_type dif_x2 = to_fpt(site2.x()) - to_fpt(site3.x()); 1000 fpt_type dif_y1 = to_fpt(site1.y()) - to_fpt(site2.y()); 1001 fpt_type dif_y2 = to_fpt(site2.y()) - to_fpt(site3.y()); 1002 fpt_type orientation = robust_cross_product( 1003 static_cast<int_x2_type>(site1.x()) - 1004 static_cast<int_x2_type>(site2.x()), 1005 static_cast<int_x2_type>(site2.x()) - 1006 static_cast<int_x2_type>(site3.x()), 1007 static_cast<int_x2_type>(site1.y()) - 1008 static_cast<int_x2_type>(site2.y()), 1009 static_cast<int_x2_type>(site2.y()) - 1010 static_cast<int_x2_type>(site3.y())); 1011 robust_fpt_type inv_orientation(to_fpt(0.5) / orientation, to_fpt(2.0)); 1012 fpt_type sum_x1 = to_fpt(site1.x()) + to_fpt(site2.x()); 1013 fpt_type sum_x2 = to_fpt(site2.x()) + to_fpt(site3.x()); 1014 fpt_type sum_y1 = to_fpt(site1.y()) + to_fpt(site2.y()); 1015 fpt_type sum_y2 = to_fpt(site2.y()) + to_fpt(site3.y()); 1016 fpt_type dif_x3 = to_fpt(site1.x()) - to_fpt(site3.x()); 1017 fpt_type dif_y3 = to_fpt(site1.y()) - to_fpt(site3.y()); 1018 robust_dif_type c_x, c_y; 1019 c_x += robust_fpt_type(dif_x1 * sum_x1 * dif_y2, to_fpt(2.0)); 1020 c_x += robust_fpt_type(dif_y1 * sum_y1 * dif_y2, to_fpt(2.0)); 1021 c_x -= robust_fpt_type(dif_x2 * sum_x2 * dif_y1, to_fpt(2.0)); 1022 c_x -= robust_fpt_type(dif_y2 * sum_y2 * dif_y1, to_fpt(2.0)); 1023 c_y += robust_fpt_type(dif_x2 * sum_x2 * dif_x1, to_fpt(2.0)); 1024 c_y += robust_fpt_type(dif_y2 * sum_y2 * dif_x1, to_fpt(2.0)); 1025 c_y -= robust_fpt_type(dif_x1 * sum_x1 * dif_x2, to_fpt(2.0)); 1026 c_y -= robust_fpt_type(dif_y1 * sum_y1 * dif_x2, to_fpt(2.0)); 1027 robust_dif_type lower_x(c_x); 1028 lower_x -= robust_fpt_type(get_sqrt( 1029 (dif_x1 * dif_x1 + dif_y1 * dif_y1) * 1030 (dif_x2 * dif_x2 + dif_y2 * dif_y2) * 1031 (dif_x3 * dif_x3 + dif_y3 * dif_y3)), to_fpt(5.0)); 1032 c_event = circle_type( 1033 c_x.dif().fpv() * inv_orientation.fpv(), 1034 c_y.dif().fpv() * inv_orientation.fpv(), 1035 lower_x.dif().fpv() * inv_orientation.fpv()); 1036 bool recompute_c_x = c_x.dif().ulp() > ULPS; 1037 bool recompute_c_y = c_y.dif().ulp() > ULPS; 1038 bool recompute_lower_x = lower_x.dif().ulp() > ULPS; 1039 if (recompute_c_x || recompute_c_y || recompute_lower_x) { 1040 exact_circle_formation_functor_.ppp( 1041 site1, site2, site3, c_event, 1042 recompute_c_x, recompute_c_y, recompute_lower_x); 1043 } 1044 } 1045 pps(const site_type & site1,const site_type & site2,const site_type & site3,int segment_index,circle_type & c_event)1046 void pps(const site_type& site1, 1047 const site_type& site2, 1048 const site_type& site3, 1049 int segment_index, 1050 circle_type& c_event) { 1051 fpt_type line_a = to_fpt(site3.y1()) - to_fpt(site3.y0()); 1052 fpt_type line_b = to_fpt(site3.x0()) - to_fpt(site3.x1()); 1053 fpt_type vec_x = to_fpt(site2.y()) - to_fpt(site1.y()); 1054 fpt_type vec_y = to_fpt(site1.x()) - to_fpt(site2.x()); 1055 robust_fpt_type teta(robust_cross_product( 1056 static_cast<int_x2_type>(site3.y1()) - 1057 static_cast<int_x2_type>(site3.y0()), 1058 static_cast<int_x2_type>(site3.x0()) - 1059 static_cast<int_x2_type>(site3.x1()), 1060 static_cast<int_x2_type>(site2.x()) - 1061 static_cast<int_x2_type>(site1.x()), 1062 static_cast<int_x2_type>(site2.y()) - 1063 static_cast<int_x2_type>(site1.y())), to_fpt(1.0)); 1064 robust_fpt_type A(robust_cross_product( 1065 static_cast<int_x2_type>(site3.y0()) - 1066 static_cast<int_x2_type>(site3.y1()), 1067 static_cast<int_x2_type>(site3.x0()) - 1068 static_cast<int_x2_type>(site3.x1()), 1069 static_cast<int_x2_type>(site3.y1()) - 1070 static_cast<int_x2_type>(site1.y()), 1071 static_cast<int_x2_type>(site3.x1()) - 1072 static_cast<int_x2_type>(site1.x())), to_fpt(1.0)); 1073 robust_fpt_type B(robust_cross_product( 1074 static_cast<int_x2_type>(site3.y0()) - 1075 static_cast<int_x2_type>(site3.y1()), 1076 static_cast<int_x2_type>(site3.x0()) - 1077 static_cast<int_x2_type>(site3.x1()), 1078 static_cast<int_x2_type>(site3.y1()) - 1079 static_cast<int_x2_type>(site2.y()), 1080 static_cast<int_x2_type>(site3.x1()) - 1081 static_cast<int_x2_type>(site2.x())), to_fpt(1.0)); 1082 robust_fpt_type denom(robust_cross_product( 1083 static_cast<int_x2_type>(site1.y()) - 1084 static_cast<int_x2_type>(site2.y()), 1085 static_cast<int_x2_type>(site1.x()) - 1086 static_cast<int_x2_type>(site2.x()), 1087 static_cast<int_x2_type>(site3.y1()) - 1088 static_cast<int_x2_type>(site3.y0()), 1089 static_cast<int_x2_type>(site3.x1()) - 1090 static_cast<int_x2_type>(site3.x0())), to_fpt(1.0)); 1091 robust_fpt_type inv_segm_len(to_fpt(1.0) / 1092 get_sqrt(line_a * line_a + line_b * line_b), to_fpt(3.0)); 1093 robust_dif_type t; 1094 if (ot::eval(denom) == ot::COLLINEAR) { 1095 t += teta / (robust_fpt_type(to_fpt(8.0)) * A); 1096 t -= A / (robust_fpt_type(to_fpt(2.0)) * teta); 1097 } else { 1098 robust_fpt_type det = ((teta * teta + denom * denom) * A * B).sqrt(); 1099 if (segment_index == 2) { 1100 t -= det / (denom * denom); 1101 } else { 1102 t += det / (denom * denom); 1103 } 1104 t += teta * (A + B) / (robust_fpt_type(to_fpt(2.0)) * denom * denom); 1105 } 1106 robust_dif_type c_x, c_y; 1107 c_x += robust_fpt_type(to_fpt(0.5) * 1108 (to_fpt(site1.x()) + to_fpt(site2.x()))); 1109 c_x += robust_fpt_type(vec_x) * t; 1110 c_y += robust_fpt_type(to_fpt(0.5) * 1111 (to_fpt(site1.y()) + to_fpt(site2.y()))); 1112 c_y += robust_fpt_type(vec_y) * t; 1113 robust_dif_type r, lower_x(c_x); 1114 r -= robust_fpt_type(line_a) * robust_fpt_type(site3.x0()); 1115 r -= robust_fpt_type(line_b) * robust_fpt_type(site3.y0()); 1116 r += robust_fpt_type(line_a) * c_x; 1117 r += robust_fpt_type(line_b) * c_y; 1118 if (r.pos().fpv() < r.neg().fpv()) 1119 r = -r; 1120 lower_x += r * inv_segm_len; 1121 c_event = circle_type( 1122 c_x.dif().fpv(), c_y.dif().fpv(), lower_x.dif().fpv()); 1123 bool recompute_c_x = c_x.dif().ulp() > ULPS; 1124 bool recompute_c_y = c_y.dif().ulp() > ULPS; 1125 bool recompute_lower_x = lower_x.dif().ulp() > ULPS; 1126 if (recompute_c_x || recompute_c_y || recompute_lower_x) { 1127 exact_circle_formation_functor_.pps( 1128 site1, site2, site3, segment_index, c_event, 1129 recompute_c_x, recompute_c_y, recompute_lower_x); 1130 } 1131 } 1132 pss(const site_type & site1,const site_type & site2,const site_type & site3,int point_index,circle_type & c_event)1133 void pss(const site_type& site1, 1134 const site_type& site2, 1135 const site_type& site3, 1136 int point_index, 1137 circle_type& c_event) { 1138 const point_type& segm_start1 = site2.point1(); 1139 const point_type& segm_end1 = site2.point0(); 1140 const point_type& segm_start2 = site3.point0(); 1141 const point_type& segm_end2 = site3.point1(); 1142 fpt_type a1 = to_fpt(segm_end1.x()) - to_fpt(segm_start1.x()); 1143 fpt_type b1 = to_fpt(segm_end1.y()) - to_fpt(segm_start1.y()); 1144 fpt_type a2 = to_fpt(segm_end2.x()) - to_fpt(segm_start2.x()); 1145 fpt_type b2 = to_fpt(segm_end2.y()) - to_fpt(segm_start2.y()); 1146 bool recompute_c_x, recompute_c_y, recompute_lower_x; 1147 robust_fpt_type orientation(robust_cross_product( 1148 static_cast<int_x2_type>(segm_end1.y()) - 1149 static_cast<int_x2_type>(segm_start1.y()), 1150 static_cast<int_x2_type>(segm_end1.x()) - 1151 static_cast<int_x2_type>(segm_start1.x()), 1152 static_cast<int_x2_type>(segm_end2.y()) - 1153 static_cast<int_x2_type>(segm_start2.y()), 1154 static_cast<int_x2_type>(segm_end2.x()) - 1155 static_cast<int_x2_type>(segm_start2.x())), to_fpt(1.0)); 1156 if (ot::eval(orientation) == ot::COLLINEAR) { 1157 robust_fpt_type a(a1 * a1 + b1 * b1, to_fpt(2.0)); 1158 robust_fpt_type c(robust_cross_product( 1159 static_cast<int_x2_type>(segm_end1.y()) - 1160 static_cast<int_x2_type>(segm_start1.y()), 1161 static_cast<int_x2_type>(segm_end1.x()) - 1162 static_cast<int_x2_type>(segm_start1.x()), 1163 static_cast<int_x2_type>(segm_start2.y()) - 1164 static_cast<int_x2_type>(segm_start1.y()), 1165 static_cast<int_x2_type>(segm_start2.x()) - 1166 static_cast<int_x2_type>(segm_start1.x())), to_fpt(1.0)); 1167 robust_fpt_type det( 1168 robust_cross_product( 1169 static_cast<int_x2_type>(segm_end1.x()) - 1170 static_cast<int_x2_type>(segm_start1.x()), 1171 static_cast<int_x2_type>(segm_end1.y()) - 1172 static_cast<int_x2_type>(segm_start1.y()), 1173 static_cast<int_x2_type>(site1.x()) - 1174 static_cast<int_x2_type>(segm_start1.x()), 1175 static_cast<int_x2_type>(site1.y()) - 1176 static_cast<int_x2_type>(segm_start1.y())) * 1177 robust_cross_product( 1178 static_cast<int_x2_type>(segm_end1.y()) - 1179 static_cast<int_x2_type>(segm_start1.y()), 1180 static_cast<int_x2_type>(segm_end1.x()) - 1181 static_cast<int_x2_type>(segm_start1.x()), 1182 static_cast<int_x2_type>(site1.y()) - 1183 static_cast<int_x2_type>(segm_start2.y()), 1184 static_cast<int_x2_type>(site1.x()) - 1185 static_cast<int_x2_type>(segm_start2.x())), 1186 to_fpt(3.0)); 1187 robust_dif_type t; 1188 t -= robust_fpt_type(a1) * robust_fpt_type(( 1189 to_fpt(segm_start1.x()) + to_fpt(segm_start2.x())) * to_fpt(0.5) - 1190 to_fpt(site1.x())); 1191 t -= robust_fpt_type(b1) * robust_fpt_type(( 1192 to_fpt(segm_start1.y()) + to_fpt(segm_start2.y())) * to_fpt(0.5) - 1193 to_fpt(site1.y())); 1194 if (point_index == 2) { 1195 t += det.sqrt(); 1196 } else { 1197 t -= det.sqrt(); 1198 } 1199 t /= a; 1200 robust_dif_type c_x, c_y; 1201 c_x += robust_fpt_type(to_fpt(0.5) * ( 1202 to_fpt(segm_start1.x()) + to_fpt(segm_start2.x()))); 1203 c_x += robust_fpt_type(a1) * t; 1204 c_y += robust_fpt_type(to_fpt(0.5) * ( 1205 to_fpt(segm_start1.y()) + to_fpt(segm_start2.y()))); 1206 c_y += robust_fpt_type(b1) * t; 1207 robust_dif_type lower_x(c_x); 1208 if (is_neg(c)) { 1209 lower_x -= robust_fpt_type(to_fpt(0.5)) * c / a.sqrt(); 1210 } else { 1211 lower_x += robust_fpt_type(to_fpt(0.5)) * c / a.sqrt(); 1212 } 1213 recompute_c_x = c_x.dif().ulp() > ULPS; 1214 recompute_c_y = c_y.dif().ulp() > ULPS; 1215 recompute_lower_x = lower_x.dif().ulp() > ULPS; 1216 c_event = 1217 circle_type(c_x.dif().fpv(), c_y.dif().fpv(), lower_x.dif().fpv()); 1218 } else { 1219 robust_fpt_type sqr_sum1(get_sqrt(a1 * a1 + b1 * b1), to_fpt(2.0)); 1220 robust_fpt_type sqr_sum2(get_sqrt(a2 * a2 + b2 * b2), to_fpt(2.0)); 1221 robust_fpt_type a(robust_cross_product( 1222 static_cast<int_x2_type>(segm_end1.x()) - 1223 static_cast<int_x2_type>(segm_start1.x()), 1224 static_cast<int_x2_type>(segm_end1.y()) - 1225 static_cast<int_x2_type>(segm_start1.y()), 1226 static_cast<int_x2_type>(segm_start2.y()) - 1227 static_cast<int_x2_type>(segm_end2.y()), 1228 static_cast<int_x2_type>(segm_end2.x()) - 1229 static_cast<int_x2_type>(segm_start2.x())), to_fpt(1.0)); 1230 if (!is_neg(a)) { 1231 a += sqr_sum1 * sqr_sum2; 1232 } else { 1233 a = (orientation * orientation) / (sqr_sum1 * sqr_sum2 - a); 1234 } 1235 robust_fpt_type or1(robust_cross_product( 1236 static_cast<int_x2_type>(segm_end1.y()) - 1237 static_cast<int_x2_type>(segm_start1.y()), 1238 static_cast<int_x2_type>(segm_end1.x()) - 1239 static_cast<int_x2_type>(segm_start1.x()), 1240 static_cast<int_x2_type>(segm_end1.y()) - 1241 static_cast<int_x2_type>(site1.y()), 1242 static_cast<int_x2_type>(segm_end1.x()) - 1243 static_cast<int_x2_type>(site1.x())), to_fpt(1.0)); 1244 robust_fpt_type or2(robust_cross_product( 1245 static_cast<int_x2_type>(segm_end2.x()) - 1246 static_cast<int_x2_type>(segm_start2.x()), 1247 static_cast<int_x2_type>(segm_end2.y()) - 1248 static_cast<int_x2_type>(segm_start2.y()), 1249 static_cast<int_x2_type>(segm_end2.x()) - 1250 static_cast<int_x2_type>(site1.x()), 1251 static_cast<int_x2_type>(segm_end2.y()) - 1252 static_cast<int_x2_type>(site1.y())), to_fpt(1.0)); 1253 robust_fpt_type det = robust_fpt_type(to_fpt(2.0)) * a * or1 * or2; 1254 robust_fpt_type c1(robust_cross_product( 1255 static_cast<int_x2_type>(segm_end1.y()) - 1256 static_cast<int_x2_type>(segm_start1.y()), 1257 static_cast<int_x2_type>(segm_end1.x()) - 1258 static_cast<int_x2_type>(segm_start1.x()), 1259 static_cast<int_x2_type>(segm_end1.y()), 1260 static_cast<int_x2_type>(segm_end1.x())), to_fpt(1.0)); 1261 robust_fpt_type c2(robust_cross_product( 1262 static_cast<int_x2_type>(segm_end2.x()) - 1263 static_cast<int_x2_type>(segm_start2.x()), 1264 static_cast<int_x2_type>(segm_end2.y()) - 1265 static_cast<int_x2_type>(segm_start2.y()), 1266 static_cast<int_x2_type>(segm_end2.x()), 1267 static_cast<int_x2_type>(segm_end2.y())), to_fpt(1.0)); 1268 robust_fpt_type inv_orientation = 1269 robust_fpt_type(to_fpt(1.0)) / orientation; 1270 robust_dif_type t, b, ix, iy; 1271 ix += robust_fpt_type(a2) * c1 * inv_orientation; 1272 ix += robust_fpt_type(a1) * c2 * inv_orientation; 1273 iy += robust_fpt_type(b1) * c2 * inv_orientation; 1274 iy += robust_fpt_type(b2) * c1 * inv_orientation; 1275 1276 b += ix * (robust_fpt_type(a1) * sqr_sum2); 1277 b += ix * (robust_fpt_type(a2) * sqr_sum1); 1278 b += iy * (robust_fpt_type(b1) * sqr_sum2); 1279 b += iy * (robust_fpt_type(b2) * sqr_sum1); 1280 b -= sqr_sum1 * robust_fpt_type(robust_cross_product( 1281 static_cast<int_x2_type>(segm_end2.x()) - 1282 static_cast<int_x2_type>(segm_start2.x()), 1283 static_cast<int_x2_type>(segm_end2.y()) - 1284 static_cast<int_x2_type>(segm_start2.y()), 1285 static_cast<int_x2_type>(-site1.y()), 1286 static_cast<int_x2_type>(site1.x())), to_fpt(1.0)); 1287 b -= sqr_sum2 * robust_fpt_type(robust_cross_product( 1288 static_cast<int_x2_type>(segm_end1.x()) - 1289 static_cast<int_x2_type>(segm_start1.x()), 1290 static_cast<int_x2_type>(segm_end1.y()) - 1291 static_cast<int_x2_type>(segm_start1.y()), 1292 static_cast<int_x2_type>(-site1.y()), 1293 static_cast<int_x2_type>(site1.x())), to_fpt(1.0)); 1294 t -= b; 1295 if (point_index == 2) { 1296 t += det.sqrt(); 1297 } else { 1298 t -= det.sqrt(); 1299 } 1300 t /= (a * a); 1301 robust_dif_type c_x(ix), c_y(iy); 1302 c_x += t * (robust_fpt_type(a1) * sqr_sum2); 1303 c_x += t * (robust_fpt_type(a2) * sqr_sum1); 1304 c_y += t * (robust_fpt_type(b1) * sqr_sum2); 1305 c_y += t * (robust_fpt_type(b2) * sqr_sum1); 1306 if (t.pos().fpv() < t.neg().fpv()) { 1307 t = -t; 1308 } 1309 robust_dif_type lower_x(c_x); 1310 if (is_neg(orientation)) { 1311 lower_x -= t * orientation; 1312 } else { 1313 lower_x += t * orientation; 1314 } 1315 recompute_c_x = c_x.dif().ulp() > ULPS; 1316 recompute_c_y = c_y.dif().ulp() > ULPS; 1317 recompute_lower_x = lower_x.dif().ulp() > ULPS; 1318 c_event = circle_type( 1319 c_x.dif().fpv(), c_y.dif().fpv(), lower_x.dif().fpv()); 1320 } 1321 if (recompute_c_x || recompute_c_y || recompute_lower_x) { 1322 exact_circle_formation_functor_.pss( 1323 site1, site2, site3, point_index, c_event, 1324 recompute_c_x, recompute_c_y, recompute_lower_x); 1325 } 1326 } 1327 sss(const site_type & site1,const site_type & site2,const site_type & site3,circle_type & c_event)1328 void sss(const site_type& site1, 1329 const site_type& site2, 1330 const site_type& site3, 1331 circle_type& c_event) { 1332 robust_fpt_type a1(to_fpt(site1.x1()) - to_fpt(site1.x0())); 1333 robust_fpt_type b1(to_fpt(site1.y1()) - to_fpt(site1.y0())); 1334 robust_fpt_type c1(robust_cross_product( 1335 site1.x0(), site1.y0(), 1336 site1.x1(), site1.y1()), to_fpt(1.0)); 1337 1338 robust_fpt_type a2(to_fpt(site2.x1()) - to_fpt(site2.x0())); 1339 robust_fpt_type b2(to_fpt(site2.y1()) - to_fpt(site2.y0())); 1340 robust_fpt_type c2(robust_cross_product( 1341 site2.x0(), site2.y0(), 1342 site2.x1(), site2.y1()), to_fpt(1.0)); 1343 1344 robust_fpt_type a3(to_fpt(site3.x1()) - to_fpt(site3.x0())); 1345 robust_fpt_type b3(to_fpt(site3.y1()) - to_fpt(site3.y0())); 1346 robust_fpt_type c3(robust_cross_product( 1347 site3.x0(), site3.y0(), 1348 site3.x1(), site3.y1()), to_fpt(1.0)); 1349 1350 robust_fpt_type len1 = (a1 * a1 + b1 * b1).sqrt(); 1351 robust_fpt_type len2 = (a2 * a2 + b2 * b2).sqrt(); 1352 robust_fpt_type len3 = (a3 * a3 + b3 * b3).sqrt(); 1353 robust_fpt_type cross_12(robust_cross_product( 1354 static_cast<int_x2_type>(site1.x1()) - 1355 static_cast<int_x2_type>(site1.x0()), 1356 static_cast<int_x2_type>(site1.y1()) - 1357 static_cast<int_x2_type>(site1.y0()), 1358 static_cast<int_x2_type>(site2.x1()) - 1359 static_cast<int_x2_type>(site2.x0()), 1360 static_cast<int_x2_type>(site2.y1()) - 1361 static_cast<int_x2_type>(site2.y0())), to_fpt(1.0)); 1362 robust_fpt_type cross_23(robust_cross_product( 1363 static_cast<int_x2_type>(site2.x1()) - 1364 static_cast<int_x2_type>(site2.x0()), 1365 static_cast<int_x2_type>(site2.y1()) - 1366 static_cast<int_x2_type>(site2.y0()), 1367 static_cast<int_x2_type>(site3.x1()) - 1368 static_cast<int_x2_type>(site3.x0()), 1369 static_cast<int_x2_type>(site3.y1()) - 1370 static_cast<int_x2_type>(site3.y0())), to_fpt(1.0)); 1371 robust_fpt_type cross_31(robust_cross_product( 1372 static_cast<int_x2_type>(site3.x1()) - 1373 static_cast<int_x2_type>(site3.x0()), 1374 static_cast<int_x2_type>(site3.y1()) - 1375 static_cast<int_x2_type>(site3.y0()), 1376 static_cast<int_x2_type>(site1.x1()) - 1377 static_cast<int_x2_type>(site1.x0()), 1378 static_cast<int_x2_type>(site1.y1()) - 1379 static_cast<int_x2_type>(site1.y0())), to_fpt(1.0)); 1380 1381 // denom = cross_12 * len3 + cross_23 * len1 + cross_31 * len2. 1382 robust_dif_type denom; 1383 denom += cross_12 * len3; 1384 denom += cross_23 * len1; 1385 denom += cross_31 * len2; 1386 1387 // denom * r = (b2 * c_x - a2 * c_y - c2 * denom) / len2. 1388 robust_dif_type r; 1389 r -= cross_12 * c3; 1390 r -= cross_23 * c1; 1391 r -= cross_31 * c2; 1392 1393 robust_dif_type c_x; 1394 c_x += a1 * c2 * len3; 1395 c_x -= a2 * c1 * len3; 1396 c_x += a2 * c3 * len1; 1397 c_x -= a3 * c2 * len1; 1398 c_x += a3 * c1 * len2; 1399 c_x -= a1 * c3 * len2; 1400 1401 robust_dif_type c_y; 1402 c_y += b1 * c2 * len3; 1403 c_y -= b2 * c1 * len3; 1404 c_y += b2 * c3 * len1; 1405 c_y -= b3 * c2 * len1; 1406 c_y += b3 * c1 * len2; 1407 c_y -= b1 * c3 * len2; 1408 1409 robust_dif_type lower_x = c_x + r; 1410 1411 robust_fpt_type denom_dif = denom.dif(); 1412 robust_fpt_type c_x_dif = c_x.dif() / denom_dif; 1413 robust_fpt_type c_y_dif = c_y.dif() / denom_dif; 1414 robust_fpt_type lower_x_dif = lower_x.dif() / denom_dif; 1415 1416 bool recompute_c_x = c_x_dif.ulp() > ULPS; 1417 bool recompute_c_y = c_y_dif.ulp() > ULPS; 1418 bool recompute_lower_x = lower_x_dif.ulp() > ULPS; 1419 c_event = circle_type(c_x_dif.fpv(), c_y_dif.fpv(), lower_x_dif.fpv()); 1420 if (recompute_c_x || recompute_c_y || recompute_lower_x) { 1421 exact_circle_formation_functor_.sss( 1422 site1, site2, site3, c_event, 1423 recompute_c_x, recompute_c_y, recompute_lower_x); 1424 } 1425 } 1426 1427 private: 1428 exact_circle_formation_functor_type exact_circle_formation_functor_; 1429 to_fpt_converter to_fpt; 1430 }; 1431 1432 template <typename Site, 1433 typename Circle, 1434 typename CEP = circle_existence_predicate<Site>, 1435 typename CFF = lazy_circle_formation_functor<Site, Circle> > 1436 class circle_formation_predicate { 1437 public: 1438 typedef Site site_type; 1439 typedef Circle circle_type; 1440 typedef CEP circle_existence_predicate_type; 1441 typedef CFF circle_formation_functor_type; 1442 1443 // Create a circle event from the given three sites. 1444 // Returns true if the circle event exists, else false. 1445 // If exists circle event is saved into the c_event variable. operator ()(const site_type & site1,const site_type & site2,const site_type & site3,circle_type & circle)1446 bool operator()(const site_type& site1, const site_type& site2, 1447 const site_type& site3, circle_type& circle) { 1448 if (!site1.is_segment()) { 1449 if (!site2.is_segment()) { 1450 if (!site3.is_segment()) { 1451 // (point, point, point) sites. 1452 if (!circle_existence_predicate_.ppp(site1, site2, site3)) 1453 return false; 1454 circle_formation_functor_.ppp(site1, site2, site3, circle); 1455 } else { 1456 // (point, point, segment) sites. 1457 if (!circle_existence_predicate_.pps(site1, site2, site3, 3)) 1458 return false; 1459 circle_formation_functor_.pps(site1, site2, site3, 3, circle); 1460 } 1461 } else { 1462 if (!site3.is_segment()) { 1463 // (point, segment, point) sites. 1464 if (!circle_existence_predicate_.pps(site1, site3, site2, 2)) 1465 return false; 1466 circle_formation_functor_.pps(site1, site3, site2, 2, circle); 1467 } else { 1468 // (point, segment, segment) sites. 1469 if (!circle_existence_predicate_.pss(site1, site2, site3, 1)) 1470 return false; 1471 circle_formation_functor_.pss(site1, site2, site3, 1, circle); 1472 } 1473 } 1474 } else { 1475 if (!site2.is_segment()) { 1476 if (!site3.is_segment()) { 1477 // (segment, point, point) sites. 1478 if (!circle_existence_predicate_.pps(site2, site3, site1, 1)) 1479 return false; 1480 circle_formation_functor_.pps(site2, site3, site1, 1, circle); 1481 } else { 1482 // (segment, point, segment) sites. 1483 if (!circle_existence_predicate_.pss(site2, site1, site3, 2)) 1484 return false; 1485 circle_formation_functor_.pss(site2, site1, site3, 2, circle); 1486 } 1487 } else { 1488 if (!site3.is_segment()) { 1489 // (segment, segment, point) sites. 1490 if (!circle_existence_predicate_.pss(site3, site1, site2, 3)) 1491 return false; 1492 circle_formation_functor_.pss(site3, site1, site2, 3, circle); 1493 } else { 1494 // (segment, segment, segment) sites. 1495 if (!circle_existence_predicate_.sss(site1, site2, site3)) 1496 return false; 1497 circle_formation_functor_.sss(site1, site2, site3, circle); 1498 } 1499 } 1500 } 1501 if (lies_outside_vertical_segment(circle, site1) || 1502 lies_outside_vertical_segment(circle, site2) || 1503 lies_outside_vertical_segment(circle, site3)) { 1504 return false; 1505 } 1506 return true; 1507 } 1508 1509 private: lies_outside_vertical_segment(const circle_type & c,const site_type & s)1510 bool lies_outside_vertical_segment( 1511 const circle_type& c, const site_type& s) { 1512 if (!s.is_segment() || !is_vertical(s)) { 1513 return false; 1514 } 1515 fpt_type y0 = to_fpt(s.is_inverse() ? s.y1() : s.y0()); 1516 fpt_type y1 = to_fpt(s.is_inverse() ? s.y0() : s.y1()); 1517 return ulp_cmp(c.y(), y0, ULPS) == ulp_cmp_type::LESS || 1518 ulp_cmp(c.y(), y1, ULPS) == ulp_cmp_type::MORE; 1519 } 1520 1521 private: 1522 to_fpt_converter to_fpt; 1523 ulp_cmp_type ulp_cmp; 1524 circle_existence_predicate_type circle_existence_predicate_; 1525 circle_formation_functor_type circle_formation_functor_; 1526 }; 1527 }; 1528 } // detail 1529 } // polygon 1530 } // boost 1531 1532 #endif // BOOST_POLYGON_DETAIL_VORONOI_PREDICATES 1533