1 // Copyright (c) 1998-2003 ETH Zurich (Switzerland).
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/Bounding_volumes/include/CGAL/pierce_rectangles_2.h $
7 // $Id: pierce_rectangles_2.h bf325bf 2021-01-06T10:54:49+01:00 Sébastien Loriot
8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
9 //
10 //
11 // Author(s) : Michael Hoffmann <hoffmann@inf.ethz.ch>
12
13 #ifndef CGAL_PIERCE_RECTANGLES_2_H
14 #define CGAL_PIERCE_RECTANGLES_2_H 1
15
16 #include <CGAL/license/Bounding_volumes.h>
17
18
19 #include <CGAL/Optimisation/assertions.h>
20 #include <CGAL/circulator.h>
21 #include <CGAL/algorithm.h>
22 #include <algorithm>
23 #include <iterator>
24 #include <vector>
25
26 #if defined(BOOST_MSVC)
27 # pragma warning(push)
28 # pragma warning(disable:4355) // complaint about using 'this' to
29 #endif // initialize a member
30
31 namespace CGAL {
32
33 //!!! STL-extensions
34 template < class T >
35 struct Wastebasket {
36 typedef std::output_iterator_tag iterator_category;
37 typedef Wastebasket< T > iterator;
38
39 iterator
40 operator=( const T&)
41 { return *this; }
42
43 iterator
44 operator*()
45 { return *this; }
46
47 iterator
48 operator++()
49 { return *this; }
50
51 iterator
52 operator++( int)
53 { return *this; }
54 };
55
56 template < class Traits_ >
57 struct Loc_domain {
58 // ---------------------------------------------
59 // types:
60
61 typedef Traits_ Traits;
62 typedef typename Traits::FT FT;
63 typedef typename Traits::Point_2 Point_2;
64
65 typedef std::vector< Point_2 > Container;
66 typedef typename Container::iterator Iterator;
67 typedef typename Container::const_iterator Citerator;
68 typedef typename Container::reverse_iterator Riterator;
69 typedef typename Container::const_reference Creference;
70 typedef typename Container::size_type size_type;
71
72 // ---------------------------------------------
73 // creation:
74
75 void
updateLoc_domain76 update(int j, Citerator i)
77 {
78 CGAL_optimisation_precondition(j >= 0 && j < 4);
79 if (j < 2)
80 if (j == 0) {
81 if (traits.less_x_2_object()(*i, minx)) minx = *i;
82 if (traits.less_y_2_object()(*i, miny)) miny = *i;
83 }
84 else {
85 if (traits.less_y_2_object()(*i, miny)) miny = *i;
86 if (traits.less_x_2_object()(maxx, *i)) maxx = *i;
87 }
88 else
89 if (j == 2) {
90 if (traits.less_x_2_object()(maxx, *i)) maxx = *i;
91 if (traits.less_y_2_object()(maxy, *i)) maxy = *i;
92 }
93 else {
94 if (traits.less_y_2_object()(maxy, *i)) maxy = *i;
95 if (traits.less_x_2_object()(*i, minx)) minx = *i;
96 }
97 }
98
99 template < class InputIC >
Loc_domainLoc_domain100 Loc_domain(InputIC b, InputIC e, Traits t)
101 : pts(b, e),
102 end_(pts.end()),
103 minx(pts.front()),
104 miny(pts.front()),
105 maxx(pts.front()),
106 maxy(pts.front()),
107 traits(t)
108 {
109 CGAL_optimisation_precondition(b != e);
110 Iterator i = pts.begin();
111 CGAL_optimisation_assertion(i != pts.end());
112 while (++i != pts.end()) {
113 if (traits.less_x_2_object()(*i, minx)) minx = *i;
114 if (traits.less_x_2_object()(maxx, *i)) maxx = *i;
115 if (traits.less_y_2_object()(*i, miny)) miny = *i;
116 if (traits.less_y_2_object()(maxy, *i)) maxy = *i;
117 }
118 }
119
120 // ---------------------------------------------
121 // access operations:
122
123 Point_2
124 operator[](int i) const
125 // return corner points (0 <-> bottom-left, 1 <-> bottom-right)
126 {
127 CGAL_optimisation_precondition(i >= 0 && i < 4);
128 if (i == 0)
129 return traits.construct_point_2_above_right_implicit_point_2_object()(
130 minx, miny, r);
131 else if (i == 1)
132 return traits.construct_point_2_above_left_implicit_point_2_object()(
133 maxx, miny, r);
134 else if (i == 2)
135 return traits.construct_point_2_below_left_implicit_point_2_object()(
136 maxx, maxy, r);
137 return traits.construct_point_2_below_right_implicit_point_2_object()(
138 minx, maxy, r);
139 }
140
141 Point_2
extremeLoc_domain142 extreme(int i) const
143 // return extreme points (0 <-> left, 1 <-> bottom)
144 {
145 CGAL_optimisation_precondition(i >= 0 && i < 4);
146 if (i > 1) return i == 2 ? maxx : maxy;
147 return i == 0 ? minx : miny;
148 }
149
150 Point_2&
extremeLoc_domain151 extreme(int i)
152 // return extreme points (0 <-> left, 1 <-> bottom)
153 {
154 CGAL_optimisation_precondition(i >= 0 && i < 4);
155 if (i > 1) return i == 2 ? maxx : maxy;
156 return i == 0 ? minx : miny;
157 }
158
beginLoc_domain159 Citerator begin() const { return pts.begin(); }
beginLoc_domain160 Iterator begin() { return pts.begin(); }
161 #ifndef _MSC_VER
endLoc_domain162 Citerator end() const { return end_; }
163 #else
164 // Yet another really great MSVC feature ...
165 // without that static_cast to itself it does not compile :-)
endLoc_domain166 Citerator end() const { return static_cast<Iterator>(end_); }
167 #endif
endLoc_domain168 Iterator& end() { return end_; }
real_endLoc_domain169 Iterator real_end() { return pts.end(); }
170
emptyLoc_domain171 bool empty() const { return begin() == end(); }
sizeLoc_domain172 size_type size() const { return end() - begin(); }
frontLoc_domain173 Creference front() const { return pts.front(); }
174
175 // ---------------------------------------------
176 // check operation:
177
178 void
checkLoc_domain179 check() const {
180 CGAL_optimisation_expensive_assertion_code(
181 Iterator i = pts.begin();
182 do {
183 CGAL_optimisation_assertion(!traits.less_x_2_object()(*i, minx));
184 CGAL_optimisation_assertion(!traits.less_x_2_object()(maxx, *i));
185 CGAL_optimisation_assertion(!traits.less_y_2_object()(*i, miny));
186 CGAL_optimisation_assertion(!traits.less_y_2_object()(maxy, *i));
187 } while (++i != end);
188 )
189 }
190
191 protected:
192 // points container
193 Container pts;
194 // const iterator to past-the-end
195 Iterator end_;
196 public:
197 // actual center radius
198 FT r;
199 // (copies of) elements with minimal/maximal x/y coordinate
200 Point_2 minx, miny, maxx, maxy;
201 // Traits class
202 Traits traits;
203
204 }; // class Loc_domain
205 template < class Traits_ >
206 struct Staircases : public Loc_domain< Traits_ > {
207 typedef Traits_ Traits;
208 typedef Loc_domain< Traits > Base;
209 typedef typename Base::Container Container;
210 typedef typename Base::Iterator Iterator;
211 typedef typename Base::Citerator Citerator;
212 typedef typename Base::Riterator Riterator;
213 typedef typename Traits::Point_2 Point_2;
214 typedef typename Traits::FT FT;
215 typedef std::pair< Point_2, Point_2 > Intervall;
216
217 template < class InputIC >
StaircasesStaircases218 Staircases(InputIC b, InputIC e, Traits t)
219 : Base(b, e, t),
220 sorted(this->pts),
221 xgy(t.signed_x_distance_2_object()(this->maxx, this->minx) >
222 t.signed_y_distance_2_object()(this->maxy, this->miny))
223 {
224 using std::sort;
225 using std::find_if;
226
227 Container& xsort = xgy ? sorted : this->pts;
228 Container& ysort = xgy ? this->pts : sorted;
229
230 // build top-left and bottom-right staircases
231 sort(ysort.begin(), ysort.end(), this->traits.less_y_2_object());
232 // bottom-right
233 Iterator i = ysort.begin();
234 do {
235 brstc.push_back(*i++);
236 i = find_if(i, ysort.end(),
237 [this](const Point_2& p)
238 { return this->traits.less_x_2_object()(this->brstc.back(), p);});
239 } while (i != ysort.end());
240 // top-left
241 Riterator j = ysort.rbegin();
242 do {
243 tlstc.push_back(*j++);
244 j = find_if(j, ysort.rend(),
245 [this](const Point_2& p)
246 { return this->traits.less_x_2_object()(p, this->tlstc.back());});
247 } while (j != ysort.rend());
248
249 // build left-bottom and right-top staircases
250 sort(xsort.begin(), xsort.end(), this->traits.less_x_2_object());
251 // left-bottom
252 i = xsort.begin();
253 do {
254 lbstc.push_back(*i++);
255 i = find_if(i, xsort.end(),
256 [this](const Point_2& p)
257 { return this->traits.less_y_2_object()(p, this->lbstc.back());});
258 } while (i != xsort.end());
259 // right-top
260 j = xsort.rbegin();
261 do {
262 rtstc.push_back(*j++);
263 j = find_if(j, xsort.rend(),
264 [this](const Point_2& p)
265 { return this->traits.less_y_2_object()(this->rtstc.back(), p);});
266 } while (j != xsort.rend());
267 } // Staircases(b, e, t)
268
is_middle_emptyStaircases269 bool is_middle_empty() const {
270 //!!! the "middle" point could be precomputed in advance
271 Citerator i = this->pts.begin();
272 FT rr = FT(2) * this->r;
273 do
274 if (this->traits.signed_x_distance_2_object()(this->maxx, *i) > rr &&
275 this->traits.signed_x_distance_2_object()(*i, this->minx) > rr &&
276 this->traits.signed_y_distance_2_object()(*i, this->miny) > rr &&
277 this->traits.signed_y_distance_2_object()(this->maxy, *i) > rr)
278 return false;
279 while (++i != this->pts.end());
280 return true;
281 } // is_middle()
282
is_x_greater_yStaircases283 bool is_x_greater_y() const { return xgy; }
284
tlstc_beginStaircases285 Citerator tlstc_begin() const { return tlstc.begin(); }
tlstc_endStaircases286 Citerator tlstc_end() const { return tlstc.end(); }
lbstc_beginStaircases287 Citerator lbstc_begin() const { return lbstc.begin(); }
lbstc_endStaircases288 Citerator lbstc_end() const { return lbstc.end(); }
brstc_beginStaircases289 Citerator brstc_begin() const { return brstc.begin(); }
brstc_endStaircases290 Citerator brstc_end() const { return brstc.end(); }
rtstc_beginStaircases291 Citerator rtstc_begin() const { return rtstc.begin(); }
rtstc_endStaircases292 Citerator rtstc_end() const { return rtstc.end(); }
293
top_intervallStaircases294 Intervall top_intervall() const {
295 Point_2 p =
296 this->traits.construct_point_2_above_right_implicit_point_2_object()(
297 this->minx, this->miny, FT(2) * this->r);
298 Point_2 q =
299 this->traits.construct_point_2_above_left_implicit_point_2_object()(
300 this->maxx, this->miny, FT(2) * this->r);
301
302 Citerator i =
303 min_element_if(
304 this->pts.begin(), this->pts.end(),
305 this->traits.less_x_2_object(),
306 [&p, this](const Point_2& pt)
307 {
308 return this->traits.less_x_2_object()(p, pt) &&
309 this->traits.less_y_2_object()(p, pt);
310 });
311 Citerator j =
312 max_element_if(
313 this->pts.begin(), this->pts.end(),
314 this->traits.less_x_2_object(),
315 [&q, this](const Point_2& pt)
316 {
317 return this->traits.less_x_2_object()(pt, q) &&
318 this->traits.less_y_2_object()(q, pt);
319 });
320 return Intervall(i == this->pts.end() ? this->maxx : *i,
321 j == this->pts.end() ? this->minx : *j);
322 } // top_intervall()
323
bottom_intervallStaircases324 Intervall bottom_intervall() const {
325 Point_2 p =
326 this->traits.construct_point_2_below_right_implicit_point_2_object()(
327 this->minx, this->maxy, FT(2) * this->r);
328 Point_2 q =
329 this->traits.construct_point_2_below_left_implicit_point_2_object()(
330 this->maxx, this->maxy, FT(2) * this->r);
331
332 Citerator i =
333 min_element_if(
334 this->pts.begin(), this->pts.end(),
335 this->traits.less_x_2_object(),
336 [&p, this](const Point_2& pt)
337 {
338 return this->traits.less_x_2_object()(p, pt) &&
339 this->traits.less_y_2_object()(pt, p);
340 });
341 Citerator j =
342 max_element_if(
343 this->pts.begin(), this->pts.end(),
344 this->traits.less_x_2_object(),
345 [&q, this](const Point_2& pt)
346 {
347 return this->traits.less_x_2_object()(pt, q) &&
348 this->traits.less_y_2_object()(pt, q);
349 });
350 return Intervall(i == this->pts.end() ? this->maxx : *i,
351 j == this->pts.end() ? this->minx : *j);
352 } // bottom_intervall()
353
left_intervallStaircases354 Intervall left_intervall() const {
355 Point_2 p =
356 this->traits.construct_point_2_above_left_implicit_point_2_object()(
357 this->maxx, this->miny, FT(2) * this->r);
358 Point_2 q =
359 this->traits.construct_point_2_below_left_implicit_point_2_object()(
360 this->maxx, this->maxy, FT(2) * this->r);
361
362 Citerator i =
363 min_element_if(
364 this->pts.begin(), this->pts.end(),
365 this->traits.less_y_2_object(),
366 [&p, this](const Point_2& pt)
367 {
368 return this->traits.less_x_2_object()(pt, p) &&
369 this->traits.less_y_2_object()(p, pt);
370 });
371 Citerator j =
372 max_element_if(
373 this->pts.begin(), this->pts.end(),
374 this->traits.less_y_2_object(),
375 [&q, this](const Point_2& pt)
376 {
377 return this->traits.less_x_2_object()(pt, q) &&
378 this->traits.less_y_2_object()(pt, q);
379 });
380 return Intervall(i == this->pts.end() ? this->maxy : *i,
381 j == this->pts.end() ? this->miny : *j);
382 } // left_intervall()
383
right_intervallStaircases384 Intervall right_intervall() const {
385 Point_2 p =
386 this->traits.construct_point_2_above_right_implicit_point_2_object()(
387 this->minx, this->miny, FT(2) * this->r);
388 Point_2 q =
389 this->traits.construct_point_2_below_right_implicit_point_2_object()(
390 this->minx, this->maxy, FT(2) * this->r);
391
392 Citerator i =
393 min_element_if(
394 this->pts.begin(), this->pts.end(),
395 this->traits.less_y_2_object(),
396 [&p, this](const Point_2& pt)
397 {
398 return this->traits.less_x_2_object()(p, pt) &&
399 this->traits.less_y_2_object()(p, pt);
400 });
401 Citerator j =
402 max_element_if(
403 this->pts.begin(), this->pts.end(),
404 this->traits.less_y_2_object(),
405 [&q, this](const Point_2& pt)
406 {
407 return this->traits.less_x_2_object()(q, pt) &&
408 this->traits.less_y_2_object()(pt, q);
409 });
410 return Intervall(i == this->pts.end() ? this->maxy : *i,
411 j == this->pts.end() ? this->miny : *j);
412 } // right_intervall()
413
414 template < class OutputIterator >
shared_intervallStaircases415 OutputIterator shared_intervall(OutputIterator o) const {
416 if (xgy) {
417 if (this->traits.signed_y_distance_2_object()(
418 this->maxy, this->miny) > FT(4) * this->r)
419 return o;
420 Point_2 p =
421 this->traits.construct_point_2_below_right_implicit_point_2_object()(
422 this->minx, this->maxy, FT(2) * this->r);
423 Point_2 q =
424 this->traits.construct_point_2_above_left_implicit_point_2_object()(
425 this->maxx, this->miny, FT(2) * this->r);
426 //!!! start with binary search
427 for (Citerator i = sorted.begin(); i != sorted.end(); ++i)
428 if (this->traits.less_x_2_object()(p, *i) &&
429 this->traits.less_x_2_object()(*i, q) &&
430 !this->traits.less_y_2_object()(*i, p) &&
431 !this->traits.less_y_2_object()(q, *i))
432 *o++ = *i;
433 } else {
434 if (this->traits.signed_x_distance_2_object()(
435 this->maxx, this->minx) > FT(4) * this->r)
436 return o;
437 Point_2 p =
438 this->traits.construct_point_2_above_left_implicit_point_2_object()(
439 this->maxx, this->miny, FT(2) * this->r);
440 Point_2 q =
441 this->traits.construct_point_2_below_right_implicit_point_2_object()(
442 this->minx, this->maxy, FT(2) * this->r);
443 //!!! start with binary search
444 for (Citerator i = sorted.begin(); i != sorted.end(); ++i)
445 if (!this->traits.less_x_2_object()(*i, p) &&
446 !this->traits.less_x_2_object()(q, *i) &&
447 this->traits.less_y_2_object()(p, *i) &&
448 this->traits.less_y_2_object()(*i, q))
449 *o++ = *i;
450 }
451 return o;
452 } // shared_intervall(o)
453
454 private:
455 Container tlstc, lbstc, brstc, rtstc, sorted;
456 // exceeds the x-dimension of the location domain its y-dimension?
457 bool xgy;
458 };
459
460
461 template < class InputIC, class OutputIterator, class Traits >
462 inline OutputIterator
two_cover_points(InputIC f,InputIC l,OutputIterator o,bool & ok,const Traits & t)463 two_cover_points(
464 InputIC f, InputIC l, OutputIterator o, bool& ok, const Traits& t)
465 {
466 CGAL_optimisation_precondition(f != l);
467
468 // compute location domain:
469 Loc_domain< Traits > d(f, l, t);
470
471 return two_cover_points(d, o, ok, t);
472 } // two_cover_points(f, l, o, ok, t)
473 template < class InputIC, class OutputIterator, class Traits >
474 inline OutputIterator
three_cover_points(InputIC f,InputIC l,OutputIterator o,bool & ok,const Traits & t)475 three_cover_points(
476 InputIC f, InputIC l, OutputIterator o, bool& ok, const Traits& t)
477 {
478 CGAL_optimisation_precondition(f != l);
479
480 // compute location domain:
481 Loc_domain< Traits > d(f, l, t);
482
483 return three_cover_points(d, o, ok, t);
484 } // three_cover_points(f, l, o, ok, t)
485 template < class InputIC, class OutputIterator, class Traits >
486 inline OutputIterator
four_cover_points(InputIC f,InputIC l,OutputIterator o,bool & ok,const Traits & t)487 four_cover_points(
488 InputIC f, InputIC l, OutputIterator o, bool& ok, const Traits& t)
489 {
490 CGAL_optimisation_precondition(f != l);
491
492 // compute location domain:
493 Loc_domain< Traits > d(f, l, t);
494
495 return four_cover_points(d, o, ok, t);
496 } // four_cover_points(f, l, o, ok, t)
497
498 template < class OutputIterator, class Traits >
499 OutputIterator
two_cover_points(const Loc_domain<Traits> & d,OutputIterator o,bool & ok)500 two_cover_points(
501 const Loc_domain< Traits >& d,
502 OutputIterator o,
503 bool& ok)
504 {
505 using std::find_if;
506 using std::less;
507
508 typedef typename Traits::FT FT;
509 typedef typename Traits::Point_2 Point_2;
510 typename Traits::Infinity_distance_2 dist =
511 d.traits.infinity_distance_2_object();
512 typename Traits::Signed_infinity_distance_2 sdist =
513 d.traits.signed_infinity_distance_2_object();
514
515 if (sdist(d[2], d[0]) <= FT(0)) {
516 // the location domain is degenerate and [f,l) is one-pierceable
517 *o++ = d[0];
518 ok = true;
519 return o;
520 }
521
522
523 // check whether {d[0], d[2]} forms a piercing set
524 if (d.end() ==
525 find_if(d.begin(),
526 d.end(),
527 [&dist,&d](const Point_2& p)
528 { return d.r < Min<FT>()(dist(d[0], p), dist(d[2],p)); }))
529 {
530 *o++ = d[0];
531 *o++ = d[2];
532 ok = true;
533 return o;
534 }
535
536 // check whether {d[1], d[3]} forms a piercing set
537 if (d.end() ==
538 find_if(d.begin(),
539 d.end(),
540 [&dist,&d](const Point_2& p)
541 { return d.r < Min<FT>()(dist(d[1], p), dist(d[3],p)); }))
542 {
543 *o++ = d[1];
544 *o++ = d[3];
545 ok = true;
546 return o;
547 }
548
549 // no piercing set exists:
550 ok = false;
551 return o;
552 } // two_cover_points(d, o, ok, t)
553
554 template < class OutputIterator, class Traits >
555 OutputIterator
three_cover_points(Loc_domain<Traits> & d,OutputIterator o,bool & ok)556 three_cover_points(
557 Loc_domain< Traits >& d,
558 OutputIterator o,
559 bool& ok)
560 {
561 using std::find_if;
562 using std::less;
563 using std::iter_swap;
564
565 CGAL_optimisation_precondition(!d.empty());
566
567 // typedefs:
568 typedef typename Traits::Point_2 Point_2;
569 typedef typename Loc_domain< Traits >::Iterator Iterator;
570 typename Traits::Infinity_distance_2 dist =
571 d.traits.infinity_distance_2_object();
572
573 // test the four corners:
574 for (int k = 0; k < 4; ++k) {
575
576 // extract all points which are close enough to d[k]
577 Point_2 corner = d[k];
578
579 // find first point not covered by the rectangle at d[k]
580 Iterator i = find_if(d.begin(), d.end(),
581 [&d,&dist, &corner](const Point_2& p)
582 { return d.r < dist(corner, p); });
583
584 // are all points already covered?
585 if (i == d.end()) {
586 CGAL_optimisation_assertion(k == 0);
587 *o++ = d[0];
588 ok = true;
589 return o;
590 } // if (i == d.end())
591
592 // save changing sides of d:
593 Point_2 save_side1 = d.extreme(k);
594 Point_2 save_side2 = d.extreme((k+1) % 4);
595 Iterator save_end = d.end();
596
597 // now run through it:
598 // initialize the two (possibly) changing sides of d
599 d.extreme(k) = d.extreme((k+1) % 4) = *i;
600
601 // is there any point covered?
602 if (i == d.begin()) {
603 while (++i != d.end() && dist(corner, *i) > d.r)
604 d.update(k, i);
605 d.end() = i;
606 } else {
607 // move *i to the begin of pts
608 iter_swap(i, d.begin());
609 d.end() = d.begin() + 1;
610 }
611
612 // [d.begin(), d.end()) shall be the range of uncovered points
613 if (i != save_end)
614 while (++i != save_end)
615 if (dist(corner, *i) > d.r) {
616 d.update(k, i);
617 iter_swap(i, d.end());
618 ++d.end();
619 } // if (dist(corner, *i) > d.r)
620
621 // check disjoint for two-pierceability:
622
623 CGAL_optimisation_expensive_assertion(
624 save_end == find_if(d.end(), save_end,
625 [&d, &dist, &corner](const Point_2& p)
626 { return d.r < dist(corner, p); }));
627 CGAL_optimisation_expensive_assertion(
628 d.end() == find_if(d.begin(), d.end(),
629 [&d,&dist, &corner](const Point_2& p)
630 { return d.r >= dist(corner, p); }));
631
632
633 two_cover_points(d, o, ok);
634
635 // restore saved sides of d:
636 d.extreme(k) = save_side1;
637 d.extreme((k+1) % 4) = save_side2;
638
639 if (ok) {
640 // does any rectangle contain the corner?
641 if (d.end() != save_end) {
642 *o++ = corner;
643 d.end() = save_end;
644 }
645 return o;
646 } // if (ok)
647
648 d.end() = save_end;
649 } // for (int k = 0; k < 4; ++k)
650
651 // no piercing set exists:
652 ok = false;
653 return o;
654
655 } // three_cover_points(d, o, ok)
656 } //namespace CGAL
657 namespace CGAL {
658 template < class OutputIterator, class Traits >
659 OutputIterator
four_cover_points(Staircases<Traits> & d,OutputIterator o,bool & ok)660 four_cover_points(Staircases< Traits >& d, OutputIterator o, bool& ok)
661 {
662
663 using std::less;
664 using std::iter_swap;
665 using std::find_if;
666 using std::back_inserter;
667
668 typedef typename Traits::Point_2 Point_2;
669 typedef typename Traits::FT FT;
670 typedef typename Traits::Less_x_2 Less_x_2;
671 typedef typename Traits::Less_y_2 Less_y_2;
672 typedef typename Traits::Signed_x_distance_2 Signed_x_distance_2;
673 typedef typename Traits::Signed_y_distance_2 Signed_y_distance_2;
674 typedef typename Traits::Infinity_distance_2 Infinity_distance_2;
675 typedef typename Staircases< Traits >::Container Container;
676 typedef typename Staircases< Traits >::Iterator Iterator;
677 typedef typename Staircases< Traits >::Citerator Citerator;
678 typedef typename Staircases< Traits >::Intervall Intervall;
679
680 Infinity_distance_2 dist = d.traits.infinity_distance_2_object();
681 Less_x_2 lessx = d.traits.less_x_2_object();
682 Less_y_2 lessy = d.traits.less_y_2_object();
683 Signed_x_distance_2 sdistx = d.traits.signed_x_distance_2_object();
684 Signed_y_distance_2 sdisty = d.traits.signed_y_distance_2_object();
685
686 typename Traits::Construct_point_2_above_right_implicit_point_2
687 cparip =
688 d.traits.construct_point_2_above_right_implicit_point_2_object();
689 typename Traits::Construct_point_2_above_left_implicit_point_2
690 cpalip =
691 d.traits.construct_point_2_above_left_implicit_point_2_object();
692 typename Traits::Construct_point_2_below_right_implicit_point_2
693 cpbrip =
694 d.traits.construct_point_2_below_right_implicit_point_2_object();
695
696
697
698 // test the four corners:
699 for (int j = 0; j < 5; ++j) {
700 const int k = j < 4 ? j : 3;
701
702 // extract all points which are close enough to this point
703 Point_2 corner = d[k];
704 if (j >= 3) {
705 if (j == 3) {
706 Citerator i = d.tlstc_begin();
707 while (sdistx(*i, d.minx) > FT(2) * d.r)
708 ++i;
709 corner = cpbrip(d.minx, *i, d.r);
710 } else {
711 Citerator i = d.tlstc_end();
712 while (sdisty(d.maxy, *--i) > FT(2) * d.r) {}
713 corner = cpbrip(*i, d.maxy, d.r);
714 }
715 }
716
717 // find first point not covered by the rectangle at d[k]
718 Iterator i = find_if(d.begin(), d.end(),
719 [&d,&dist,&corner](const Point_2& p)
720 { return d.r < dist(corner, p); });
721
722 // are all points already covered?
723 if (i == d.end()) {
724 CGAL_optimisation_assertion(k == 0);
725 *o++ = d[0];
726 ok = true;
727 return o;
728 } // if (i == d.end())
729
730 // save changing sides of d:
731 Point_2 save_side1 = d.extreme(k);
732 Point_2 save_side2 = d.extreme((k+1) % 4);
733 Iterator save_end = d.end();
734
735 // now run through it:
736 // initialize the two (possibly) changing sides of d
737 d.extreme(k) = d.extreme((k+1) % 4) = *i;
738
739 // is there any point covered?
740 if (i == d.begin()) {
741 while (++i != d.end() && dist(corner, *i) > d.r)
742 d.update(k, i);
743 d.end() = i;
744 } else {
745 // move *i to the begin of pts
746 iter_swap(i, d.begin());
747 d.end() = d.begin() + 1;
748 }
749
750 // [d.begin(), d.end()) shall be the range of uncovered points
751 if (i != save_end)
752 while (++i != save_end)
753 if (dist(corner, *i) > d.r) {
754 d.update(k, i);
755 iter_swap(i, d.end());
756 ++d.end();
757 } // if (dist(corner, *i) > d.r)
758
759 // check disjoint for two-pierceability:
760
761 CGAL_optimisation_expensive_assertion(
762 save_end == find_if(d.end(), save_end,
763 [&d,&dist,&corner](const Point_2& p)
764 { return d.r < dist(corner, p); }));
765 CGAL_optimisation_expensive_assertion(
766 d.end() == find_if(d.begin(), d.end(),
767 [&d,&dist,&corner](const Point_2& p)
768 { return d.r >= dist(corner, p); }));
769
770
771 three_cover_points(d, o, ok);
772
773 // restore saved sides of d:
774 d.extreme(k) = save_side1;
775 d.extreme((k+1) % 4) = save_side2;
776
777 if (ok) {
778 // does any rectangle contain the corner?
779 if (d.end() != save_end) {
780 *o++ = corner;
781 d.end() = save_end;
782 }
783 return o;
784 } // if (ok)
785
786 d.end() = save_end;
787 } // for (int k = 0; k < 4; ++k)
788
789
790 // test if four covering rectangles can be placed
791 // on the boundary of d, one on each side
792
793 // if there is any point that cannot be covered in this way, stop here
794 if (d.is_middle_empty()) {
795
796 // now try to position the bottom piercing point in each
797 // of the intervalls formed by S_bt and S_br
798 // (no need to consider S_bl, since we move from left
799 // to right and leaving a rectangle won't make piercing easier)
800
801
802 Intervall top_i = d.top_intervall();
803 Intervall left_i = d.left_intervall();
804 Intervall bottom_i = d.bottom_intervall();
805 Intervall right_i = d.right_intervall();
806 Container share;
807 d.shared_intervall(back_inserter(share));
808 Citerator shf = share.end();
809 Citerator shl = share.end();
810
811 Citerator tl = d.tlstc_begin();
812 Citerator lb = d.lbstc_begin();
813 Citerator br = d.brstc_begin();
814 Citerator rt = d.rtstc_begin();
815
816 // make sure the top intervall is covered (left endpoint)
817 // (it might be that top_i.first determines the placement of
818 // the top square)
819 Point_2 top = top_i.first;
820 if (lessx(top, *tl))
821 while (++tl != d.tlstc_end() && !lessx(*tl, top)) {}
822 else
823 top = *tl++;
824
825
826 if (tl != d.tlstc_end()) {
827 for (;;) {
828
829 // make sure the top intervall is covered (right endpoint)
830 if (sdistx(top_i.second, top) > FT(2) * d.r)
831 break;
832
833 // compute position of left square
834 Point_2 left = lessy(left_i.second, *tl) ? *tl : left_i.second;
835
836 // make sure the left intervall is covered
837 if (sdisty(left, left_i.first) <= FT(2) * d.r) {
838
839 // compute position of bottom square
840 while (lb != d.lbstc_end() && sdisty(left, *lb) <= FT(2) * d.r)
841 ++lb;
842 // has the left square reached the bottom-left corner?
843 if (lb == d.lbstc_end())
844 break;
845 Point_2 bottom = lessx(bottom_i.first, *lb) ? bottom_i.first : *lb;
846
847 // check the shared x-intervall
848 if (!share.empty() && d.is_x_greater_y()) {
849 // compute position of top in share
850 #ifndef _MSC_VER
851 while (shf != share.begin() && !lessx(*(shf - 1), top))
852 #else
853 while (shf != Citerator(share.begin()) &&
854 !lessx(*(shf - 1), top))
855 #endif
856 --shf;
857 #ifndef _MSC_VER
858 while (shl != share.begin() &&
859 #else
860 while (shl != Citerator(share.begin()) &&
861 #endif
862 sdistx(*(shl - 1), top) > FT(2) * d.r)
863 --shl;
864
865 // make sure shared intervall is covered (left endpoint)
866 #ifndef _MSC_VER
867 if ((shf != share.begin() || shl == share.end()) &&
868 lessx(share.front(), bottom))
869 bottom = share.front();
870 else if (shl != share.end() && lessx(*shl, bottom))
871 bottom = *shl;
872 #else
873 if ((shf != Citerator(share.begin()) ||
874 shl == Citerator(share.end())) &&
875 lessx(share.front(), bottom))
876 bottom = share.front();
877 else if (shl != Citerator(share.end()) && lessx(*shl, bottom))
878 bottom = *shl;
879 #endif
880
881 }
882
883
884 // make sure the bottom and the shared intervall (right endpoint)
885 // are covered
886 #ifndef _MSC_VER
887 if (sdistx(bottom_i.second, bottom) <= FT(2) * d.r &&
888 (!d.is_x_greater_y() ||
889 ((shl == share.end() ||
890 sdistx(share.back(), bottom) <= FT(2) * d.r) &&
891 (shf == share.begin() ||
892 sdistx(*(shf - 1), bottom) <= FT(2) * d.r))))
893 #else
894 if (sdistx(bottom_i.second, bottom) <= FT(2) * d.r &&
895 ((!d.is_x_greater_y() ||
896 (shl == Citerator(share.end()) ||
897 sdistx(share.back(), bottom) <= FT(2) * d.r) &&
898 (shf == Citerator(share.begin()) ||
899 sdistx(*(shf - 1), bottom) <= FT(2) * d.r))))
900 #endif
901 {
902 // compute position of right square
903 while (br != d.brstc_end() && sdistx(*br, bottom) <= FT(2) * d.r)
904 ++br;
905 // has the bottom square reached the bottom-right corner?
906 if (br == d.brstc_end())
907 break;
908 Point_2 right = lessy(right_i.first, *br) ? right_i.first : *br;
909
910 // check the shared y-intervall
911 if (!share.empty() && !d.is_x_greater_y()) {
912 // compute position of left in share
913 #ifndef _MSC_VER
914 while (shf != share.begin() &&
915 sdisty(left, *(shf - 1)) <= FT(2) * d.r)
916 #else
917 while (shf != Citerator(share.begin()) &&
918 sdisty(left, *(shf - 1)) <= FT(2) * d.r)
919 #endif
920 --shf;
921 #ifndef _MSC_VER
922 while (shl != share.begin() &&
923 #else
924 while (shl != Citerator(share.begin()) &&
925 #endif
926 lessy(left, *(shl - 1)))
927 --shl;
928
929 // make sure shared intervall is covered (bottom endpoint)
930 #ifndef _MSC_VER
931 if ((shf != share.begin() || shl == share.end()) &&
932 lessy(share.front(), right))
933 right = share.front();
934 else if (shl != share.end() && lessy(*shl, right))
935 right = *shl;
936 #else
937 if ((shf != Citerator(share.begin()) ||
938 shl == Citerator(share.end())) &&
939 lessy(share.front(), right))
940 right = share.front();
941 else if (shl != Citerator(share.end()) && lessy(*shl, right))
942 right = *shl;
943 #endif
944
945 }
946
947
948 // make sure the right intervall and the shared intervall
949 // (top endpoint) are covered
950 #ifndef _MSC_VER
951 if (sdisty(right_i.second, right) <= FT(2) * d.r &&
952 (d.is_x_greater_y() ||
953 ((shl == share.end() ||
954 sdisty(share.back(), right) <= FT(2) * d.r) &&
955 (shf == share.begin() ||
956 sdisty(*(shf - 1), right) <= FT(2) * d.r))))
957 #else
958 if (sdisty(right_i.second, right) <= FT(2) * d.r &&
959 (d.is_x_greater_y() ||
960 ((shl == Citerator(share.end()) ||
961 sdisty(share.back(), right) <= FT(2) * d.r) &&
962 (shf == Citerator(share.begin()) ||
963 sdisty(*(shf - 1), right) <= FT(2) * d.r))))
964 #endif
965 {
966 // compute right bound for top square
967 while (rt != d.rtstc_end() &&
968 sdisty(*rt, right) <= FT(2) * d.r)
969 ++rt;
970 // has the right square reached the top-right corner?
971 if (rt == d.rtstc_end())
972 break;
973
974 // Finally: Do we have a covering?
975 if (sdistx(*rt, top) <= FT(2) * d.r) {
976 *o++ = cpbrip(d.minx, left, d.r);
977 *o++ = cparip(bottom, d.miny, d.r);
978 *o++ = cpalip(d.maxx, right, d.r);
979 *o++ = cpbrip(top, d.maxy, d.r);
980 ok = true;
981
982
983 return o;
984 } // if (covering)
985
986 } // if (sdisty(right_i.second, right) <= FT(2) * d.r)
987
988 } // if (bottom and shared intervall are covered)
989
990 } // if (sdisty(left, left_i.first) <= FT(2) * d.r)
991
992 top = *tl;
993 if (!lessy(left_i.second, *tl) || ++tl == d.tlstc_end() ||
994 lessx(bottom_i.first, *lb) || lessy(right_i.first, *br))
995 break;
996
997 } // for (;;)
998 } // if (tl != d.tlstc_end())
999
1000
1001 } // if (!d.is_middle_empty())
1002
1003 ok = false;
1004 return o;
1005
1006 } // four_cover_points(d, o, ok)
1007
1008 struct Two_covering_algorithm {
1009 template < class Traits, class OutputIterator >
1010 OutputIterator
operatorTwo_covering_algorithm1011 operator()(Staircases< Traits >& d,
1012 OutputIterator o,
1013 bool& ok) const
1014 { return two_cover_points(d, o, ok); }
1015 }; // class Two_covering_algorithm
1016 struct Three_covering_algorithm {
1017 template < class Traits, class OutputIterator >
1018 OutputIterator
operatorThree_covering_algorithm1019 operator()(Staircases< Traits >& d,
1020 OutputIterator o,
1021 bool& ok) const
1022 { return three_cover_points(d, o, ok); }
1023 }; // class Three_covering_algorithm
1024 struct Four_covering_algorithm {
1025 template < class Traits, class OutputIterator >
1026 OutputIterator
operatorFour_covering_algorithm1027 operator()(Staircases< Traits >& d,
1028 OutputIterator o,
1029 bool& ok) const
1030 { return four_cover_points(d, o, ok); }
1031 }; // class Four_covering_algorithm
1032
1033 } //namespace CGAL
1034
1035 #if defined(BOOST_MSVC)
1036 # pragma warning(pop)
1037 #endif
1038
1039 #endif // ! (CGAL_PIERCE_RECTANGLES_2_H)
1040