1 // Copyright (c) 2006-2008 Tel-Aviv University (Israel). 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/Minkowski_sum_2/include/CGAL/Minkowski_sum_2/Approx_offset_base_2.h $ 7 // $Id: Approx_offset_base_2.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 // Author(s) : Ron Wein <wein_r@yahoo.com> 11 // Andreas Fabri <Andreas.Fabri@geometryfactory.com> 12 // Laurent Rineau <Laurent.Rineau@geometryfactory.com> 13 // Efi Fogel <efif@post.tau.ac.il> 14 15 #ifndef CGAL_APPROXIMATED_OFFSET_BASE_H 16 #define CGAL_APPROXIMATED_OFFSET_BASE_H 17 18 #include <CGAL/license/Minkowski_sum_2.h> 19 20 21 #include <CGAL/Polygon_2.h> 22 #include <CGAL/Polygon_with_holes_2.h> 23 #include <CGAL/Gps_circle_segment_traits_2.h> 24 #include <CGAL/Minkowski_sum_2/Labels.h> 25 #include <CGAL/Minkowski_sum_2/Arr_labeled_traits_2.h> 26 #include <CGAL/use.h> 27 28 namespace CGAL { 29 30 /*! \class 31 * A base class for approximating the offset of a given polygon by a given 32 * radius. 33 */ 34 template <class Kernel_, class Container_> 35 class Approx_offset_base_2 36 { 37 private: 38 39 typedef Kernel_ Kernel; 40 typedef typename Kernel::FT NT; 41 42 protected: 43 44 typedef Kernel Basic_kernel; 45 typedef NT Basic_NT; 46 typedef CGAL::Polygon_2<Kernel, Container_> Polygon_2; 47 typedef CGAL::Polygon_with_holes_2<Kernel, Container_> Polygon_with_holes_2; 48 49 private: 50 51 // Kernel types: 52 typedef typename Kernel::Point_2 Point_2; 53 typedef typename Kernel::Line_2 Line_2; 54 55 // Polygon-related types: 56 typedef typename Polygon_2::Vertex_circulator Vertex_circulator; 57 58 // Traits-class types: 59 typedef Gps_circle_segment_traits_2<Kernel> Traits_2; 60 typedef typename Traits_2::CoordNT CoordNT; 61 typedef typename Traits_2::Point_2 Tr_point_2; 62 typedef typename Traits_2::Curve_2 Curve_2; 63 typedef typename Traits_2::X_monotone_curve_2 X_monotone_curve_2; 64 65 protected: 66 67 typedef typename Traits_2::Polygon_2 Offset_polygon_2; 68 69 typedef Arr_labeled_traits_2<Traits_2> Labeled_traits_2; 70 typedef typename Labeled_traits_2::X_monotone_curve_2 Labeled_curve_2; 71 72 // Data members: 73 double _eps; // An upper bound on the approximation error. 74 int _inv_sqrt_eps; // The inverse squared root of _eps. 75 76 public: 77 78 /*! 79 * Constructor. 80 * \param eps An upper bound on the approximation error. 81 */ Approx_offset_base_2(const double & eps)82 Approx_offset_base_2 (const double& eps) : 83 _eps (eps) 84 { 85 CGAL_precondition (CGAL::sign (eps) == POSITIVE); 86 87 _inv_sqrt_eps = static_cast<int> (1.0 / CGAL::sqrt (_eps)); 88 if (_inv_sqrt_eps <= 0) 89 _inv_sqrt_eps = 1; 90 } 91 92 protected: 93 94 /*! 95 * Compute curves that constitute the offset of a simple polygon by a given 96 * radius, with a given approximation error. 97 * \param pgn The polygon. 98 * \param orient The orientation to traverse the vertices. 99 * \param r The offset radius. 100 * \param cycle_id The index of the cycle. 101 * \param oi An output iterator for the offset curves. 102 * \pre The value type of the output iterator is Labeled_curve_2. 103 * \return A past-the-end iterator for the holes container. 104 */ 105 template <class OutputIterator> _offset_polygon(const Polygon_2 & pgn,CGAL::Orientation orient,const Basic_NT & r,unsigned int cycle_id,OutputIterator oi)106 OutputIterator _offset_polygon (const Polygon_2& pgn, 107 CGAL::Orientation orient, 108 const Basic_NT& r, 109 unsigned int cycle_id, 110 OutputIterator oi) const 111 { 112 // Prepare circulators over the polygon vertices. 113 const bool forward = (pgn.orientation() == orient); 114 Vertex_circulator first, curr, next, prev; 115 116 first = pgn.vertices_circulator(); 117 curr = first; 118 next = first; 119 prev = first; 120 121 if (forward) 122 --prev; 123 else 124 ++prev; 125 126 // Traverse the polygon vertices and edges and approximate the arcs that 127 // constitute the single convolution cycle. 128 NT x1, y1; // The source of the current edge. 129 NT x2, y2; // The target of the current edge. 130 NT delta_x, delta_y; // (x2 - x1) and (y2 - y1), resp. 131 NT abs_delta_x; 132 NT abs_delta_y; 133 CGAL::Sign sign_delta_x; // The sign of (x2 - x1). 134 CGAL::Sign sign_delta_y; // The sign of (y2 - y1). 135 NT sqr_d; // The squared length of the edge. 136 NT err_bound; // An approximation bound for d. 137 NT app_d; // The apporximated edge length. 138 NT app_err; // The approximation error. 139 CGAL::Sign sign_app_err; // Its sign. 140 NT lower_tan_half_phi; 141 NT upper_tan_half_phi; 142 NT sqr_tan_half_phi; 143 NT sin_phi, cos_phi; 144 Point_2 op1; // The approximated offset point 145 // corresponding to (x1, y1). 146 Point_2 op2; // The approximated offset point 147 // corresponding to (x2, y2). 148 Line_2 l1, l2; // Lines tangent at op1 and op2. 149 Object obj; 150 bool assign_success; 151 Point_2 mid_p; // The intersection of l1 and l2. 152 Point_2 first_op; // op1 for the first edge visited. 153 Point_2 prev_op; // op2 for the previous edge. 154 155 unsigned int curve_index = 0; 156 X_monotone_curve_2 seg1, seg2; 157 bool dir_right1 = false, dir_right2 = false; 158 X_monotone_curve_2 seg_short; 159 bool dir_right_short; 160 int n_segments; 161 162 Kernel ker; 163 typename Kernel::Intersect_2 f_intersect = ker.intersect_2_object(); 164 typename Kernel::Construct_line_2 f_line = ker.construct_line_2_object(); 165 typename Kernel::Construct_perpendicular_line_2 166 f_perp_line = ker.construct_perpendicular_line_2_object(); 167 typename Kernel::Compare_xy_2 f_comp_xy = ker.compare_xy_2_object(); 168 typename Kernel::Orientation_2 f_orient = ker.orientation_2_object(); 169 170 Traits_2 traits; 171 std::list<Object> xobjs; 172 std::list<Object>::iterator xobj_it; 173 typename Traits_2::Make_x_monotone_2 174 f_make_x_monotone = traits.make_x_monotone_2_object(); 175 Curve_2 arc; 176 X_monotone_curve_2 xarc; 177 178 do 179 { 180 // Get a circulator for the next vertex (in the proper orientation). 181 if (forward) 182 ++next; 183 else 184 --next; 185 186 // Compute the vector v = (delta_x, delta_y) of the current edge, 187 // and compute the squared edge length. 188 x1 = curr->x(); 189 y1 = curr->y(); 190 x2 = next->x(); 191 y2 = next->y(); 192 193 delta_x = x2 - x1; 194 delta_y = y2 - y1; 195 sqr_d = CGAL::square (delta_x) + CGAL::square (delta_y); 196 197 sign_delta_x = CGAL::sign (delta_x); 198 sign_delta_y = CGAL::sign (delta_y); 199 200 if (sign_delta_x == CGAL::ZERO) 201 { 202 CGAL_assertion (sign_delta_y != CGAL::ZERO); 203 204 // The edge [(x1, y1) -> (x2, y2)] is vertical. The offset edge lies 205 // at a distance r to the right if y2 > y1, and to the left if y2 < y1. 206 if (sign_delta_y == CGAL::POSITIVE) 207 { 208 op1 = Point_2 (x1 + r, y1); 209 op2 = Point_2 (x2 + r, y2); 210 } 211 else 212 { 213 op1 = Point_2 (x1 - r, y1); 214 op2 = Point_2 (x2 - r, y2); 215 } 216 217 // Create the offset segment [op1 -> op2]. 218 seg1 = X_monotone_curve_2 (op1, op2); 219 dir_right1 = (sign_delta_y == CGAL::POSITIVE); 220 221 n_segments = 1; 222 } 223 else if (sign_delta_y == CGAL::ZERO) 224 { 225 // The edge [(x1, y1) -> (x2, y2)] is horizontal. The offset edge lies 226 // at a distance r to the bottom if x2 > x1, and to the top if x2 < x1. 227 if (sign_delta_x == CGAL::POSITIVE) 228 { 229 op1 = Point_2 (x1, y1 - r); 230 op2 = Point_2 (x2, y2 - r); 231 } 232 else 233 { 234 op1 = Point_2 (x1, y1 + r); 235 op2 = Point_2 (x2, y2 + r); 236 } 237 238 // Create the offset segment [op1 -> op2]. 239 seg1 = X_monotone_curve_2 (op1, op2); 240 dir_right1 = (sign_delta_x == CGAL::POSITIVE); 241 242 n_segments = 1; 243 } 244 else 245 { 246 abs_delta_x = (sign_delta_x == POSITIVE) ? delta_x : NT(-delta_x); 247 abs_delta_y = (sign_delta_y == POSITIVE) ? delta_y : NT(-delta_y); 248 249 // In this general case, the length d of the current edge is usually 250 // an irrational number. 251 // Compute the upper bound for the approximation error for d. 252 // This bound is given by: 253 // 254 // d - |delta_y| 255 // bound = 2 * d * eps * --------------- 256 // |delta_x| 257 // 258 // As we use floating-point arithmetic, if |delta_x| is small, then 259 // it might be that to_double(|delta_y|) == to_double(d), hence we 260 // have a 0 tolerance in the approximation bound. Luckily, because 261 // of symmetry, we can rotate the scene by pi/2, and swap roles of 262 // x and y. In fact, we do that in order to get a larger approximation 263 // bound if possible. 264 const double dd = CGAL::sqrt (CGAL::to_double (sqr_d)); 265 const double dabs_dx = CGAL::to_double (abs_delta_x); 266 const double dabs_dy = CGAL::to_double (abs_delta_y); 267 double derr_bound; 268 269 if (dabs_dy < dabs_dx) 270 { 271 derr_bound = 2 * dd * _eps * (dd - dabs_dy) / dabs_dx; 272 } 273 else 274 { 275 derr_bound = 2 * dd * _eps * (dd - dabs_dx) / dabs_dy; 276 } 277 278 CGAL_assertion (derr_bound > 0); 279 err_bound = NT (derr_bound); 280 281 // Compute an approximation for d (the squared root of sqr_d). 282 int numer; 283 int denom = _inv_sqrt_eps; 284 const int max_int = (1 << (8*sizeof(int) - 2)); 285 286 numer = static_cast<int> (dd * denom + 0.5); 287 if (numer > 0) 288 { 289 while (static_cast<double>(max_int) / denom < dd && 290 numer > 0) 291 { 292 denom >>= 1; 293 numer = static_cast<int> (dd * denom + 0.5); 294 } 295 } 296 else if (numer == 0) 297 { 298 while (numer == 0) 299 { 300 denom <<= 1; 301 if (denom > 0) 302 { 303 numer = static_cast<int> (dd * denom + 0.5); 304 } 305 else 306 { 307 // In case of overflow of denom 308 numer = 1; 309 denom = max_int; 310 } 311 } 312 } 313 else {// if numer < 0 (overflow) 314 numer = max_int; 315 denom = 1; 316 } 317 318 319 app_d = NT (numer) / NT (denom); 320 app_err = sqr_d - CGAL::square (app_d); 321 322 while (CGAL::compare (CGAL::abs (app_err), 323 err_bound) == CGAL::LARGER || 324 CGAL::compare (app_d, abs_delta_x) != LARGER || 325 CGAL::compare (app_d, abs_delta_y) != LARGER) 326 { 327 app_d = (app_d + sqr_d/app_d) / 2; 328 app_err = sqr_d - CGAL::square (app_d); 329 } 330 331 sign_app_err = CGAL::sign (app_err); 332 333 if (sign_app_err == CGAL::ZERO) 334 { 335 // In this case d is a rational number, and we should shift the 336 // both edge endpoints by (r * delta_y / d, -r * delta_x / d) to 337 // obtain the offset points op1 and op2. 338 const NT trans_x = r * delta_y / app_d; 339 const NT trans_y = r * (-delta_x) / app_d; 340 341 op1 = Point_2 (x1 + trans_x, y1 + trans_y); 342 op2 = Point_2 (x2 + trans_x, y2 + trans_y); 343 344 seg1 = X_monotone_curve_2 (op1, op2); 345 dir_right1 = (sign_delta_x == CGAL::POSITIVE); 346 347 n_segments = 1; 348 } 349 else 350 { 351 // In case |x2 - x1| < |y2 - y1| (and phi is small) it is possible 352 // that the approximation t' of t = tan(phi/2) is of opposite sign. 353 // To avoid this problem, we symbolically rotate the scene by pi/2, 354 // swapping roles between x and y. Thus, t is not close to zero, and 355 // we are guaranteed to have: phi- < phi < phi+ . 356 bool rotate_pi2 = false; 357 358 if (CGAL::compare (CGAL::abs(delta_x), 359 CGAL::abs(delta_y)) == SMALLER) 360 { 361 rotate_pi2 = true; 362 363 // We use the rotation matrix by pi/2: 364 // 365 // +- -+ 366 // | 0 -1 | 367 // | 1 0 | 368 // +- -+ 369 // 370 // Thus, the point (x, y) is converted to (-y, x): 371 NT tmp = x1; 372 373 x1 = -y1; 374 y1 = tmp; 375 376 tmp = x2; 377 x2 = -y2; 378 y2 = tmp; 379 380 // Swap the delta_x and delta_y values. 381 tmp = delta_x; 382 delta_x = -delta_y; 383 delta_y = tmp; 384 385 CGAL::Sign tmp_sign = sign_delta_x; 386 sign_delta_x = CGAL::opposite (sign_delta_y); 387 sign_delta_y = tmp_sign; 388 } 389 390 // Act according to the sign of delta_x. 391 if (sign_delta_x == CGAL::NEGATIVE) 392 { 393 // x1 > x2, so we take a lower approximation for the squared root. 394 if (sign_app_err == CGAL::NEGATIVE) 395 app_d = sqr_d / app_d; 396 } 397 else 398 { 399 // x1 < x2, so we take an upper approximation for the squared root. 400 if (sign_app_err == CGAL::POSITIVE) 401 app_d = sqr_d / app_d; 402 } 403 404 // If theta is the angle that the vector (delta_x, delta_y) forms 405 // with the x-axis, the perpendicular vector forms an angle of 406 // phi = theta - PI/2, and we can approximate tan(phi/2) from below 407 // and from above using: 408 lower_tan_half_phi = (app_d - delta_y) / (-delta_x); 409 upper_tan_half_phi = (-delta_x) / (app_d + delta_y); 410 411 // Translate (x1, y1) by (r*cos(phi-), r*sin(phi-)) and create the 412 // first offset point. 413 // If tan(phi/2) = t is rational, then sin(phi) = 2t/(1 + t^2) 414 // and cos(phi) = (1 - t^2)/(1 + t^2) are also rational. 415 sqr_tan_half_phi = CGAL::square (lower_tan_half_phi); 416 sin_phi = 2 * lower_tan_half_phi / (1 + sqr_tan_half_phi); 417 cos_phi = (1 - sqr_tan_half_phi) / (1 + sqr_tan_half_phi); 418 419 if (! rotate_pi2) 420 { 421 op1 = Point_2 (x1 + r*cos_phi, y1 + r*sin_phi); 422 } 423 else 424 { 425 // In case of symbolic rotation by pi/2, we have to rotate the 426 // translated point by -(pi/2), transforming (x, y) to (y, -x). 427 op1 = Point_2 (y1 + r*sin_phi, -(x1 + r*cos_phi)); 428 } 429 430 // Translate (x2, y2) by (r*cos(phi+), r*sin(phi+)) and create the 431 // second offset point. 432 sqr_tan_half_phi = CGAL::square (upper_tan_half_phi); 433 sin_phi = 2 * upper_tan_half_phi / (1 + sqr_tan_half_phi); 434 cos_phi = (1 - sqr_tan_half_phi) / (1 + sqr_tan_half_phi); 435 436 if (! rotate_pi2) 437 { 438 op2 = Point_2 (x2 + r*cos_phi, y2 + r*sin_phi); 439 } 440 else 441 { 442 // In case of symbolic rotation by pi/2, we have to rotate the 443 // translated point by -(pi/2), transforming (x, y) to (y, -x). 444 op2 = Point_2 (y2 + r*sin_phi, -(x2 + r*cos_phi)); 445 } 446 447 // Compute the line l1 tangent to the circle centered at (x1, y1) 448 // with radius r at the approximated point op1. 449 l1 = f_perp_line (f_line (*curr, op1), op1); 450 451 // Compute the line l2 tangent to the circle centered at (x2, y2) 452 // with radius r at the approximated point op2. 453 l2 = f_perp_line (f_line (*next, op2), op2); 454 455 // Intersect the two lines. The intersection point serves as a common 456 // end point for the two line segments we are about to introduce. 457 obj = f_intersect (l1, l2); 458 459 assign_success = CGAL::assign (mid_p, obj); 460 CGAL_assertion (assign_success); 461 CGAL_USE(assign_success); 462 463 // Andreas's assertions: 464 CGAL_assertion( right_turn(*curr, *next, op2) ); 465 CGAL_assertion( angle(*curr, *next, op2) != ACUTE); 466 CGAL_assertion( angle(op1, *curr, *next) != ACUTE); 467 CGAL_assertion( right_turn(op1, *curr, *next) ); 468 469 // Create the two segments [op1 -> p_mid] and [p_min -> op2]. 470 seg1 = X_monotone_curve_2 (op1, mid_p); 471 dir_right1 = (f_comp_xy (op1, mid_p) == CGAL::SMALLER); 472 473 seg2 = X_monotone_curve_2 (mid_p, op2); 474 dir_right2 = (f_comp_xy (mid_p, op2) == CGAL::SMALLER); 475 476 n_segments = 2; 477 } 478 } 479 480 if (curr == first) { 481 // This is the first edge we visit -- store op1 for future use. 482 first_op = op1; 483 } 484 else { 485 CGAL::Orientation orient = f_orient (*prev, *curr, *next); 486 if (orient == CGAL::COLLINEAR) { 487 /* If the orientation is collinear, figure out whether it's a 180 488 * turn. If so, assume that it is an antena that generates 489 * a positive spike, and treat it as a left turn. 490 * A complete solution would need to distinguish between positive 491 * and negative spikes, and treat them as left and right turns, 492 * respectively. 493 */ 494 typename Kernel::Compare_x_2 f_comp_x = ker.compare_x_2_object(); 495 Comparison_result res1, res2; 496 res1 = f_comp_x(*prev, *curr); 497 if (res1 != CGAL::EQUAL) 498 res2 = f_comp_x(*curr, *next); 499 else { 500 typename Kernel::Compare_y_2 f_comp_y = ker.compare_y_2_object(); 501 res1 = f_comp_y(*prev, *curr); 502 res2 = f_comp_y(*curr, *next); 503 } 504 if (res1 != res2) orient = CGAL::LEFT_TURN; 505 } 506 507 // Connect the offset target point of the previous edge to the 508 // offset source of the current edge. 509 if (orient == CGAL::LEFT_TURN) { 510 // Connect prev_op and op1 with a circular arc, whose supporting 511 // circle is (x1, x2) with radius r. 512 arc = Curve_2 (*curr, r, CGAL::COUNTERCLOCKWISE, 513 Tr_point_2 (prev_op.x(), prev_op.y()), 514 Tr_point_2 (op1.x(), op1.y())); 515 516 // Subdivide the arc into x-monotone subarcs and insert them into the 517 // convolution cycle. 518 xobjs.clear(); 519 f_make_x_monotone (arc, std::back_inserter(xobjs)); 520 for (xobj_it = xobjs.begin(); xobj_it != xobjs.end(); ++xobj_it) { 521 assign_success = CGAL::assign (xarc, *xobj_it); 522 CGAL_assertion (assign_success); 523 CGAL_USE(assign_success); 524 *oi++ = Labeled_curve_2 (xarc, 525 X_curve_label (xarc.is_directed_right(), 526 cycle_id, curve_index++)); 527 } 528 } 529 else if (orient == CGAL::RIGHT_TURN) { 530 // In case the current angle between the previous and the current 531 // edge is larger than pi/2, it not necessary to connect prev_op 532 // and op1 by a circular arc (as the case above): it is sufficient 533 // to shortcut the circular arc using a segment, whose sole purpose 534 // is to guarantee the continuity of the convolution cycle (we know 535 // this segment will not be part of the output offset or inset). 536 seg_short = X_monotone_curve_2(prev_op, op1); 537 538 dir_right_short = (f_comp_xy (prev_op, op1) == CGAL::SMALLER); 539 *oi++ = Labeled_curve_2 (seg_short, 540 X_curve_label (dir_right_short, 541 cycle_id, curve_index++)); 542 } 543 } 544 545 // Append the offset segment(s) to the convolution cycle. 546 CGAL_assertion (n_segments == 1 || n_segments == 2); 547 *oi++ = Labeled_curve_2 (seg1, X_curve_label (dir_right1, 548 cycle_id, curve_index++)); 549 550 if (n_segments == 2) 551 { 552 *oi++ = Labeled_curve_2 (seg2, X_curve_label (dir_right2, 553 cycle_id, curve_index++)); 554 } 555 556 // Proceed to the next polygon vertex. 557 prev_op = op2; 558 prev = curr; 559 curr = next; 560 561 } while (curr != first); 562 563 // Close the convolution cycle by creating the final circular arc, 564 // centered at the first vertex. 565 arc = Curve_2 (*first, r, CGAL::COUNTERCLOCKWISE, 566 Tr_point_2 (op2.x(), op2.y()), 567 Tr_point_2 (first_op.x(), first_op.y())); 568 569 // Subdivide the arc into x-monotone subarcs and insert them to the 570 // convolution cycle. 571 xobjs.clear(); 572 f_make_x_monotone (arc, std::back_inserter(xobjs)); 573 574 xobj_it = xobjs.begin(); 575 while (xobj_it != xobjs.end()) 576 { 577 assign_success = CGAL::assign (xarc, *xobj_it); 578 CGAL_assertion (assign_success); 579 CGAL_USE(assign_success); 580 581 ++xobj_it; 582 bool is_last = (xobj_it == xobjs.end()); 583 584 *oi++ = Labeled_curve_2 (xarc, 585 X_curve_label (xarc.is_directed_right(), 586 cycle_id, curve_index++, is_last)); 587 } 588 589 return (oi); 590 } 591 592 }; 593 594 } //namespace CGAL 595 596 #endif 597