1 /*
2   Copyright 2008 Intel Corporation
3 
4   Use, modification and distribution are subject to the Boost Software License,
5   Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
6   http://www.boost.org/LICENSE_1_0.txt).
7 */
8 #ifndef BOOST_POLYGON_SCAN_ARBITRARY_HPP
9 #define BOOST_POLYGON_SCAN_ARBITRARY_HPP
10 #ifdef BOOST_POLYGON_DEBUG_FILE
11 #include <fstream>
12 #endif
13 #include "polygon_sort_adaptor.hpp"
14 namespace boost { namespace polygon{
15 
16   template <typename Unit>
17   class line_intersection : public scanline_base<Unit> {
18   private:
19     typedef typename scanline_base<Unit>::Point Point;
20 
21     //the first point is the vertex and and second point establishes the slope of an edge eminating from the vertex
22     //typedef std::pair<Point, Point> half_edge;
23     typedef typename scanline_base<Unit>::half_edge half_edge;
24 
25     //scanline comparator functor
26     typedef typename scanline_base<Unit>::less_half_edge less_half_edge;
27     typedef typename scanline_base<Unit>::less_point less_point;
28 
29     //when parallel half edges are encounterd the set of segments is expanded
30     //when a edge leaves the scanline it is removed from the set
31     //when the set is empty the element is removed from the map
32     typedef int segment_id;
33     typedef std::pair<half_edge, std::set<segment_id> > scanline_element;
34     typedef std::map<half_edge, std::set<segment_id>, less_half_edge> edge_scanline;
35     typedef typename edge_scanline::iterator iterator;
36 
37 //     std::map<Unit, std::set<segment_id> > vertical_data_;
38 //     edge_scanline edge_scanline_;
39 //     Unit x_;
40 //     int just_before_;
41 //     segment_id segment_id_;
42 //     std::vector<std::pair<half_edge, int> > event_edges_;
43 //     std::set<Point> intersection_queue_;
44   public:
45 //     inline line_intersection() : vertical_data_(), edge_scanline_(), x_((std::numeric_limits<Unit>::max)()), just_before_(0), segment_id_(0), event_edges_(), intersection_queue_() {
46 //       less_half_edge lessElm(&x_, &just_before_);
47 //       edge_scanline_ = edge_scanline(lessElm);
48 //     }
49 //     inline line_intersection(const line_intersection& that) : vertical_data_(), edge_scanline_(), x_(), just_before_(), segment_id_(), event_edges_(), intersection_queue_() { (*this) = that; }
50 //     inline line_intersection& operator=(const line_intersection& that) {
51 //       x_ = that.x_;
52 //       just_before_ = that.just_before_;
53 //       segment_id_ = that.segment_id_;
54 
55 //       //I cannot simply copy that.edge_scanline_ to this edge_scanline_ becuase the functor store pointers to other members!
56 //       less_half_edge lessElm(&x_, &just_before_);
57 //       edge_scanline_ = edge_scanline(lessElm);
58 
59 //       edge_scanline_.insert(that.edge_scanline_.begin(), that.edge_scanline_.end());
60 //       return *this;
61 //     }
62 
63 //     static inline void between(Point pt, Point pt1, Point pt2) {
64 //       less_point lp;
65 //       if(lp(pt1, pt2))
66 //         return lp(pt, pt2) && lp(pt1, pt);
67 //       return lp(pt, pt1) && lp(pt2, pt);
68 //     }
69 
70     template <typename iT>
compute_histogram_in_y(iT begin,iT end,std::size_t size,std::vector<std::pair<Unit,std::pair<std::size_t,std::size_t>>> & histogram)71     static inline void compute_histogram_in_y(iT begin, iT end, std::size_t size, std::vector<std::pair<Unit, std::pair<std::size_t, std::size_t> > >& histogram) {
72       std::vector<std::pair<Unit, int> > ends;
73       ends.reserve(size * 2);
74       for(iT itr = begin ; itr != end; ++itr) {
75         int count = (*itr).first.first.y() < (*itr).first.second.y() ? 1 : -1;
76         ends.push_back(std::make_pair((*itr).first.first.y(), count));
77         ends.push_back(std::make_pair((*itr).first.second.y(), -count));
78       }
79       polygon_sort(ends.begin(), ends.end());
80       histogram.reserve(ends.size());
81       histogram.push_back(std::make_pair(ends.front().first, std::make_pair(0, 0)));
82       for(typename std::vector<std::pair<Unit, int> >::iterator itr = ends.begin(); itr != ends.end(); ++itr) {
83         if((*itr).first != histogram.back().first) {
84           histogram.push_back(std::make_pair((*itr).first, histogram.back().second));
85         }
86         if((*itr).second < 0)
87           histogram.back().second.second -= (*itr).second;
88         histogram.back().second.first += (*itr).second;
89       }
90     }
91 
92     template <typename iT>
compute_y_cuts(std::vector<Unit> & y_cuts,iT begin,iT end,std::size_t size)93     static inline void compute_y_cuts(std::vector<Unit>& y_cuts, iT begin, iT end, std::size_t size) {
94       if(begin == end) return;
95       if(size < 30) return; //30 is empirically chosen, but the algorithm is not sensitive to this constant
96       std::size_t min_cut = size;
97       iT cut = begin;
98       std::size_t position = 0;
99       std::size_t cut_size = 0;
100       std::size_t histogram_size = std::distance(begin, end);
101       for(iT itr = begin; itr != end; ++itr, ++position) {
102         if(position < histogram_size / 3)
103           continue;
104         if(histogram_size - position < histogram_size / 3) break;
105         if((*itr).second.first < min_cut) {
106           cut = itr;
107           min_cut = (*cut).second.first;
108           cut_size = position;
109         }
110       }
111       if(cut_size == 0 || (*cut).second.first > size / 9) //nine is empirically chosen
112         return;
113       compute_y_cuts(y_cuts, begin, cut, (*cut).second.first + (*cut).second.second);
114       y_cuts.push_back((*cut).first);
115       compute_y_cuts(y_cuts, cut, end, size - (*cut).second.second);
116     }
117 
118     template <typename iT>
validate_scan_divide_and_conquer(std::vector<std::set<Point>> & intersection_points,iT begin,iT end)119     static inline void validate_scan_divide_and_conquer(std::vector<std::set<Point> >& intersection_points,
120                                                         iT begin, iT end) {
121       std::vector<std::pair<Unit, std::pair<std::size_t, std::size_t> > > histogram;
122       compute_histogram_in_y(begin, end, std::distance(begin, end), histogram);
123       std::vector<Unit> y_cuts;
124       compute_y_cuts(y_cuts, histogram.begin(), histogram.end(), std::distance(begin, end));
125       std::map<Unit, std::vector<std::pair<half_edge, segment_id> > > bins;
126       bins[histogram.front().first] = std::vector<std::pair<half_edge, segment_id> >();
127       for(typename std::vector<Unit>::iterator itr = y_cuts.begin(); itr != y_cuts.end(); ++itr) {
128         bins[*itr] = std::vector<std::pair<half_edge, segment_id> >();
129       }
130       for(iT itr = begin; itr != end; ++itr) {
131         typename std::map<Unit, std::vector<std::pair<half_edge, segment_id> > >::iterator lb =
132           bins.lower_bound((std::min)((*itr).first.first.y(), (*itr).first.second.y()));
133         if(lb != bins.begin())
134           --lb;
135         typename std::map<Unit, std::vector<std::pair<half_edge, segment_id> > >::iterator ub =
136           bins.upper_bound((std::max)((*itr).first.first.y(), (*itr).first.second.y()));
137         for( ; lb != ub; ++lb) {
138           (*lb).second.push_back(*itr);
139         }
140       }
141       validate_scan(intersection_points, bins[histogram.front().first].begin(), bins[histogram.front().first].end());
142       for(typename std::vector<Unit>::iterator itr = y_cuts.begin(); itr != y_cuts.end(); ++itr) {
143         validate_scan(intersection_points, bins[*itr].begin(), bins[*itr].end(), *itr);
144       }
145     }
146 
147     template <typename iT>
validate_scan(std::vector<std::set<Point>> & intersection_points,iT begin,iT end)148     static inline void validate_scan(std::vector<std::set<Point> >& intersection_points,
149                                      iT begin, iT end) {
150       validate_scan(intersection_points, begin, end, (std::numeric_limits<Unit>::min)());
151     }
152     //quadratic algorithm to do same work as optimal scan for cross checking
153     template <typename iT>
validate_scan(std::vector<std::set<Point>> & intersection_points,iT begin,iT end,Unit min_y)154     static inline void validate_scan(std::vector<std::set<Point> >& intersection_points,
155                                      iT begin, iT end, Unit min_y) {
156       std::vector<Point> pts;
157       std::vector<std::pair<half_edge, segment_id> > data(begin, end);
158       for(std::size_t i = 0; i < data.size(); ++i) {
159         if(data[i].first.second < data[i].first.first) {
160           std::swap(data[i].first.first, data[i].first.second);
161         }
162       }
163       typename scanline_base<Unit>::compute_intersection_pack pack_;
164       polygon_sort(data.begin(), data.end());
165       //find all intersection points
166       for(typename std::vector<std::pair<half_edge, segment_id> >::iterator outer = data.begin();
167           outer != data.end(); ++outer) {
168         const half_edge& he1 = (*outer).first;
169         //its own end points
170         pts.push_back(he1.first);
171         pts.push_back(he1.second);
172         std::set<Point>& segmentpts = intersection_points[(*outer).second];
173         for(typename std::set<Point>::iterator itr = segmentpts.begin(); itr != segmentpts.end(); ++itr) {
174           if ((*itr).y() >= min_y) {
175             pts.push_back(*itr);
176           }
177         }
178         bool have_first_y = he1.first.y() >= min_y && he1.second.y() >= min_y;
179         for(typename std::vector<std::pair<half_edge, segment_id> >::iterator inner = outer;
180             inner != data.end(); ++inner) {
181           const half_edge& he2 = (*inner).first;
182           if(have_first_y || (he2.first.y() >= min_y && he2.second.y() >= min_y)) {
183             //at least one segment has a low y value within the range
184             if(he1 == he2) continue;
185             if((std::min)(he2. first.get(HORIZONTAL),
186                           he2.second.get(HORIZONTAL)) >=
187                (std::max)(he1.second.get(HORIZONTAL),
188                           he1.first.get(HORIZONTAL)))
189               break;
190             if(he1.first == he2.first || he1.second == he2.second)
191               continue;
192             Point intersection;
193             if(pack_.compute_intersection(intersection, he1, he2)) {
194               //their intersection point
195               pts.push_back(intersection);
196               intersection_points[(*inner).second].insert(intersection);
197               intersection_points[(*outer).second].insert(intersection);
198             }
199           }
200         }
201       }
202       polygon_sort(pts.begin(), pts.end());
203       typename std::vector<Point>::iterator newend = std::unique(pts.begin(), pts.end());
204       typename std::vector<Point>::iterator lfinger = pts.begin();
205       //find all segments that interact with intersection points
206       for(typename std::vector<std::pair<half_edge, segment_id> >::iterator outer = data.begin();
207           outer != data.end(); ++outer) {
208         const half_edge& he1 = (*outer).first;
209         segment_id id1 = (*outer).second;
210         //typedef rectangle_data<Unit> Rectangle;
211         //Rectangle rect1;
212         //set_points(rect1, he1.first, he1.second);
213         //typename std::vector<Point>::iterator itr = lower_bound(pts.begin(), newend, (std::min)(he1.first, he1.second));
214         //typename std::vector<Point>::iterator itr2 = upper_bound(pts.begin(), newend, (std::max)(he1.first, he1.second));
215         Point startpt = (std::min)(he1.first, he1.second);
216         Point stoppt = (std::max)(he1.first, he1.second);
217         //while(itr != newend && itr != pts.begin() && (*itr).get(HORIZONTAL) >= (std::min)(he1.first.get(HORIZONTAL), he1.second.get(HORIZONTAL))) --itr;
218         //while(itr2 != newend && (*itr2).get(HORIZONTAL) <= (std::max)(he1.first.get(HORIZONTAL), he1.second.get(HORIZONTAL))) ++itr2;
219         //itr = pts.begin();
220         //itr2 = pts.end();
221         while(lfinger != newend && (*lfinger).x() < startpt.x()) ++lfinger;
222         for(typename std::vector<Point>::iterator itr = lfinger ; itr != newend && (*itr).x() <= stoppt.x(); ++itr) {
223           if(scanline_base<Unit>::intersects_grid(*itr, he1))
224             intersection_points[id1].insert(*itr);
225         }
226       }
227     }
228 
229     template <typename iT, typename property_type>
validate_scan(std::vector<std::pair<half_edge,std::pair<property_type,int>>> & output_segments,iT begin,iT end)230     static inline void validate_scan(std::vector<std::pair<half_edge, std::pair<property_type, int> > >& output_segments,
231                                      iT begin, iT end) {
232       std::vector<std::pair<property_type, int> > input_properties;
233       std::vector<std::pair<half_edge, int> > input_segments, intermediate_segments;
234       int index = 0;
235       for( ; begin != end; ++begin) {
236         input_properties.push_back((*begin).second);
237         input_segments.push_back(std::make_pair((*begin).first, index++));
238       }
239       validate_scan(intermediate_segments, input_segments.begin(), input_segments.end());
240       for(std::size_t i = 0; i < intermediate_segments.size(); ++i) {
241         output_segments.push_back(std::make_pair(intermediate_segments[i].first,
242                                                  input_properties[intermediate_segments[i].second]));
243         less_point lp;
244         if(lp(output_segments.back().first.first, output_segments.back().first.second) !=
245            lp(input_segments[intermediate_segments[i].second].first.first,
246               input_segments[intermediate_segments[i].second].first.second)) {
247           //edge changed orientation, invert count on edge
248           output_segments.back().second.second *= -1;
249         }
250         if(!scanline_base<Unit>::is_vertical(input_segments[intermediate_segments[i].second].first) &&
251            scanline_base<Unit>::is_vertical(output_segments.back().first)) {
252           output_segments.back().second.second *= -1;
253         }
254         if(lp(output_segments.back().first.second, output_segments.back().first.first)) {
255           std::swap(output_segments.back().first.first, output_segments.back().first.second);
256         }
257       }
258     }
259 
260     template <typename iT>
validate_scan(std::vector<std::pair<half_edge,int>> & output_segments,iT begin,iT end)261     static inline void validate_scan(std::vector<std::pair<half_edge, int> >& output_segments,
262                                      iT begin, iT end) {
263       std::vector<std::set<Point> > intersection_points(std::distance(begin, end));
264       validate_scan_divide_and_conquer(intersection_points, begin, end);
265       //validate_scan(intersection_points, begin, end);
266       segment_intersections(output_segments, intersection_points, begin, end);
267 //       std::pair<segment_id, segment_id> offenders;
268 //       if(!verify_scan(offenders, output_segments.begin(), output_segments.end())) {
269 //         std::cout << "break here!\n";
270 //         for(typename std::set<Point>::iterator itr = intersection_points[offenders.first].begin();
271 //             itr != intersection_points[offenders.first].end(); ++itr) {
272 //           std::cout << (*itr).x() << " " << (*itr).y() << " ";
273 //         } std::cout << "\n";
274 //         for(typename std::set<Point>::iterator itr = intersection_points[offenders.second].begin();
275 //             itr != intersection_points[offenders.second].end(); ++itr) {
276 //           std::cout << (*itr).x() << " " << (*itr).y() << " ";
277 //         } std::cout << "\n";
278 //         exit(1);
279 //       }
280     }
281 
282     //quadratic algorithm to find intersections
283     template <typename iT, typename segment_id>
verify_scan(std::pair<segment_id,segment_id> & offenders,iT begin,iT end)284     static inline bool verify_scan(std::pair<segment_id, segment_id>& offenders,
285                                    iT begin, iT end) {
286 
287       std::vector<std::pair<half_edge, segment_id> > data(begin, end);
288       for(std::size_t i = 0; i < data.size(); ++i) {
289         if(data[i].first.second < data[i].first.first) {
290           std::swap(data[i].first.first, data[i].first.second);
291         }
292       }
293       polygon_sort(data.begin(), data.end());
294       for(typename std::vector<std::pair<half_edge, segment_id> >::iterator outer = data.begin();
295           outer != data.end(); ++outer) {
296         const half_edge& he1 = (*outer).first;
297         segment_id id1 = (*outer).second;
298         for(typename std::vector<std::pair<half_edge, segment_id> >::iterator inner = outer;
299             inner != data.end(); ++inner) {
300           const half_edge& he2 = (*inner).first;
301           if(he1 == he2) continue;
302           if((std::min)(he2. first.get(HORIZONTAL),
303                         he2.second.get(HORIZONTAL)) >
304              (std::max)(he1.second.get(HORIZONTAL),
305                         he1.first.get(HORIZONTAL)))
306             break;
307           segment_id id2 = (*inner).second;
308           if(scanline_base<Unit>::intersects(he1, he2)) {
309             offenders.first = id1;
310             offenders.second = id2;
311             //std::cout << he1.first.x() << " " << he1.first.y() << " " << he1.second.x() << " " << he1.second.y() << " " << he2.first.x() << " " << he2.first.y() << " " << he2.second.x() << " " << he2.second.y() << "\n";
312             return false;
313           }
314         }
315       }
316       return true;
317     }
318 
319     class less_point_down_slope : public std::binary_function<Point, Point, bool> {
320     public:
less_point_down_slope()321       inline less_point_down_slope() {}
operator ()(const Point & pt1,const Point & pt2) const322       inline bool operator () (const Point& pt1, const Point& pt2) const {
323         if(pt1.get(HORIZONTAL) < pt2.get(HORIZONTAL)) return true;
324         if(pt1.get(HORIZONTAL) == pt2.get(HORIZONTAL)) {
325           if(pt1.get(VERTICAL) > pt2.get(VERTICAL)) return true;
326         }
327         return false;
328       }
329     };
330 
331     template <typename iT>
segment_edge(std::vector<std::pair<half_edge,int>> & output_segments,const half_edge &,segment_id id,iT begin,iT end)332     static inline void segment_edge(std::vector<std::pair<half_edge, int> >& output_segments,
333                                     const half_edge& , segment_id id, iT begin, iT end) {
334       iT current = begin;
335       iT next = begin;
336       ++next;
337       while(next != end) {
338         output_segments.push_back(std::make_pair(half_edge(*current, *next), id));
339         current = next;
340         ++next;
341       }
342     }
343 
344     template <typename iT>
segment_intersections(std::vector<std::pair<half_edge,int>> & output_segments,std::vector<std::set<Point>> & intersection_points,iT begin,iT end)345     static inline void segment_intersections(std::vector<std::pair<half_edge, int> >& output_segments,
346                                              std::vector<std::set<Point> >& intersection_points,
347                                              iT begin, iT end) {
348       for(iT iter = begin; iter != end; ++iter) {
349         //less_point lp;
350         const half_edge& he = (*iter).first;
351         //if(lp(he.first, he.second)) {
352         //  //it is the begin event
353           segment_id id = (*iter).second;
354           const std::set<Point>& pts = intersection_points[id];
355           Point hpt(he.first.get(HORIZONTAL)+1, he.first.get(VERTICAL));
356           if(!scanline_base<Unit>::is_vertical(he) && scanline_base<Unit>::less_slope(he.first.get(HORIZONTAL), he.first.get(VERTICAL),
357                                             he.second, hpt)) {
358             //slope is below horizontal
359             std::vector<Point> tmpPts;
360             tmpPts.reserve(pts.size());
361             tmpPts.insert(tmpPts.end(), pts.begin(), pts.end());
362             less_point_down_slope lpds;
363             polygon_sort(tmpPts.begin(), tmpPts.end(), lpds);
364             segment_edge(output_segments, he, id, tmpPts.begin(), tmpPts.end());
365           } else {
366             segment_edge(output_segments, he, id, pts.begin(), pts.end());
367           }
368           //}
369       }
370     }
371 
372 //     //iT iterator over unsorted pair<Point> representing line segments of input
373 //     //output_segments is populated with fully intersected output line segment half
374 //     //edges and the index of the input segment that they are assoicated with
375 //     //duplicate output half edges with different ids will be generated in the case
376 //     //that parallel input segments intersection
377 //     //outputs are in sorted order and include both begin and end events for
378 //     //each segment
379 //     template <typename iT>
380 //     inline void scan(std::vector<std::pair<half_edge, int> >& output_segments,
381 //                      iT begin, iT end) {
382 //       std::map<segment_id, std::set<Point> > intersection_points;
383 //       scan(intersection_points, begin, end);
384 //       segment_intersections(output_segments, intersection_points, begin, end);
385 //     }
386 
387 //     //iT iterator over sorted sequence of half edge, segment id pairs representing segment begin and end points
388 //     //intersection points provides a mapping from input segment id (vector index) to the set
389 //     //of intersection points assocated with that input segment
390 //     template <typename iT>
391 //     inline void scan(std::map<segment_id, std::set<Point> >& intersection_points,
392 //                      iT begin, iT end) {
393 //       for(iT iter = begin; iter != end; ++iter) {
394 //         const std::pair<half_edge, int>& elem = *iter;
395 //         const half_edge& he = elem.first;
396 //         Unit current_x = he.first.get(HORIZONTAL);
397 //         if(current_x != x_) {
398 //           process_scan_event(intersection_points);
399 //           while(!intersection_queue_.empty() &&
400 //                 (*(intersection_queue_.begin()).get(HORIZONTAL) < current_x)) {
401 //             x_ = *(intersection_queue_.begin()).get(HORIZONTAL);
402 //             process_intersections_at_scan_event(intersection_points);
403 //           }
404 //           x_ = current_x;
405 //         }
406 //         event_edges_.push_back(elem);
407 //       }
408 //       process_scan_event(intersection_points);
409 //     }
410 
411 //     inline iterator lookup(const half_edge& he) {
412 //       return edge_scanline_.find(he);
413 //     }
414 
415 //     inline void insert_into_scanline(const half_edge& he, int id) {
416 //       edge_scanline_[he].insert(id);
417 //     }
418 
419 //     inline void lookup_and_remove(const half_edge& he, int id) {
420 //       iterator remove_iter = lookup(he);
421 //       if(remove_iter == edge_scanline_.end()) {
422 //         //std::cout << "failed to find removal segment in scanline\n";
423 //         return;
424 //       }
425 //       std::set<segment_id>& ids = (*remove_iter).second;
426 //       std::set<segment_id>::iterator id_iter = ids.find(id);
427 //       if(id_iter == ids.end()) {
428 //         //std::cout << "failed to find removal segment id in scanline set\n";
429 //         return;
430 //       }
431 //       ids.erase(id_iter);
432 //       if(ids.empty())
433 //         edge_scanline_.erase(remove_iter);
434 //     }
435 
436 //     static inline void update_segments(std::map<segment_id, std::set<Point> >& intersection_points,
437 //                                        const std::set<segment_id>& segments, Point pt) {
438 //       for(std::set<segment_id>::const_iterator itr = segments.begin(); itr != segments.end(); ++itr) {
439 //         intersection_points[*itr].insert(pt);
440 //       }
441 //     }
442 
443 //     inline void process_intersections_at_scan_event(std::map<segment_id, std::set<Point> >& intersection_points) {
444 //       //there may be additional intersection points at this x location that haven't been
445 //       //found yet if vertical or near vertical line segments intersect more than
446 //       //once before the next x location
447 //       just_before_ = true;
448 //       std::set<iterator> intersecting_elements;
449 //       std::set<Unit> intersection_locations;
450 //       typedef typename std::set<Point>::iterator intersection_iterator;
451 //       intersection_iterator iter;
452 //       //first find all secondary intersection locations and all scanline iterators
453 //       //that are intersecting
454 //       for(iter = intersection_queue_.begin();
455 //           iter != intersection_queue_.end() && (*iter).get(HORIZONTAL) == x_; ++iter) {
456 //         Point pt = *iter;
457 //         Unit y = pt.get(VERTICAL);
458 //         intersection_locations.insert(y);
459 //         //if x_ is max there can be only end events and no sloping edges
460 //         if(x_ != (std::numeric_limits<Unit>::max)()) {
461 //           //deal with edges that project to the right of scanline
462 //           //first find the edges in the scanline adjacent to primary intersectin points
463 //           //lookup segment in scanline at pt
464 //           iterator itr = edge_scanline_.lower_bound(half_edge(pt, Point(x_+1, y)));
465 //           //look above pt in scanline until reaching end or segment that doesn't intersect
466 //           //1x1 grid upper right of pt
467 //           //look below pt in scanline until reaching begin or segment that doesn't interset
468 //           //1x1 grid upper right of pt
469 
470 //           //second find edges in scanline on the y interval of each edge found in the previous
471 //           //step for x_ to x_ + 1
472 
473 //           //third find overlaps in the y intervals of all found edges to find all
474 //           //secondary intersection points
475 
476 //         }
477 //       }
478 //       //erase the intersection points from the queue
479 //       intersection_queue_.erase(intersection_queue_.begin(), iter);
480 //       std::vector<scanline_element> insertion_edges;
481 //       insertion_edges.reserve(intersecting_elements.size());
482 //       std::vector<std::pair<Unit, iterator> > sloping_ends;
483 //       //do all the work of updating the output of all intersecting
484 //       for(typename std::set<iterator>::iterator inter_iter = intersecting_elements.begin();
485 //           inter_iter != intersecting_elements.end(); ++inter_iter) {
486 //         //if it is horizontal update it now and continue
487 //         if(is_horizontal((*inter_iter).first)) {
488 //           update_segments(intersection_points, (*inter_iter).second, Point(x_, (*inter_iter).first.get(VERTICAL)));
489 //         } else {
490 //           //if x_ is max there can be only end events and no sloping edges
491 //           if(x_ != (std::numeric_limits<Unit>::max)()) {
492 //             //insert its end points into the vector of sloping ends
493 //             const half_edge& he = (*inter_iter).first;
494 //             Unit y = evalAtXforY(x_, he.first, he.second);
495 //             Unit y2 = evalAtXforY(x_+1, he.first, he.second);
496 //             if(y2 >= y) y2 +=1; //we round up, in exact case we don't worry about overbite of one
497 //             else y += 1; //downward sloping round up
498 //             sloping_ends.push_back(std::make_pair(y, inter_iter));
499 //             sloping_ends.push_back(std::make_pair(y2, inter_iter));
500 //           }
501 //         }
502 //       }
503 
504 //       //merge sloping element data
505 //       polygon_sort(sloping_ends.begin(), sloping_ends.end());
506 //       std::map<Unit, std::set<iterator> > sloping_elements;
507 //       std::set<iterator> merge_elements;
508 //       for(typename std::vector<std::pair<Unit, iterator> >::iterator slop_iter = sloping_ends.begin();
509 //           slop_iter == sloping_ends.end(); ++slop_iter) {
510 //         //merge into sloping elements
511 //         typename std::set<iterator>::iterator merge_iterator = merge_elements.find((*slop_iter).second);
512 //         if(merge_iterator == merge_elements.end()) {
513 //           merge_elements.insert((*slop_iter).second);
514 //         } else {
515 //           merge_elements.erase(merge_iterator);
516 //         }
517 //         sloping_elements[(*slop_iter).first] = merge_elements;
518 //       }
519 
520 //       //scan intersection points
521 //       typename std::map<Unit, std::set<segment_id> >::iterator vertical_iter = vertical_data_.begin();
522 //       typename std::map<Unit, std::set<iterator> >::iterator sloping_iter = sloping_elements.begin();
523 //       for(typename std::set<Unit>::iterator position_iter = intersection_locations.begin();
524 //           position_iter == intersection_locations.end(); ++position_iter) {
525 //         //look for vertical segments that intersect this point and update them
526 //         Unit y = *position_iter;
527 //         Point pt(x_, y);
528 //         //handle vertical segments
529 //         if(vertical_iter != vertical_data_.end()) {
530 //           typename std::map<Unit, std::set<segment_id> >::iterator next_vertical = vertical_iter;
531 //           for(++next_vertical; next_vertical != vertical_data_.end() &&
532 //                 (*next_vertical).first < y; ++next_vertical) {
533 //             vertical_iter = next_vertical;
534 //           }
535 //           if((*vertical_iter).first < y && !(*vertical_iter).second.empty()) {
536 //             update_segments(intersection_points, (*vertical_iter).second, pt);
537 //             ++vertical_iter;
538 //             if(vertical_iter != vertical_data_.end() && (*vertical_iter).first == y)
539 //               update_segments(intersection_points, (*vertical_iter).second, pt);
540 //           }
541 //         }
542 //         //handle sloping segments
543 //         if(sloping_iter != sloping_elements.end()) {
544 //           typename std::map<Unit, std::set<iterator> >::iterator next_sloping = sloping_iter;
545 //           for(++next_sloping; next_sloping != sloping_elements.end() &&
546 //                 (*next_sloping).first < y; ++next_sloping) {
547 //             sloping_iter = next_sloping;
548 //           }
549 //           if((*sloping_iter).first < y && !(*sloping_iter).second.empty()) {
550 //             for(typename std::set<iterator>::iterator element_iter = (*sloping_iter).second.begin();
551 //                 element_iter != (*sloping_iter).second.end(); ++element_iter) {
552 //               const half_edge& he = (*element_iter).first;
553 //               if(intersects_grid(pt, he)) {
554 //                 update_segments(intersection_points, (*element_iter).second, pt);
555 //               }
556 //             }
557 //             ++sloping_iter;
558 //             if(sloping_iter != sloping_elements.end() && (*sloping_iter).first == y &&
559 //                !(*sloping_iter).second.empty()) {
560 //               for(typename std::set<iterator>::iterator element_iter = (*sloping_iter).second.begin();
561 //                   element_iter != (*sloping_iter).second.end(); ++element_iter) {
562 //                 const half_edge& he = (*element_iter).first;
563 //                 if(intersects_grid(pt, he)) {
564 //                   update_segments(intersection_points, (*element_iter).second, pt);
565 //                 }
566 //               }
567 //             }
568 //           }
569 //         }
570 //       }
571 
572 //       //erase and reinsert edges into scanline with check for future intersection
573 //     }
574 
575 //     inline void process_scan_event(std::map<segment_id, std::set<Point> >& intersection_points) {
576 //       just_before_ = true;
577 
578 //       //process end events by removing those segments from the scanline
579 //       //and insert vertices of all events into intersection queue
580 //       Point prev_point((std::numeric_limits<Unit>::min)(), (std::numeric_limits<Unit>::min)());
581 //       less_point lp;
582 //       std::set<segment_id> vertical_ids;
583 //       vertical_data_.clear();
584 //       for(std::size_t i = 0; i < event_edges_.size(); ++i) {
585 //         segment_id id = event_edges_[i].second;
586 //         const half_edge& he = event_edges_[i].first;
587 //         //vertical half edges are handled during intersection processing because
588 //         //they cannot be inserted into the scanline
589 //         if(!is_vertical(he)) {
590 //           if(lp(he.second, he.first)) {
591 //             //half edge is end event
592 //             lookup_and_remove(he, id);
593 //           } else {
594 //             //half edge is begin event
595 //             insert_into_scanline(he, id);
596 //             //note that they will be immediately removed and reinserted after
597 //             //handling their intersection (vertex)
598 //             //an optimization would allow them to be processed specially to avoid the redundant
599 //             //removal and reinsertion
600 //           }
601 //         } else {
602 //           //common case if you are lucky
603 //           //update the map of y to set of segment id
604 //           if(lp(he.second, he.first)) {
605 //             //half edge is end event
606 //             std::set<segment_id>::iterator itr = vertical_ids.find(id);
607 //             if(itr == vertical_ids.end()) {
608 //               //std::cout << "Failed to find end event id in vertical ids\n";
609 //             } else {
610 //               vertical_ids.erase(itr);
611 //               vertical_data_[he.first.get(HORIZONTAL)] = vertical_ids;
612 //             }
613 //           } else {
614 //             //half edge is a begin event
615 //             vertical_ids.insert(id);
616 //             vertical_data_[he.first.get(HORIZONTAL)] = vertical_ids;
617 //           }
618 //         }
619 //         //prevent repeated insertion of same vertex into intersection queue
620 //         if(prev_point != he.first)
621 //           intersection_queue_.insert(he.first);
622 //         else
623 //           prev_point = he.first;
624 //         // process intersections at scan event
625 //         process_intersections_at_scan_event(intersection_points);
626 //       }
627 //       event_edges_.clear();
628 //     }
629 
630   public:
631     template <typename stream_type>
test_validate_scan(stream_type & stdcout)632     static inline bool test_validate_scan(stream_type& stdcout) {
633       std::vector<std::pair<half_edge, segment_id> > input, edges;
634       input.push_back(std::make_pair(half_edge(Point(0, 0), Point(0, 10)), 0));
635       input.push_back(std::make_pair(half_edge(Point(0, 0), Point(10, 10)), 1));
636       std::pair<segment_id, segment_id> result;
637       validate_scan(edges, input.begin(), input.end());
638       if(!verify_scan(result, edges.begin(), edges.end())) {
639         stdcout << "s fail1 " << result.first << " " << result.second << "\n";
640         return false;
641       }
642       input.push_back(std::make_pair(half_edge(Point(0, 5), Point(5, 5)), 2));
643       edges.clear();
644       validate_scan(edges, input.begin(), input.end());
645       if(!verify_scan(result, edges.begin(), edges.end())) {
646         stdcout << "s fail2 " << result.first << " " << result.second << "\n";
647         return false;
648       }
649       input.pop_back();
650       input.push_back(std::make_pair(half_edge(Point(1, 0), Point(11, 11)), input.size()));
651       edges.clear();
652       validate_scan(edges, input.begin(), input.end());
653       if(!verify_scan(result, edges.begin(), edges.end())) {
654         stdcout << "s fail3 " << result.first << " " << result.second << "\n";
655         return false;
656       }
657       input.push_back(std::make_pair(half_edge(Point(1, 0), Point(10, 11)), input.size()));
658       edges.clear();
659       validate_scan(edges, input.begin(), input.end());
660       if(!verify_scan(result, edges.begin(), edges.end())) {
661         stdcout << "s fail4 " << result.first << " " << result.second << "\n";
662         return false;
663       }
664       input.pop_back();
665       input.push_back(std::make_pair(half_edge(Point(1, 2), Point(11, 11)), input.size()));
666       edges.clear();
667       validate_scan(edges, input.begin(), input.end());
668       if(!verify_scan(result, edges.begin(), edges.end())) {
669         stdcout << "s fail5 " << result.first << " " << result.second << "\n";
670         return false;
671       }
672       input.push_back(std::make_pair(half_edge(Point(0, 5), Point(0, 11)), input.size()));
673       edges.clear();
674       validate_scan(edges, input.begin(), input.end());
675       if(!verify_scan(result, edges.begin(), edges.end())) {
676         stdcout << "s fail6 " << result.first << " " << result.second << "\n";
677         return false;
678       }
679       input.pop_back();
680       for(std::size_t i = 0; i < input.size(); ++i) {
681         std::swap(input[i].first.first, input[i].first.second);
682       }
683       edges.clear();
684       validate_scan(edges, input.begin(), input.end());
685       if(!verify_scan(result, edges.begin(), edges.end())) {
686         stdcout << "s fail5 2 " << result.first << " " << result.second << "\n";
687         return false;
688       }
689       for(std::size_t i = 0; i < input.size(); ++i) {
690         input[i].first.first = Point(input[i].first.first.get(HORIZONTAL) * -1,
691                                      input[i].first.first.get(VERTICAL) * -1);
692         input[i].first.second = Point(input[i].first.second.get(HORIZONTAL) * -1,
693                                      input[i].first.second.get(VERTICAL) * -1);
694       }
695       edges.clear();
696       validate_scan(edges, input.begin(), input.end());
697       stdcout << edges.size() << "\n";
698       if(!verify_scan(result, edges.begin(), edges.end())) {
699         stdcout << "s fail5 3 " << result.first << " " << result.second << "\n";
700         return false;
701       }
702       input.clear();
703       edges.clear();
704       input.push_back(std::make_pair(half_edge(Point(5, 7), Point(7, 6)), 0));
705       input.push_back(std::make_pair(half_edge(Point(2, 4), Point(6, 7)), 1));
706             validate_scan(edges, input.begin(), input.end());
707       if(!verify_scan(result, edges.begin(), edges.end())) {
708         stdcout << "s fail2 1 " << result.first << " " << result.second << "\n";
709         print(input);
710         print(edges);
711         return false;
712       }
713       input.clear();
714       edges.clear();
715       input.push_back(std::make_pair(half_edge(Point(3, 2), Point(1, 7)), 0));
716       input.push_back(std::make_pair(half_edge(Point(0, 6), Point(7, 4)), 1));
717             validate_scan(edges, input.begin(), input.end());
718       if(!verify_scan(result, edges.begin(), edges.end())) {
719         stdcout << "s fail2 2 " << result.first << " " << result.second << "\n";
720         print(input);
721         print(edges);
722         return false;
723       }
724       input.clear();
725       edges.clear();
726       input.push_back(std::make_pair(half_edge(Point(6, 6), Point(1, 0)), 0));
727       input.push_back(std::make_pair(half_edge(Point(3, 6), Point(2, 3)), 1));
728             validate_scan(edges, input.begin(), input.end());
729       if(!verify_scan(result, edges.begin(), edges.end())) {
730         stdcout << "s fail2 3 " << result.first << " " << result.second << "\n";
731         print(input);
732         print(edges);
733         return false;
734       }
735       input.clear();
736       edges.clear();
737       input.push_back(std::make_pair(half_edge(Point(0, 0), Point(7, 0)), 0));
738       input.push_back(std::make_pair(half_edge(Point(6, 0), Point(2, 0)), 1));
739             validate_scan(edges, input.begin(), input.end());
740       if(!verify_scan(result, edges.begin(), edges.end())) {
741         stdcout << "s fail2 4 " << result.first << " " << result.second << "\n";
742         print(input);
743         print(edges);
744         return false;
745       }
746       input.clear();
747       edges.clear();
748       input.push_back(std::make_pair(half_edge(Point(-17333131 - -17208131, -10316869 - -10191869), Point(0, 0)), 0));
749       input.push_back(std::make_pair(half_edge(Point(-17291260 - -17208131, -10200000 - -10191869), Point(-17075000 - -17208131, -10200000 - -10191869)), 1));
750       validate_scan(edges, input.begin(), input.end());
751       if(!verify_scan(result, edges.begin(), edges.end())) {
752         stdcout << "s fail2 5 " << result.first << " " << result.second << "\n";
753         print(input);
754         print(edges);
755         return false;
756       }
757       input.clear();
758       edges.clear();
759       input.push_back(std::make_pair(half_edge(Point(-17333131, -10316869), Point(-17208131, -10191869)), 0));
760       input.push_back(std::make_pair(half_edge(Point(-17291260, -10200000), Point(-17075000, -10200000)), 1));
761       validate_scan(edges, input.begin(), input.end());
762       if(!verify_scan(result, edges.begin(), edges.end())) {
763         stdcout << "s fail2 6 " << result.first << " " << result.second << "\n";
764         print(input);
765         print(edges);
766         return false;
767       }
768       input.clear();
769       edges.clear();
770       input.push_back(std::make_pair(half_edge(Point(-9850009+9853379, -286971+290340), Point(-12777869+9853379, -3214831+290340)), 0));
771       input.push_back(std::make_pair(half_edge(Point(-5223510+9853379, -290340+290340), Point(-9858140+9853379, -290340+290340)), 1));
772       validate_scan(edges, input.begin(), input.end());
773       print(edges);
774       if(!verify_scan(result, edges.begin(), edges.end())) {
775         stdcout << "s fail2 7 " << result.first << " " << result.second << "\n";
776         print(input);
777         print(edges);
778         return false;
779       }
780       input.clear();
781       edges.clear();
782       input.push_back(std::make_pair(half_edge(Point(-9850009, -286971), Point(-12777869, -3214831)), 0));
783       input.push_back(std::make_pair(half_edge(Point(-5223510, -290340), Point(-9858140, -290340)), 1));
784       validate_scan(edges, input.begin(), input.end());
785       if(!verify_scan(result, edges.begin(), edges.end())) {
786         stdcout << "s fail2 8 " << result.first << " " << result.second << "\n";
787         print(input);
788         print(edges);
789         return false;
790       }
791       //3 3 2 2: 0; 4 2 0 6: 1; 0 3 6 3: 2; 4 1 5 5: 3;
792       input.clear();
793       edges.clear();
794       input.push_back(std::make_pair(half_edge(Point(3, 3), Point(2, 2)), 0));
795       input.push_back(std::make_pair(half_edge(Point(4, 2), Point(0, 6)), 1));
796       input.push_back(std::make_pair(half_edge(Point(0, 3), Point(6, 3)), 2));
797       input.push_back(std::make_pair(half_edge(Point(4, 1), Point(5, 5)), 3));
798             validate_scan(edges, input.begin(), input.end());
799       if(!verify_scan(result, edges.begin(), edges.end())) {
800         stdcout << "s fail4 1 " << result.first << " " << result.second << "\n";
801         print(input);
802         print(edges);
803         return false;
804       }
805       //5 7 1 3: 0; 4 5 2 1: 1; 2 5 2 1: 2; 4 1 5 3: 3;
806       input.clear();
807       edges.clear();
808       input.push_back(std::make_pair(half_edge(Point(5, 7), Point(1, 3)), 0));
809       input.push_back(std::make_pair(half_edge(Point(4, 5), Point(2, 1)), 1));
810       input.push_back(std::make_pair(half_edge(Point(2, 5), Point(2, 1)), 2));
811       input.push_back(std::make_pair(half_edge(Point(4, 1), Point(5, 3)), 3));
812             validate_scan(edges, input.begin(), input.end());
813       if(!verify_scan(result, edges.begin(), edges.end())) {
814         stdcout << "s fail4 2 " << result.first << " " << result.second << "\n";
815         print(input);
816         print(edges);
817         return false;
818       }
819       //1 0 -4 -1: 0; 0 0 2 -1: 1;
820       input.clear();
821       edges.clear();
822       input.push_back(std::make_pair(half_edge(Point(1, 0), Point(-4, -1)), 0));
823       input.push_back(std::make_pair(half_edge(Point(0, 0), Point(2, -1)), 1));
824             validate_scan(edges, input.begin(), input.end());
825       if(!verify_scan(result, edges.begin(), edges.end())) {
826         stdcout << "s fail2 5 " << result.first << " " << result.second << "\n";
827         print(input);
828         print(edges);
829         return false;
830       }
831       Unit min_c =0;
832       Unit max_c =0;
833       for(unsigned int outer = 0; outer < 1000; ++outer) {
834         input.clear();
835         for(unsigned int i = 0; i < 4; ++i) {
836           Unit x1 = rand();
837           Unit x2 = rand();
838           Unit y1 = rand();
839           Unit y2 = rand();
840           int neg1 = rand() % 2;
841           if(neg1) x1 *= -1;
842           int neg2 = rand() % 2;
843           if(neg2) x2 *= -1;
844           int neg3 = rand() % 2;
845           if(neg3) y1 *= -1;
846           int neg4 = rand() % 2;
847           if(neg4) y2 *= -1;
848           if(x1 < min_c) min_c = x1;
849           if(x2 < min_c) min_c = x2;
850           if(y1 < min_c) min_c = y1;
851           if(y2 < min_c) min_c = y2;
852           if(x1 > max_c) max_c = x1;
853           if(x2 > max_c) max_c = x2;
854           if(y1 > max_c) max_c = y1;
855           if(y2 > max_c) max_c = y2;
856           Point pt1(x1, y1);
857           Point pt2(x2, y2);
858           if(pt1 != pt2)
859             input.push_back(std::make_pair(half_edge(pt1, pt2), i));
860         }
861         edges.clear();
862         validate_scan(edges, input.begin(), input.end());
863         if(!verify_scan(result, edges.begin(), edges.end())) {
864           stdcout << "s fail9 " << outer << ": " << result.first << " " << result.second << "\n";
865           print(input);
866           print(edges);
867           return false;
868         }
869       }
870       return true;
871     }
872 
873     //static void print(const std::pair<half_edge, segment_id>& segment) {
874       //std::cout << segment.first.first << " " << segment.first.second << ": " << segment.second << "; ";
875     //}
print(const std::vector<std::pair<half_edge,segment_id>> & vec)876     static void print(const std::vector<std::pair<half_edge, segment_id> >& vec) {
877       for(std::size_t i = 0; i < vec.size(); ++ i) {
878       //  print(vec[i]);
879       }
880       //std::cout << "\n";
881     }
882 
883     template <typename stream_type>
test_verify_scan(stream_type & stdcout)884     static inline bool test_verify_scan(stream_type& stdcout) {
885       std::vector<std::pair<half_edge, segment_id> > edges;
886       edges.push_back(std::make_pair(half_edge(Point(0, 0), Point(0, 10)), 0));
887       edges.push_back(std::make_pair(half_edge(Point(0, 0), Point(10, 10)), 1));
888       std::pair<segment_id, segment_id> result;
889       if(!verify_scan(result, edges.begin(), edges.end())) {
890         stdcout << "fail1\n";
891         return false;
892       }
893       edges.push_back(std::make_pair(half_edge(Point(0, 5), Point(5, 5)), 2));
894       if(verify_scan(result, edges.begin(), edges.end())) {
895         stdcout << "fail2\n";
896         return false;
897       }
898       edges.pop_back();
899       edges.push_back(std::make_pair(half_edge(Point(1, 0), Point(11, 11)), (segment_id)edges.size()));
900       if(!verify_scan(result, edges.begin(), edges.end())) {
901         stdcout << "fail3\n";
902         return false;
903       }
904       edges.push_back(std::make_pair(half_edge(Point(1, 0), Point(10, 11)), (segment_id)edges.size()));
905       if(verify_scan(result, edges.begin(), edges.end())) {
906         stdcout << "fail4\n";
907         return false;
908       }
909       edges.pop_back();
910       edges.push_back(std::make_pair(half_edge(Point(1, 2), Point(11, 11)), (segment_id)edges.size()));
911       if(!verify_scan(result, edges.begin(), edges.end())) {
912         stdcout << "fail5 " << result.first << " " << result.second << "\n";
913         return false;
914       }
915       edges.push_back(std::make_pair(half_edge(Point(0, 5), Point(0, 11)), (segment_id)edges.size()));
916       if(verify_scan(result, edges.begin(), edges.end())) {
917         stdcout << "fail6 " << result.first << " " << result.second << "\n";
918         return false;
919       }
920       edges.pop_back();
921       for(std::size_t i = 0; i < edges.size(); ++i) {
922         std::swap(edges[i].first.first, edges[i].first.second);
923       }
924       if(!verify_scan(result, edges.begin(), edges.end())) {
925         stdcout << "fail5 2 " << result.first << " " << result.second << "\n";
926         return false;
927       }
928       for(std::size_t i = 0; i < edges.size(); ++i) {
929         edges[i].first.first = Point(edges[i].first.first.get(HORIZONTAL) * -1,
930                                      edges[i].first.first.get(VERTICAL) * -1);
931         edges[i].first.second = Point(edges[i].first.second.get(HORIZONTAL) * -1,
932                                      edges[i].first.second.get(VERTICAL) * -1);
933       }
934       if(!verify_scan(result, edges.begin(), edges.end())) {
935         stdcout << "fail5 3 " << result.first << " " << result.second << "\n";
936         return false;
937       }
938       return true;
939     }
940 
941   };
942 
943   //scanline consumes the "flattened" fully intersected line segments produced by
944   //a pass of line_intersection along with property and count information and performs a
945   //useful operation like booleans or property merge or connectivity extraction
946   template <typename Unit, typename property_type, typename keytype = std::set<property_type> >
947   class scanline : public scanline_base<Unit> {
948   public:
949     //definitions
950     typedef typename scanline_base<Unit>::Point Point;
951 
952     //the first point is the vertex and and second point establishes the slope of an edge eminating from the vertex
953     //typedef std::pair<Point, Point> half_edge;
954     typedef typename scanline_base<Unit>::half_edge half_edge;
955 
956     //scanline comparator functor
957     typedef typename scanline_base<Unit>::less_half_edge less_half_edge;
958     typedef typename scanline_base<Unit>::less_point less_point;
959 
960     typedef keytype property_set;
961     //this is the data type used internally to store the combination of property counts at a given location
962     typedef std::vector<std::pair<property_type, int> > property_map;
963     //this data structure assocates a property and count to a half edge
964     typedef std::pair<half_edge, std::pair<property_type, int> > vertex_property;
965     //this data type is used internally to store the combined property data for a given half edge
966     typedef std::pair<half_edge, property_map> vertex_data;
967     //this data type stores the combination of many half edges
968     typedef std::vector<vertex_property> property_merge_data;
969     //this data structure stores end points of edges in the scanline
970     typedef std::set<Point, less_point> end_point_queue;
971 
972     //this is the output data type that is created by the scanline before it is post processed based on content of property sets
973     typedef std::pair<half_edge, std::pair<property_set, property_set> > half_edge_property;
974 
975     //this is the scanline data structure
976     typedef std::map<half_edge, property_map, less_half_edge> scanline_type;
977     typedef std::pair<half_edge, property_map> scanline_element;
978     typedef typename scanline_type::iterator iterator;
979     typedef typename scanline_type::const_iterator const_iterator;
980 
981     //data
982     scanline_type scan_data_;
983     std::vector<iterator> removal_set_; //edges to be removed at the current scanline stop
984     std::vector<scanline_element> insertion_set_; //edge to be inserted after current scanline stop
985     end_point_queue end_point_queue_;
986     Unit x_;
987     Unit y_;
988     int just_before_;
989     typename scanline_base<Unit>::evalAtXforYPack evalAtXforYPack_;
990   public:
scanline()991     inline scanline() : scan_data_(), removal_set_(), insertion_set_(), end_point_queue_(),
992                         x_((std::numeric_limits<Unit>::max)()), y_((std::numeric_limits<Unit>::max)()), just_before_(false), evalAtXforYPack_() {
993       less_half_edge lessElm(&x_, &just_before_, &evalAtXforYPack_);
994       scan_data_ = scanline_type(lessElm);
995     }
scanline(const scanline & that)996     inline scanline(const scanline& that) : scan_data_(), removal_set_(), insertion_set_(), end_point_queue_(),
997                                             x_((std::numeric_limits<Unit>::max)()), y_((std::numeric_limits<Unit>::max)()), just_before_(false), evalAtXforYPack_() {
998       (*this) = that; }
operator =(const scanline & that)999     inline scanline& operator=(const scanline& that) {
1000       x_ = that.x_;
1001       y_ = that.y_;
1002       just_before_ = that.just_before_;
1003       end_point_queue_ = that.end_point_queue_;
1004       //I cannot simply copy that.scanline_type to this scanline_type becuase the functor store pointers to other members!
1005       less_half_edge lessElm(&x_, &just_before_);
1006       scan_data_ = scanline_type(lessElm);
1007 
1008       scan_data_.insert(that.scan_data_.begin(), that.scan_data_.end());
1009       return *this;
1010     }
1011 
1012     template <typename result_type, typename result_functor>
write_out(result_type & result,result_functor rf,const half_edge & he,const property_map & pm_left,const property_map & pm_right)1013     void write_out(result_type& result, result_functor rf, const half_edge& he,
1014                    const property_map& pm_left, const property_map& pm_right) {
1015       //std::cout << "write out ";
1016       //std::cout << he.first << ", " << he.second << "\n";
1017       property_set ps_left, ps_right;
1018       set_unique_property(ps_left, pm_left);
1019       set_unique_property(ps_right, pm_right);
1020       if(ps_left != ps_right) {
1021         //std::cout << "!equivalent\n";
1022         rf(result, he, ps_left, ps_right);
1023       }
1024     }
1025 
1026     template <typename result_type, typename result_functor, typename iT>
handle_input_events(result_type & result,result_functor rf,iT begin,iT end)1027     iT handle_input_events(result_type& result, result_functor rf, iT begin, iT end) {
1028       //typedef typename high_precision_type<Unit>::type high_precision;
1029       //for each event
1030       property_map vertical_properties_above;
1031       property_map vertical_properties_below;
1032       half_edge vertical_edge_above;
1033       half_edge vertical_edge_below;
1034       std::vector<scanline_element> insertion_elements;
1035       //current_iter should increase monotonically toward end as we process scanline stop
1036       iterator current_iter = scan_data_.begin();
1037       just_before_ = true;
1038       Unit y = (std::numeric_limits<Unit>::min)();
1039       bool first_iteration = true;
1040       //we want to return from inside the loop when we hit end or new x
1041 #ifdef BOOST_POLYGON_MSVC
1042 #pragma warning (push)
1043 #pragma warning (disable: 4127)
1044 #endif
1045       while(true) {
1046         if(begin == end || (!first_iteration && ((*begin).first.first.get(VERTICAL) != y ||
1047                                                  (*begin).first.first.get(HORIZONTAL) != x_))) {
1048           //lookup iterator range in scanline for elements coming in from the left
1049           //that end at this y
1050           Point pt(x_, y);
1051           //grab the properties coming in from below
1052           property_map properties_below;
1053           if(current_iter != scan_data_.end()) {
1054             //make sure we are looking at element in scanline just below y
1055             //if(evalAtXforY(x_, (*current_iter).first.first, (*current_iter).first.second) != y) {
1056             if(scanline_base<Unit>::on_above_or_below(Point(x_, y), (*current_iter).first) != 0) {
1057               Point e2(pt);
1058               if(e2.get(VERTICAL) != (std::numeric_limits<Unit>::max)())
1059                 e2.set(VERTICAL, e2.get(VERTICAL) + 1);
1060               else
1061                 e2.set(VERTICAL, e2.get(VERTICAL) - 1);
1062               half_edge vhe(pt, e2);
1063               current_iter = scan_data_.lower_bound(vhe);
1064             }
1065             if(current_iter != scan_data_.end()) {
1066               //get the bottom iterator for elements at this point
1067               //while(evalAtXforY(x_, (*current_iter).first.first, (*current_iter).first.second) >= (high_precision)y &&
1068               while(scanline_base<Unit>::on_above_or_below(Point(x_, y), (*current_iter).first) != 1 &&
1069                     current_iter != scan_data_.begin()) {
1070                 --current_iter;
1071               }
1072               //if(evalAtXforY(x_, (*current_iter).first.first, (*current_iter).first.second) >= (high_precision)y) {
1073               if(scanline_base<Unit>::on_above_or_below(Point(x_, y), (*current_iter).first) != 1) {
1074                 properties_below.clear();
1075               } else {
1076                 properties_below = (*current_iter).second;
1077                 //move back up to y or one past y
1078                 ++current_iter;
1079               }
1080             }
1081           }
1082           std::vector<iterator> edges_from_left;
1083           while(current_iter != scan_data_.end() &&
1084                 //can only be true if y is integer
1085                 //evalAtXforY(x_, (*current_iter).first.first, (*current_iter).first.second) == y) {
1086                 scanline_base<Unit>::on_above_or_below(Point(x_, y), (*current_iter).first) == 0) {
1087             //removal_set_.push_back(current_iter);
1088             ++current_iter;
1089           }
1090           //merge vertical count with count from below
1091           if(!vertical_properties_below.empty()) {
1092             merge_property_maps(vertical_properties_below, properties_below);
1093             //write out vertical edge
1094             write_out(result, rf, vertical_edge_below, properties_below, vertical_properties_below);
1095           } else {
1096             merge_property_maps(vertical_properties_below, properties_below);
1097           }
1098           //iteratively add intertion element counts to count from below
1099           //and write them to insertion set
1100           for(std::size_t i = 0; i < insertion_elements.size(); ++i) {
1101             if(i == 0) {
1102               merge_property_maps(insertion_elements[i].second, vertical_properties_below);
1103               write_out(result, rf, insertion_elements[i].first, insertion_elements[i].second, vertical_properties_below);
1104             } else {
1105               merge_property_maps(insertion_elements[i].second, insertion_elements[i-1].second);
1106               write_out(result, rf, insertion_elements[i].first, insertion_elements[i].second, insertion_elements[i-1].second);
1107             }
1108             insertion_set_.push_back(insertion_elements[i]);
1109           }
1110           if((begin == end || (*begin).first.first.get(HORIZONTAL) != x_)) {
1111             if(vertical_properties_above.empty()) {
1112               return begin;
1113             } else {
1114               y = vertical_edge_above.second.get(VERTICAL);
1115               vertical_properties_below.clear();
1116               vertical_properties_above.swap(vertical_properties_below);
1117               vertical_edge_below = vertical_edge_above;
1118               insertion_elements.clear();
1119               continue;
1120             }
1121           }
1122           vertical_properties_below.clear();
1123           vertical_properties_above.swap(vertical_properties_below);
1124           vertical_edge_below = vertical_edge_above;
1125           insertion_elements.clear();
1126         }
1127         if(begin != end) {
1128           const vertex_property& vp = *begin;
1129           const half_edge& he = vp.first;
1130           y = he.first.get(VERTICAL);
1131           first_iteration = false;
1132           if(! vertical_properties_below.empty() &&
1133              vertical_edge_below.second.get(VERTICAL) < y) {
1134             y = vertical_edge_below.second.get(VERTICAL);
1135             continue;
1136           }
1137           if(scanline_base<Unit>::is_vertical(he)) {
1138             update_property_map(vertical_properties_above, vp.second);
1139             vertical_edge_above = he;
1140           } else {
1141             if(insertion_elements.empty() ||
1142                insertion_elements.back().first != he) {
1143               insertion_elements.push_back(scanline_element(he, property_map()));
1144             }
1145             update_property_map(insertion_elements.back().second, vp.second);
1146           }
1147           ++begin;
1148         }
1149       }
1150 #ifdef BOOST_POLYGON_MSVC
1151 #pragma warning (pop)
1152 #endif
1153 
1154     }
1155 
erase_end_events(typename end_point_queue::iterator epqi)1156     inline void erase_end_events(typename end_point_queue::iterator epqi) {
1157       end_point_queue_.erase(end_point_queue_.begin(), epqi);
1158       for(typename std::vector<iterator>::iterator retire_itr = removal_set_.begin();
1159           retire_itr != removal_set_.end(); ++retire_itr) {
1160         scan_data_.erase(*retire_itr);
1161       }
1162       removal_set_.clear();
1163     }
1164 
1165 
remove_retired_edges_from_scanline()1166     inline void remove_retired_edges_from_scanline() {
1167       just_before_ = true;
1168       typename end_point_queue::iterator epqi = end_point_queue_.begin();
1169       Unit current_x = x_;
1170       Unit previous_x = x_;
1171       while(epqi != end_point_queue_.end() &&
1172             (*epqi).get(HORIZONTAL) <= current_x) {
1173         x_ = (*epqi).get(HORIZONTAL);
1174         if(x_ != previous_x) erase_end_events(epqi);
1175         previous_x = x_;
1176         //lookup elements
1177         Point e2(*epqi);
1178         if(e2.get(VERTICAL) != (std::numeric_limits<Unit>::max)())
1179           e2.set(VERTICAL, e2.get(VERTICAL) + 1);
1180         else
1181           e2.set(VERTICAL, e2.get(VERTICAL) - 1);
1182         half_edge vhe_e(*epqi, e2);
1183         iterator current_iter = scan_data_.lower_bound(vhe_e);
1184         while(current_iter != scan_data_.end() && (*current_iter).first.second == (*epqi)) {
1185           //evalAtXforY(x_, (*current_iter).first.first, (*current_iter).first.second) == (*epqi).get(VERTICAL)) {
1186           removal_set_.push_back(current_iter);
1187           ++current_iter;
1188         }
1189         ++epqi;
1190       }
1191       x_ = current_x;
1192       erase_end_events(epqi);
1193     }
1194 
insert_new_edges_into_scanline()1195     inline void insert_new_edges_into_scanline() {
1196       just_before_ = false;
1197       for(typename std::vector<scanline_element>::iterator insert_itr = insertion_set_.begin();
1198           insert_itr != insertion_set_.end(); ++insert_itr) {
1199         scan_data_.insert(*insert_itr);
1200         end_point_queue_.insert((*insert_itr).first.second);
1201       }
1202       insertion_set_.clear();
1203     }
1204 
1205     //iterator over range of vertex property elements and call result functor
1206     //passing edge to be output, the merged data on both sides and the result
1207     template <typename result_type, typename result_functor, typename iT>
scan(result_type & result,result_functor rf,iT begin,iT end)1208     void scan(result_type& result, result_functor rf, iT begin, iT end) {
1209       while(begin != end) {
1210         x_ = (*begin).first.first.get(HORIZONTAL); //update scanline stop location
1211         //print_scanline();
1212         --x_;
1213         remove_retired_edges_from_scanline();
1214         ++x_;
1215         begin = handle_input_events(result, rf, begin, end);
1216         remove_retired_edges_from_scanline();
1217         //print_scanline();
1218         insert_new_edges_into_scanline();
1219       }
1220       //print_scanline();
1221       x_ = (std::numeric_limits<Unit>::max)();
1222       remove_retired_edges_from_scanline();
1223     }
1224 
1225     //inline void print_scanline() {
1226     //  std::cout << "scanline at " << x_ << ": ";
1227     //  for(iterator itr = scan_data_.begin(); itr != scan_data_.end(); ++itr) {
1228     //    const scanline_element& se = *itr;
1229     //    const half_edge& he = se.first;
1230     //    const property_map& mp = se.second;
1231     //    std::cout << he.first << ", " << he.second << " ( ";
1232     //    for(std::size_t i = 0; i < mp.size(); ++i) {
1233     //      std::cout << mp[i].first << ":" << mp[i].second << " ";
1234     //    } std::cout << ") ";
1235     //  } std::cout << "\n";
1236     //}
1237 
merge_property_maps(property_map & mp,const property_map & mp2)1238     static inline void merge_property_maps(property_map& mp, const property_map& mp2) {
1239       property_map newmp;
1240       newmp.reserve(mp.size() + mp2.size());
1241       unsigned int i = 0;
1242       unsigned int j = 0;
1243       while(i != mp.size() && j != mp2.size()) {
1244         if(mp[i].first < mp2[j].first) {
1245           newmp.push_back(mp[i]);
1246           ++i;
1247         } else if(mp[i].first > mp2[j].first) {
1248           newmp.push_back(mp2[j]);
1249           ++j;
1250         } else {
1251           int count = mp[i].second;
1252           count += mp2[j].second;
1253           if(count) {
1254             newmp.push_back(mp[i]);
1255             newmp.back().second = count;
1256           }
1257           ++i;
1258           ++j;
1259         }
1260       }
1261       while(i != mp.size()) {
1262         newmp.push_back(mp[i]);
1263         ++i;
1264       }
1265       while(j != mp2.size()) {
1266         newmp.push_back(mp2[j]);
1267         ++j;
1268       }
1269       mp.swap(newmp);
1270     }
1271 
update_property_map(property_map & mp,const std::pair<property_type,int> & prop_data)1272     static inline void update_property_map(property_map& mp, const std::pair<property_type, int>& prop_data) {
1273       property_map newmp;
1274       newmp.reserve(mp.size() +1);
1275       bool consumed = false;
1276       for(std::size_t i = 0; i < mp.size(); ++i) {
1277         if(!consumed && prop_data.first == mp[i].first) {
1278           consumed = true;
1279           int count = prop_data.second + mp[i].second;
1280           if(count)
1281             newmp.push_back(std::make_pair(prop_data.first, count));
1282         } else if(!consumed && prop_data.first < mp[i].first) {
1283           consumed = true;
1284           newmp.push_back(prop_data);
1285           newmp.push_back(mp[i]);
1286         } else {
1287           newmp.push_back(mp[i]);
1288         }
1289       }
1290       if(!consumed) newmp.push_back(prop_data);
1291       mp.swap(newmp);
1292     }
1293 
set_unique_property(property_set & unqiue_property,const property_map & property)1294     static inline void set_unique_property(property_set& unqiue_property, const property_map& property) {
1295       unqiue_property.clear();
1296       for(typename property_map::const_iterator itr = property.begin(); itr != property.end(); ++itr) {
1297         if((*itr).second > 0)
1298           unqiue_property.insert(unqiue_property.end(), (*itr).first);
1299       }
1300     }
1301 
common_vertex(const half_edge & he1,const half_edge & he2)1302     static inline bool common_vertex(const half_edge& he1, const half_edge& he2) {
1303       return he1.first == he2.first ||
1304         he1.first == he2.second ||
1305         he1.second == he2.first ||
1306         he1.second == he2.second;
1307     }
1308 
1309     typedef typename scanline_base<Unit>::vertex_half_edge vertex_half_edge;
1310     template <typename iT>
convert_segments_to_vertex_half_edges(std::vector<vertex_half_edge> & output,iT begin,iT end)1311     static inline void convert_segments_to_vertex_half_edges(std::vector<vertex_half_edge>& output, iT begin, iT end) {
1312       for( ; begin != end; ++begin) {
1313         const half_edge& he = (*begin).first;
1314         int count = (*begin).second;
1315         output.push_back(vertex_half_edge(he.first, he.second, count));
1316         output.push_back(vertex_half_edge(he.second, he.first, -count));
1317       }
1318       polygon_sort(output.begin(), output.end());
1319     }
1320 
1321     class test_functor {
1322     public:
test_functor()1323       inline test_functor() {}
operator ()(std::vector<std::pair<half_edge,std::pair<property_set,property_set>>> & result,const half_edge & he,const property_set & ps_left,const property_set & ps_right)1324       inline void operator()(std::vector<std::pair<half_edge, std::pair<property_set, property_set> > >& result,
1325                              const half_edge& he, const property_set& ps_left, const property_set& ps_right) {
1326         result.push_back(std::make_pair(he, std::make_pair(ps_left, ps_right)));
1327       }
1328     };
1329     template <typename stream_type>
test_scanline(stream_type & stdcout)1330     static inline bool test_scanline(stream_type& stdcout) {
1331       std::vector<std::pair<half_edge, std::pair<property_set, property_set> > > result;
1332       std::vector<std::pair<half_edge, std::pair<property_type, int> > > input;
1333       input.push_back(std::make_pair(half_edge(Point(0, 0), Point(0, 10)), std::make_pair(0, 1)));
1334       input.push_back(std::make_pair(half_edge(Point(0, 0), Point(10, 0)), std::make_pair(0, 1)));
1335       input.push_back(std::make_pair(half_edge(Point(0, 10), Point(10, 10)), std::make_pair(0, -1)));
1336       input.push_back(std::make_pair(half_edge(Point(10, 0), Point(10, 10)), std::make_pair(0, -1)));
1337       scanline sl;
1338       test_functor tf;
1339       sl.scan(result, tf, input.begin(), input.end());
1340       stdcout << "scanned\n";
1341       for(std::size_t i = 0; i < result.size(); ++i) {
1342         stdcout << result[i].first.first << ", " << result[i].first.second << "; ";
1343       } stdcout << "\n";
1344       input.clear();
1345       result.clear();
1346       input.push_back(std::make_pair(half_edge(Point(-1, -1), Point(10, 0)), std::make_pair(0, 1)));
1347       input.push_back(std::make_pair(half_edge(Point(-1, -1), Point(0, 10)), std::make_pair(0, -1)));
1348       input.push_back(std::make_pair(half_edge(Point(0, 10), Point(11, 11)), std::make_pair(0, -1)));
1349       input.push_back(std::make_pair(half_edge(Point(10, 0), Point(11, 11)), std::make_pair(0, 1)));
1350       scanline sl2;
1351       sl2.scan(result, tf, input.begin(), input.end());
1352       stdcout << "scanned\n";
1353       for(std::size_t i = 0; i < result.size(); ++i) {
1354         stdcout << result[i].first.first << ", " << result[i].first.second << "; ";
1355       } stdcout << "\n";
1356       input.clear();
1357       result.clear();
1358       input.push_back(std::make_pair(half_edge(Point(0, 0), Point(0, 10)), std::make_pair(0, 1)));
1359       input.push_back(std::make_pair(half_edge(Point(0, 0), Point(10, 0)), std::make_pair(0, 1)));
1360       input.push_back(std::make_pair(half_edge(Point(0, 10), Point(10, 10)), std::make_pair(0, -1)));
1361       input.push_back(std::make_pair(half_edge(Point(1, 1), Point(8, 2)), std::make_pair(1, 1)));
1362       input.push_back(std::make_pair(half_edge(Point(1, 1), Point(2, 8)), std::make_pair(1, -1)));
1363       input.push_back(std::make_pair(half_edge(Point(2, 8), Point(9, 9)), std::make_pair(1, -1)));
1364       input.push_back(std::make_pair(half_edge(Point(8, 2), Point(9, 9)), std::make_pair(1, 1)));
1365       input.push_back(std::make_pair(half_edge(Point(10, 0), Point(10, 10)), std::make_pair(0, -1)));
1366       scanline sl3;
1367       sl3.scan(result, tf, input.begin(), input.end());
1368       stdcout << "scanned\n";
1369       for(std::size_t i = 0; i < result.size(); ++i) {
1370         stdcout << result[i].first.first << ", " << result[i].first.second << "; ";
1371       } stdcout << "\n";
1372       input.clear();
1373       result.clear();
1374       input.push_back(std::make_pair(half_edge(Point(0, 0), Point(0, 10)), std::make_pair(0, 1)));
1375       input.push_back(std::make_pair(half_edge(Point(0, 0), Point(10, 0)), std::make_pair(0, 1)));
1376       input.push_back(std::make_pair(half_edge(Point(0, 10), Point(10, 10)), std::make_pair(0, -1)));
1377       input.push_back(std::make_pair(half_edge(Point(1, 1), Point(8, 2)), std::make_pair(0, 1)));
1378       input.push_back(std::make_pair(half_edge(Point(1, 1), Point(2, 8)), std::make_pair(0, -1)));
1379       input.push_back(std::make_pair(half_edge(Point(2, 8), Point(9, 9)), std::make_pair(0, -1)));
1380       input.push_back(std::make_pair(half_edge(Point(8, 2), Point(9, 9)), std::make_pair(0, 1)));
1381       input.push_back(std::make_pair(half_edge(Point(10, 0), Point(10, 10)), std::make_pair(0, -1)));
1382       scanline sl4;
1383       sl4.scan(result, tf, input.begin(), input.end());
1384       stdcout << "scanned\n";
1385       for(std::size_t i = 0; i < result.size(); ++i) {
1386         stdcout << result[i].first.first << ", " << result[i].first.second << "; ";
1387       } stdcout << "\n";
1388       input.clear();
1389       result.clear();
1390       input.push_back(std::make_pair(half_edge(Point(0, 0), Point(10, 0)), std::make_pair(0, 1)));
1391       input.push_back(std::make_pair(half_edge(Point(0, 0), Point(9, 1)), std::make_pair(0, 1)));
1392       input.push_back(std::make_pair(half_edge(Point(0, 0), Point(1, 9)), std::make_pair(0, -1)));
1393       input.push_back(std::make_pair(half_edge(Point(0, 0), Point(0, 10)), std::make_pair(0, 1)));
1394       input.push_back(std::make_pair(half_edge(Point(0, 10), Point(10, 10)), std::make_pair(0, -1)));
1395       input.push_back(std::make_pair(half_edge(Point(1, 9), Point(10, 10)), std::make_pair(0, -1)));
1396       input.push_back(std::make_pair(half_edge(Point(9, 1), Point(10, 10)), std::make_pair(0, 1)));
1397       input.push_back(std::make_pair(half_edge(Point(10, 0), Point(10, 10)), std::make_pair(0, -1)));
1398       scanline sl5;
1399       sl5.scan(result, tf, input.begin(), input.end());
1400       stdcout << "scanned\n";
1401       for(std::size_t i = 0; i < result.size(); ++i) {
1402         stdcout << result[i].first.first << ", " << result[i].first.second << "; ";
1403       } stdcout << "\n";
1404       input.clear();
1405       result.clear();
1406       input.push_back(std::make_pair(half_edge(Point(0, 0), Point(10, 0)), std::make_pair(0, 1)));
1407       input.push_back(std::make_pair(half_edge(Point(0, 0), Point(9, 1)), std::make_pair(1, 1)));
1408       input.push_back(std::make_pair(half_edge(Point(0, 0), Point(1, 9)), std::make_pair(1, -1)));
1409       input.push_back(std::make_pair(half_edge(Point(0, 0), Point(0, 10)), std::make_pair(0, 1)));
1410       input.push_back(std::make_pair(half_edge(Point(0, 10), Point(10, 10)), std::make_pair(0, -1)));
1411       input.push_back(std::make_pair(half_edge(Point(1, 9), Point(10, 10)), std::make_pair(1, -1)));
1412       input.push_back(std::make_pair(half_edge(Point(9, 1), Point(10, 10)), std::make_pair(1, 1)));
1413       input.push_back(std::make_pair(half_edge(Point(10, 0), Point(10, 10)), std::make_pair(0, -1)));
1414       scanline sl6;
1415       sl6.scan(result, tf, input.begin(), input.end());
1416       stdcout << "scanned\n";
1417       for(std::size_t i = 0; i < result.size(); ++i) {
1418         stdcout << result[i].first.first << ", " << result[i].first.second << "; ";
1419       } stdcout << "\n";
1420       input.clear();
1421       result.clear();
1422       input.push_back(std::make_pair(half_edge(Point(0, 0), Point(10, 0)), std::make_pair(0, 1)));
1423       input.push_back(std::make_pair(half_edge(Point(0, 0), Point(9, 1)), std::make_pair(1, 1)));
1424       input.push_back(std::make_pair(half_edge(Point(0, 0), Point(1, 9)), std::make_pair(1, -1)));
1425       input.push_back(std::make_pair(half_edge(Point(0, 0), Point(0, 10)), std::make_pair(0, 1)));
1426       input.push_back(std::make_pair(half_edge(Point(0, 10), Point(10, 10)), std::make_pair(0, -1)));
1427       input.push_back(std::make_pair(half_edge(Point(0, 20), Point(10, 20)), std::make_pair(0, 1)));
1428       input.push_back(std::make_pair(half_edge(Point(0, 20), Point(9, 21)), std::make_pair(1, 1)));
1429       input.push_back(std::make_pair(half_edge(Point(0, 20), Point(1, 29)), std::make_pair(1, -1)));
1430       input.push_back(std::make_pair(half_edge(Point(0, 20), Point(0, 30)), std::make_pair(0, 1)));
1431       input.push_back(std::make_pair(half_edge(Point(0, 30), Point(10, 30)), std::make_pair(0, -1)));
1432       input.push_back(std::make_pair(half_edge(Point(1, 9), Point(10, 10)), std::make_pair(1, -1)));
1433       input.push_back(std::make_pair(half_edge(Point(1, 29), Point(10, 30)), std::make_pair(1, -1)));
1434       input.push_back(std::make_pair(half_edge(Point(9, 1), Point(10, 10)), std::make_pair(1, 1)));
1435       input.push_back(std::make_pair(half_edge(Point(9, 21), Point(10, 30)), std::make_pair(1, 1)));
1436       input.push_back(std::make_pair(half_edge(Point(10, 20), Point(10, 30)), std::make_pair(0, -1)));
1437       input.push_back(std::make_pair(half_edge(Point(10, 20), Point(10, 30)), std::make_pair(0, -1)));
1438       scanline sl7;
1439       sl7.scan(result, tf, input.begin(), input.end());
1440       stdcout << "scanned\n";
1441       for(std::size_t i = 0; i < result.size(); ++i) {
1442         stdcout << result[i].first.first << ", " << result[i].first.second << "; ";
1443       } stdcout << "\n";
1444       input.clear();
1445       result.clear();
1446       input.push_back(std::make_pair(half_edge(Point(-1, -1), Point(10, 0)), std::make_pair(0, 1))); //a
1447       input.push_back(std::make_pair(half_edge(Point(-1, -1), Point(0, 10)), std::make_pair(0, -1))); //a
1448       input.push_back(std::make_pair(half_edge(Point(0, 10), Point(11, 11)), std::make_pair(0, -1))); //a
1449       input.push_back(std::make_pair(half_edge(Point(10, 0), Point(20, 0)), std::make_pair(0, 1))); //b
1450       input.push_back(std::make_pair(half_edge(Point(10, 0), Point(11, 11)), std::make_pair(0, -1))); //b
1451       input.push_back(std::make_pair(half_edge(Point(10, 0), Point(11, 11)), std::make_pair(0, 1))); //a
1452       input.push_back(std::make_pair(half_edge(Point(11, 11), Point(20, 10)), std::make_pair(0, -1))); //b
1453       input.push_back(std::make_pair(half_edge(Point(20, 0), Point(30, 0)), std::make_pair(0, 1))); //c
1454       input.push_back(std::make_pair(half_edge(Point(20, 0), Point(20, 10)), std::make_pair(0, -1))); //b
1455       input.push_back(std::make_pair(half_edge(Point(20, 0), Point(20, 10)), std::make_pair(0, 1))); //c
1456       input.push_back(std::make_pair(half_edge(Point(20, 10), Point(30, 10)), std::make_pair(0, -1))); //c
1457       input.push_back(std::make_pair(half_edge(Point(30, 0), Point(30, 10)), std::make_pair(0, -1))); //c
1458       scanline sl8;
1459       sl8.scan(result, tf, input.begin(), input.end());
1460       stdcout << "scanned\n";
1461       for(std::size_t i = 0; i < result.size(); ++i) {
1462         stdcout << result[i].first.first << ", " << result[i].first.second << "; ";
1463       } stdcout << "\n";
1464       return true;
1465     }
1466 
1467   };
1468 
1469   template <typename Unit>
1470   class merge_output_functor {
1471   public:
1472     typedef typename scanline_base<Unit>::half_edge half_edge;
merge_output_functor()1473     merge_output_functor() {}
1474     template <typename result_type, typename key_type>
operator ()(result_type & result,const half_edge & edge,const key_type & left,const key_type & right)1475     void operator()(result_type& result, const half_edge& edge, const key_type& left, const key_type& right) {
1476       typename std::pair<half_edge, int> elem;
1477       elem.first = edge;
1478       elem.second = 1;
1479       if(edge.second < edge.first) elem.second *= -1;
1480       if(scanline_base<Unit>::is_vertical(edge)) elem.second *= -1;
1481       if(!left.empty())
1482         result[left].insert_clean(elem);
1483       elem.second *= -1;
1484       if(!right.empty())
1485         result[right].insert_clean(elem);
1486     }
1487   };
1488 
1489   template <typename Unit, typename property_type, typename key_type = std::set<property_type>,
1490             typename output_functor_type = merge_output_functor<Unit> >
1491   class property_merge : public scanline_base<Unit> {
1492   protected:
1493     typedef typename scanline_base<Unit>::Point Point;
1494 
1495     //the first point is the vertex and and second point establishes the slope of an edge eminating from the vertex
1496     //typedef std::pair<Point, Point> half_edge;
1497     typedef typename scanline_base<Unit>::half_edge half_edge;
1498 
1499     //scanline comparator functor
1500     typedef typename scanline_base<Unit>::less_half_edge less_half_edge;
1501     typedef typename scanline_base<Unit>::less_point less_point;
1502 
1503     //this data structure assocates a property and count to a half edge
1504     typedef std::pair<half_edge, std::pair<property_type, int> > vertex_property;
1505     //this data type stores the combination of many half edges
1506     typedef std::vector<vertex_property> property_merge_data;
1507 
1508     //this is the data type used internally to store the combination of property counts at a given location
1509     typedef std::vector<std::pair<property_type, int> > property_map;
1510     //this data type is used internally to store the combined property data for a given half edge
1511     typedef std::pair<half_edge, property_map> vertex_data;
1512 
1513     property_merge_data pmd;
1514     typename scanline_base<Unit>::evalAtXforYPack evalAtXforYPack_;
1515 
1516     template<typename vertex_data_type>
1517     class less_vertex_data {
1518       typename scanline_base<Unit>::evalAtXforYPack* pack_;
1519     public:
less_vertex_data()1520       less_vertex_data() : pack_() {}
less_vertex_data(typename scanline_base<Unit>::evalAtXforYPack * pack)1521       less_vertex_data(typename scanline_base<Unit>::evalAtXforYPack* pack) : pack_(pack) {}
operator ()(const vertex_data_type & lvalue,const vertex_data_type & rvalue) const1522       bool operator() (const vertex_data_type& lvalue, const vertex_data_type& rvalue) const {
1523         less_point lp;
1524         if(lp(lvalue.first.first, rvalue.first.first)) return true;
1525         if(lp(rvalue.first.first, lvalue.first.first)) return false;
1526         Unit x = lvalue.first.first.get(HORIZONTAL);
1527         int just_before_ = 0;
1528         less_half_edge lhe(&x, &just_before_, pack_);
1529         return lhe(lvalue.first, rvalue.first);
1530       }
1531     };
1532 
1533 
sort_property_merge_data()1534     inline void sort_property_merge_data() {
1535       less_vertex_data<vertex_property> lvd(&evalAtXforYPack_);
1536       polygon_sort(pmd.begin(), pmd.end(), lvd);
1537     }
1538   public:
get_property_merge_data()1539     inline property_merge_data& get_property_merge_data() { return pmd; }
property_merge()1540     inline property_merge() : pmd(), evalAtXforYPack_() {}
property_merge(const property_merge & pm)1541     inline property_merge(const property_merge& pm) : pmd(pm.pmd), evalAtXforYPack_(pm.evalAtXforYPack_) {}
operator =(const property_merge & pm)1542     inline property_merge& operator=(const property_merge& pm) { pmd = pm.pmd; return *this; }
1543 
1544     template <typename polygon_type>
insert(const polygon_type & polygon_object,const property_type & property_value,bool is_hole=false)1545     void insert(const polygon_type& polygon_object, const property_type& property_value, bool is_hole = false) {
1546       insert(polygon_object, property_value, is_hole, typename geometry_concept<polygon_type>::type());
1547     }
1548 
1549     //result type should be std::map<std::set<property_type>, polygon_set_type>
1550     //or std::map<std::vector<property_type>, polygon_set_type>
1551     template <typename result_type>
merge(result_type & result)1552     void merge(result_type& result) {
1553       if(pmd.empty()) return;
1554       //intersect data
1555       property_merge_data tmp_pmd;
1556       line_intersection<Unit>::validate_scan(tmp_pmd, pmd.begin(), pmd.end());
1557       pmd.swap(tmp_pmd);
1558       sort_property_merge_data();
1559       scanline<Unit, property_type, key_type> sl;
1560       output_functor_type mof;
1561       sl.scan(result, mof, pmd.begin(), pmd.end());
1562     }
1563 
verify1()1564     inline bool verify1() {
1565       std::pair<int, int> offenders;
1566       std::vector<std::pair<half_edge, int> > lines;
1567       int count = 0;
1568       for(std::size_t i = 0; i < pmd.size(); ++i) {
1569         lines.push_back(std::make_pair(pmd[i].first, count++));
1570       }
1571       if(!line_intersection<Unit>::verify_scan(offenders, lines.begin(), lines.end())) {
1572         //stdcout << "Intersection failed!\n";
1573         //stdcout << offenders.first << " " << offenders.second << "\n";
1574         return false;
1575       }
1576       std::vector<Point> pts;
1577       for(std::size_t i = 0; i < lines.size(); ++i) {
1578         pts.push_back(lines[i].first.first);
1579         pts.push_back(lines[i].first.second);
1580       }
1581       polygon_sort(pts.begin(), pts.end());
1582       for(std::size_t i = 0; i < pts.size(); i+=2) {
1583         if(pts[i] != pts[i+1]) {
1584           //stdcout << "Non-closed figures after line intersection!\n";
1585           return false;
1586         }
1587       }
1588       return true;
1589     }
1590 
clear()1591     void clear() {*this = property_merge();}
1592 
1593   protected:
1594     template <typename polygon_type>
insert(const polygon_type & polygon_object,const property_type & property_value,bool is_hole,polygon_concept)1595     void insert(const polygon_type& polygon_object, const property_type& property_value, bool is_hole,
1596                 polygon_concept ) {
1597       bool first_iteration = true;
1598       bool second_iteration = true;
1599       Point first_point;
1600       Point second_point;
1601       Point previous_previous_point;
1602       Point previous_point;
1603       Point current_point;
1604       direction_1d winding_dir = winding(polygon_object);
1605       for(typename polygon_traits<polygon_type>::iterator_type itr = begin_points(polygon_object);
1606           itr != end_points(polygon_object); ++itr) {
1607         assign(current_point, *itr);
1608         if(first_iteration) {
1609           first_iteration = false;
1610           first_point = previous_point = current_point;
1611         } else if(second_iteration) {
1612           if(previous_point != current_point) {
1613             second_iteration = false;
1614             previous_previous_point = previous_point;
1615             second_point = previous_point = current_point;
1616           }
1617         } else {
1618           if(previous_point != current_point) {
1619             create_vertex(pmd, previous_point, current_point, winding_dir,
1620                           is_hole, property_value);
1621             previous_previous_point = previous_point;
1622             previous_point = current_point;
1623           }
1624         }
1625       }
1626       current_point = first_point;
1627       if(!first_iteration && !second_iteration) {
1628         if(previous_point != current_point) {
1629           create_vertex(pmd, previous_point, current_point, winding_dir,
1630                         is_hole, property_value);
1631           previous_previous_point = previous_point;
1632           previous_point = current_point;
1633         }
1634         current_point = second_point;
1635         create_vertex(pmd, previous_point, current_point, winding_dir,
1636                       is_hole, property_value);
1637         previous_previous_point = previous_point;
1638         previous_point = current_point;
1639       }
1640     }
1641 
1642     template <typename polygon_with_holes_type>
insert(const polygon_with_holes_type & polygon_with_holes_object,const property_type & property_value,bool is_hole,polygon_with_holes_concept)1643     void insert(const polygon_with_holes_type& polygon_with_holes_object, const property_type& property_value, bool is_hole,
1644                 polygon_with_holes_concept) {
1645       insert(polygon_with_holes_object, property_value, is_hole, polygon_concept());
1646       for(typename polygon_with_holes_traits<polygon_with_holes_type>::iterator_holes_type itr =
1647             begin_holes(polygon_with_holes_object);
1648           itr != end_holes(polygon_with_holes_object); ++itr) {
1649         insert(*itr, property_value, !is_hole, polygon_concept());
1650       }
1651     }
1652 
1653     template <typename rectangle_type>
insert(const rectangle_type & rectangle_object,const property_type & property_value,bool is_hole,rectangle_concept)1654     void insert(const rectangle_type& rectangle_object, const property_type& property_value, bool is_hole,
1655                 rectangle_concept ) {
1656       polygon_90_data<Unit> poly;
1657       assign(poly, rectangle_object);
1658       insert(poly, property_value, is_hole, polygon_concept());
1659     }
1660 
1661   public: //change to private when done testing
1662 
create_vertex(property_merge_data & pmd,const Point & current_point,const Point & next_point,direction_1d winding,bool is_hole,const property_type & property)1663     static inline void create_vertex(property_merge_data& pmd,
1664                                      const Point& current_point,
1665                                      const Point& next_point,
1666                                      direction_1d winding,
1667                                      bool is_hole, const property_type& property) {
1668       if(current_point == next_point) return;
1669       vertex_property current_vertex;
1670       current_vertex.first.first = current_point;
1671       current_vertex.first.second = next_point;
1672       current_vertex.second.first = property;
1673       int multiplier = 1;
1674       if(winding == CLOCKWISE)
1675         multiplier = -1;
1676       if(is_hole)
1677         multiplier *= -1;
1678       if(current_point < next_point) {
1679         multiplier *= -1;
1680         std::swap(current_vertex.first.first, current_vertex.first.second);
1681       }
1682       current_vertex.second.second = multiplier * (euclidean_distance(next_point, current_point, HORIZONTAL) == 0 ? -1: 1);
1683       pmd.push_back(current_vertex);
1684       //current_vertex.first.second = previous_point;
1685       //current_vertex.second.second *= -1;
1686       //pmd.push_back(current_vertex);
1687     }
1688 
sort_vertex_half_edges(vertex_data & vertex)1689     static inline void sort_vertex_half_edges(vertex_data& vertex) {
1690       less_half_edge_pair lessF(vertex.first);
1691       polygon_sort(vertex.second.begin(), vertex.second.end(), lessF);
1692     }
1693 
1694     class less_half_edge_pair {
1695     private:
1696       Point pt_;
1697     public:
less_half_edge_pair(const Point & pt)1698       less_half_edge_pair(const Point& pt) : pt_(pt) {}
operator ()(const half_edge & e1,const half_edge & e2)1699       bool operator()(const half_edge& e1, const half_edge& e2) {
1700         const Point& pt1 = e1.first;
1701         const Point& pt2 = e2.first;
1702         if(get(pt1, HORIZONTAL) ==
1703            get(pt_, HORIZONTAL)) {
1704           //vertical edge is always largest
1705           return false;
1706         }
1707         if(get(pt2, HORIZONTAL) ==
1708            get(pt_, HORIZONTAL)) {
1709           //if half edge 1 is not vertical its slope is less than that of half edge 2
1710           return get(pt1, HORIZONTAL) != get(pt2, HORIZONTAL);
1711         }
1712         return scanline_base<Unit>::less_slope(get(pt_, HORIZONTAL),
1713                                                get(pt_, VERTICAL), pt1, pt2);
1714       }
1715     };
1716 
1717   public:
1718     //test functions
1719     template <typename stream_type>
print(stream_type & o,const property_map & c)1720     static stream_type& print (stream_type& o, const property_map& c)
1721     {
1722       o << "count: {";
1723       for(typename property_map::const_iterator itr = c.begin(); itr != c.end(); ++itr) {
1724         o << ((*itr).first) << ":" << ((*itr).second) << " ";
1725       }
1726       return o << "} ";
1727     }
1728 
1729 
1730     template <typename stream_type>
print(stream_type & o,const half_edge & he)1731     static stream_type& print (stream_type& o, const half_edge& he)
1732     {
1733       o << "half edge: (";
1734       o << (he.first);
1735       return o << ", " << (he.second) << ") ";
1736     }
1737 
1738     template <typename stream_type>
print(stream_type & o,const vertex_property & c)1739     static stream_type& print (stream_type& o, const vertex_property& c)
1740     {
1741       o << "vertex property: {";
1742       print(o, c.first);
1743       o << ", " << c.second.first << ":" << c.second.second << " ";
1744       return o;
1745     }
1746 
1747     template <typename stream_type>
print(stream_type & o,const std::vector<vertex_property> & hev)1748     static stream_type& print (stream_type& o, const std::vector<vertex_property>& hev)
1749     {
1750       o << "vertex properties: {";
1751       for(std::size_t i = 0; i < hev.size(); ++i) {
1752         print(o, (hev[i])) << " ";
1753       }
1754       return o << "} ";
1755     }
1756 
1757     template <typename stream_type>
print(stream_type & o,const std::vector<half_edge> & hev)1758     static stream_type& print (stream_type& o, const std::vector<half_edge>& hev)
1759     {
1760       o << "half edges: {";
1761       for(std::size_t i = 0; i < hev.size(); ++i) {
1762         print(o, (hev[i])) << " ";
1763       }
1764       return o << "} ";
1765     }
1766 
1767     template <typename stream_type>
print(stream_type & o,const vertex_data & v)1768     static stream_type& print (stream_type& o, const vertex_data& v)
1769     {
1770       return print(o << "vertex: <" << (v.first) << ", ", (v.second)) << "> ";
1771     }
1772 
1773     template <typename stream_type>
print(stream_type & o,const std::vector<vertex_data> & vv)1774     static stream_type& print (stream_type& o, const std::vector<vertex_data>& vv)
1775     {
1776       o << "vertices: {";
1777       for(std::size_t i = 0; i < vv.size(); ++i) {
1778         print(o, (vv[i])) << " ";
1779       }
1780       return o << "} ";
1781     }
1782 
1783 
1784 
1785     template <typename stream_type>
test_insertion(stream_type & stdcout)1786     static inline bool test_insertion(stream_type& stdcout) {
1787       property_merge si;
1788       rectangle_data<Unit> rect;
1789       xl(rect, 0);
1790       yl(rect, 1);
1791       xh(rect, 10);
1792       yh(rect, 11);
1793 
1794       si.insert(rect, 333);
1795       print(stdcout, si.pmd) << "\n";
1796 
1797       Point pts[4] = {Point(0, 0), Point(10,-3), Point(13, 8), Point(0, 0) };
1798       polygon_data<Unit> poly;
1799       property_merge si2;
1800       poly.set(pts, pts+3);
1801       si2.insert(poly, 444);
1802       si2.sort_property_merge_data();
1803       print(stdcout, si2.pmd) << "\n";
1804       property_merge si3;
1805       poly.set(pts, pts+4);
1806       si3.insert(poly, 444);
1807       si3.sort_property_merge_data();
1808       stdcout << (si2.pmd == si3.pmd) << "\n";
1809       std::reverse(pts, pts+4);
1810       property_merge si4;
1811       poly.set(pts, pts+4);
1812       si4.insert(poly, 444);
1813       si4.sort_property_merge_data();
1814       print(stdcout, si4.pmd) << "\n";
1815       stdcout << (si2.pmd == si4.pmd) << "\n";
1816       std::reverse(pts, pts+3);
1817       property_merge si5;
1818       poly.set(pts, pts+4);
1819       si5.insert(poly, 444);
1820       si5.sort_property_merge_data();
1821       stdcout << (si2.pmd == si5.pmd) << "\n";
1822 
1823       return true;
1824     }
1825 
1826     template <typename stream_type>
test_merge(stream_type & stdcout)1827     static inline bool test_merge(stream_type& stdcout) {
1828       property_merge si;
1829       rectangle_data<Unit> rect;
1830       xl(rect, 0);
1831       yl(rect, 1);
1832       xh(rect, 10);
1833       yh(rect, 11);
1834 
1835       si.insert(rect, 333);
1836       std::map<std::set<property_type>, polygon_set_data<Unit> > result;
1837       si.merge(result);
1838       print(stdcout, si.pmd) << "\n";
1839       polygon_set_data<Unit> psd = (*(result.begin())).second;
1840       std::vector<polygon_data<Unit> > polys;
1841       psd.get(polys);
1842       if(polys.size() != 1) {
1843         stdcout << "fail merge 1\n";
1844         return false;
1845       }
1846       stdcout << (polys[0]) << "\n";
1847       si.clear();
1848       std::vector<Point> pts;
1849       pts.push_back(Point(0, 0));
1850       pts.push_back(Point(10, -10));
1851       pts.push_back(Point(10, 10));
1852       polygon_data<Unit> poly;
1853       poly.set(pts.begin(), pts.end());
1854       si.insert(poly, 444);
1855       pts.clear();
1856       pts.push_back(Point(5, 0));
1857       pts.push_back(Point(-5, -10));
1858       pts.push_back(Point(-5, 10));
1859       poly.set(pts.begin(), pts.end());
1860       si.insert(poly, 444);
1861       result.clear();
1862       si.merge(result);
1863       print(stdcout, si.pmd) << "\n";
1864       psd = (*(result.begin())).second;
1865       stdcout << psd << "\n";
1866       polys.clear();
1867       psd.get(polys);
1868       if(polys.size() != 1) {
1869         stdcout << "fail merge 2\n";
1870         return false;
1871       }
1872       //Polygon { -4 -1, 3 3, -2 3 }
1873       //Polygon { 0 -4, -4 -2, -2 1 }
1874       si.clear();
1875       pts.clear();
1876       pts.push_back(Point(-4, -1));
1877       pts.push_back(Point(3, 3));
1878       pts.push_back(Point(-2, 3));
1879       poly.set(pts.begin(), pts.end());
1880       si.insert(poly, 444);
1881       pts.clear();
1882       pts.push_back(Point(0, -4));
1883       pts.push_back(Point(-4, -2));
1884       pts.push_back(Point(-2, 1));
1885       poly.set(pts.begin(), pts.end());
1886       si.insert(poly, 444);
1887       result.clear();
1888       si.merge(result);
1889       print(stdcout, si.pmd) << "\n";
1890       psd = (*(result.begin())).second;
1891       stdcout << psd << "\n";
1892       polys.clear();
1893       psd.get(polys);
1894       if(polys.size() != 1) {
1895         stdcout << "fail merge 3\n";
1896         return false;
1897       }
1898       stdcout << "Polygon { -2 2, -2 2, 1 4 } \n";
1899       stdcout << "Polygon { 2 4, 2 -4, -3 1 } \n";
1900       si.clear();
1901       pts.clear();
1902       pts.push_back(Point(-2, 2));
1903       pts.push_back(Point(-2, 2));
1904       pts.push_back(Point(1, 4));
1905       poly.set(pts.begin(), pts.end());
1906       si.insert(poly, 444);
1907       pts.clear();
1908       pts.push_back(Point(2, 4));
1909       pts.push_back(Point(2, -4));
1910       pts.push_back(Point(-3, 1));
1911       poly.set(pts.begin(), pts.end());
1912       si.insert(poly, 444);
1913       result.clear();
1914       si.merge(result);
1915       print(stdcout, si.pmd) << "\n";
1916       psd = (*(result.begin())).second;
1917       stdcout << psd << "\n";
1918       polys.clear();
1919       psd.get(polys);
1920       if(polys.size() != 1) {
1921         stdcout << "fail merge 4\n";
1922         return false;
1923       }
1924       stdcout << (polys[0]) << "\n";
1925       stdcout << "Polygon { -4 0, -2 -3, 3 -4 } \n";
1926       stdcout << "Polygon { -1 1, 1 -2, -4 -3 } \n";
1927       si.clear();
1928       pts.clear();
1929       pts.push_back(Point(-4, 0));
1930       pts.push_back(Point(-2, -3));
1931       pts.push_back(Point(3, -4));
1932       poly.set(pts.begin(), pts.end());
1933       si.insert(poly, 444);
1934       pts.clear();
1935       pts.push_back(Point(-1, 1));
1936       pts.push_back(Point(1, -2));
1937       pts.push_back(Point(-4, -3));
1938       poly.set(pts.begin(), pts.end());
1939       si.insert(poly, 444);
1940       result.clear();
1941       si.merge(result);
1942       print(stdcout, si.pmd) << "\n";
1943       psd = (*(result.begin())).second;
1944       stdcout << psd << "\n";
1945       polys.clear();
1946       psd.get(polys);
1947       if(polys.size() != 1) {
1948         stdcout << "fail merge 5\n";
1949         return false;
1950       }
1951       stdcout << "Polygon { 2 2, -2 0, 0 1 }  \n";
1952       stdcout << "Polygon { 4 -2, 3 -1, 2 3 }  \n";
1953       si.clear();
1954       pts.clear();
1955       pts.push_back(Point(2, 2));
1956       pts.push_back(Point(-2, 0));
1957       pts.push_back(Point(0, 1));
1958       poly.set(pts.begin(), pts.end());
1959       si.insert(poly, 444);
1960       pts.clear();
1961       pts.push_back(Point(4, -2));
1962       pts.push_back(Point(3, -1));
1963       pts.push_back(Point(2, 3));
1964       poly.set(pts.begin(), pts.end());
1965       si.insert(poly, 444);
1966       result.clear();
1967       si.merge(result);
1968       print(stdcout, si.pmd) << "\n";
1969       if(!result.empty()) {
1970         psd = (*(result.begin())).second;
1971         stdcout << psd << "\n";
1972         polys.clear();
1973         psd.get(polys);
1974         if(polys.size() != 1) {
1975           stdcout << "fail merge 6\n";
1976           return false;
1977         }
1978         stdcout << (polys[0]) << "\n";
1979       }
1980       stdcout << "Polygon { 0 2, 3 -1, 4 1 }  \n";
1981       stdcout << "Polygon { -4 3, 3 3, 4 2 }  \n";
1982       si.clear();
1983       pts.clear();
1984       pts.push_back(Point(0, 2));
1985       pts.push_back(Point(3, -1));
1986       pts.push_back(Point(4, 1));
1987       poly.set(pts.begin(), pts.end());
1988       si.insert(poly, 444);
1989       pts.clear();
1990       pts.push_back(Point(-4, 3));
1991       pts.push_back(Point(3, 3));
1992       pts.push_back(Point(4, 2));
1993       poly.set(pts.begin(), pts.end());
1994       si.insert(poly, 444);
1995       result.clear();
1996       si.merge(result);
1997       print(stdcout, si.pmd) << "\n";
1998       if(!result.empty()) {
1999         psd = (*(result.begin())).second;
2000         stdcout << psd << "\n";
2001         polys.clear();
2002         psd.get(polys);
2003         if(polys.size() == 0) {
2004           stdcout << "fail merge 7\n";
2005           return false;
2006         }
2007         stdcout << (polys[0]) << "\n";
2008       }
2009 stdcout << "Polygon { 1 -2, -1 4, 3 -2 }   \n";
2010 stdcout << "Polygon { 0 -3, 3 1, -3 -4 }   \n";
2011       si.clear();
2012       pts.clear();
2013       pts.push_back(Point(1, -2));
2014       pts.push_back(Point(-1, 4));
2015       pts.push_back(Point(3, -2));
2016       poly.set(pts.begin(), pts.end());
2017       si.insert(poly, 444);
2018       pts.clear();
2019       pts.push_back(Point(0, -3));
2020       pts.push_back(Point(3, 1));
2021       pts.push_back(Point(-3, -4));
2022       poly.set(pts.begin(), pts.end());
2023       si.insert(poly, 444);
2024       result.clear();
2025       si.merge(result);
2026       print(stdcout, si.pmd) << "\n";
2027       if(!result.empty()) {
2028         psd = (*(result.begin())).second;
2029         stdcout << psd << "\n";
2030         polys.clear();
2031         psd.get(polys);
2032         if(polys.size() == 0) {
2033           stdcout << "fail merge 8\n";
2034           return false;
2035         }
2036         stdcout << (polys[0]) << "\n";
2037       }
2038 stdcout << "Polygon { 2 2, 3 0, -3 4 }  \n";
2039 stdcout << "Polygon { -2 -2, 0 0, -1 -1 }  \n";
2040       si.clear();
2041       pts.clear();
2042       pts.push_back(Point(2, 2));
2043       pts.push_back(Point(3, 0));
2044       pts.push_back(Point(-3, 4));
2045       poly.set(pts.begin(), pts.end());
2046       si.insert(poly, 444);
2047       pts.clear();
2048       pts.push_back(Point(-2, -2));
2049       pts.push_back(Point(0, 0));
2050       pts.push_back(Point(-1, -1));
2051       poly.set(pts.begin(), pts.end());
2052       si.insert(poly, 444);
2053       result.clear();
2054       si.merge(result);
2055       print(stdcout, si.pmd) << "\n";
2056       if(!result.empty()) {
2057         psd = (*(result.begin())).second;
2058         stdcout << psd << "\n";
2059         polys.clear();
2060         psd.get(polys);
2061         if(polys.size() == 0) {
2062           stdcout << "fail merge 9\n";
2063           return false;
2064         }
2065         stdcout << (polys[0]) << "\n";
2066       }
2067       si.clear();
2068       pts.clear();
2069       //5624841,17616200,75000,9125000
2070       //pts.push_back(Point(5624841,75000));
2071       //pts.push_back(Point(5624841,9125000));
2072       //pts.push_back(Point(17616200,9125000));
2073       //pts.push_back(Point(17616200,75000));
2074 pts.push_back(Point(12262940, 6652520 )); pts.push_back(Point(12125750, 6652520 )); pts.push_back(Point(12121272, 6652961 )); pts.push_back(Point(12112981, 6656396 )); pts.push_back(Point(12106636, 6662741 )); pts.push_back(Point(12103201, 6671032 )); pts.push_back(Point(12103201, 6680007 )); pts.push_back(Point(12106636, 6688298 ));
2075 pts.push_back(Point(12109500, 6691780 )); pts.push_back(Point(12748600, 7330890 )); pts.push_back(Point(15762600, 7330890 )); pts.push_back(Point(15904620, 7472900 )); pts.push_back(Point(15909200, 7473030 )); pts.push_back(Point(15935830, 7476006 )); pts.push_back(Point(15992796, 7499602 )); pts.push_back(Point(16036397, 7543203 ));
2076 pts.push_back(Point(16059993, 7600169 )); pts.push_back(Point(16059993, 7661830 )); pts.push_back(Point(16036397, 7718796 )); pts.push_back(Point(15992796, 7762397 )); pts.push_back(Point(15935830, 7785993 )); pts.push_back(Point(15874169, 7785993 )); pts.push_back(Point(15817203, 7762397 )); pts.push_back(Point(15773602, 7718796 ));
2077 pts.push_back(Point(15750006, 7661830 )); pts.push_back(Point(15747030, 7635200 )); pts.push_back(Point(15746900, 7630620 )); pts.push_back(Point(15670220, 7553930 )); pts.push_back(Point(14872950, 7553930 )); pts.push_back(Point(14872950, 7626170 ));
2078 pts.push_back(Point(14869973, 7661280 )); pts.push_back(Point(14846377, 7718246 )); pts.push_back(Point(14802776, 7761847 )); pts.push_back(Point(14745810, 7785443 )); pts.push_back(Point(14684149, 7785443 )); pts.push_back(Point(14627183, 7761847 )); pts.push_back(Point(14583582, 7718246 ));
2079 pts.push_back(Point(14559986, 7661280 )); pts.push_back(Point(14557070, 7636660 )); pts.push_back(Point(14556670, 7625570 )); pts.push_back(Point(13703330, 7625570 )); pts.push_back(Point(13702930, 7636660 )); pts.push_back(Point(13699993, 7661830 )); pts.push_back(Point(13676397, 7718796 ));
2080 pts.push_back(Point(13632796, 7762397 )); pts.push_back(Point(13575830, 7785993 )); pts.push_back(Point(13514169, 7785993 )); pts.push_back(Point(13457203, 7762397 )); pts.push_back(Point(13436270, 7745670 )); pts.push_back(Point(13432940, 7742520 )); pts.push_back(Point(12963760, 7742520 ));
2081 pts.push_back(Point(12959272, 7742961 )); pts.push_back(Point(12950981, 7746396 )); pts.push_back(Point(12944636, 7752741 )); pts.push_back(Point(12941201, 7761032 )); pts.push_back(Point(12941201, 7770007 )); pts.push_back(Point(12944636, 7778298 )); pts.push_back(Point(12947490, 7781780 ));
2082 pts.push_back(Point(13425330, 8259620 )); pts.push_back(Point(15601330, 8259620 )); pts.push_back(Point(15904620, 8562900 )); pts.push_back(Point(15909200, 8563030 )); pts.push_back(Point(15935830, 8566006 )); pts.push_back(Point(15992796, 8589602 )); pts.push_back(Point(16036397, 8633203 ));
2083 pts.push_back(Point(16059993, 8690169 )); pts.push_back(Point(16059993, 8751830 )); pts.push_back(Point(16036397, 8808796 )); pts.push_back(Point(15992796, 8852397 )); pts.push_back(Point(15935830, 8875993 )); pts.push_back(Point(15874169, 8875993 )); pts.push_back(Point(15817203, 8852397 )); pts.push_back(Point(15773602, 8808796 ));
2084 pts.push_back(Point(15750006, 8751830 )); pts.push_back(Point(15747030, 8725200 )); pts.push_back(Point(15746900, 8720620 )); pts.push_back(Point(15508950, 8482660 )); pts.push_back(Point(14689890, 8482660 )); pts.push_back(Point(14685412, 8483101 )); pts.push_back(Point(14677121, 8486536 ));
2085 pts.push_back(Point(14670776, 8492881 )); pts.push_back(Point(14667341, 8501172 )); pts.push_back(Point(14667341, 8510147 )); pts.push_back(Point(14670776, 8518438 )); pts.push_back(Point(14673630, 8521920 )); pts.push_back(Point(14714620, 8562900 )); pts.push_back(Point(14719200, 8563030 )); pts.push_back(Point(14745830, 8566006 ));
2086 pts.push_back(Point(14802796, 8589602 )); pts.push_back(Point(14846397, 8633203 )); pts.push_back(Point(14869993, 8690169 )); pts.push_back(Point(14869993, 8751830 )); pts.push_back(Point(14846397, 8808796 )); pts.push_back(Point(14802796, 8852397 )); pts.push_back(Point(14745830, 8875993 )); pts.push_back(Point(14684169, 8875993 ));
2087 pts.push_back(Point(14627203, 8852397 )); pts.push_back(Point(14583602, 8808796 )); pts.push_back(Point(14560006, 8751830 )); pts.push_back(Point(14557030, 8725200 )); pts.push_back(Point(14556900, 8720620 )); pts.push_back(Point(14408270, 8571980 )); pts.push_back(Point(13696320, 8571980 )); pts.push_back(Point(13696320, 8675520 ));
2088 pts.push_back(Point(13699963, 8690161 )); pts.push_back(Point(13699963, 8751818 )); pts.push_back(Point(13676368, 8808781 )); pts.push_back(Point(13632771, 8852378 )); pts.push_back(Point(13575808, 8875973 )); pts.push_back(Point(13514151, 8875973 )); pts.push_back(Point(13457188, 8852378 )); pts.push_back(Point(13436270, 8835670 )); pts.push_back(Point(13432940, 8832520 ));
2089 pts.push_back(Point(13281760, 8832520 )); pts.push_back(Point(13277272, 8832961 )); pts.push_back(Point(13268981, 8836396 )); pts.push_back(Point(13262636, 8842741 )); pts.push_back(Point(13259201, 8851032 )); pts.push_back(Point(13259201, 8860007 )); pts.push_back(Point(13262636, 8868298 )); pts.push_back(Point(13265500, 8871780 ));
2090 pts.push_back(Point(13518710, 9125000 )); pts.push_back(Point(16270720, 9125000 )); pts.push_back(Point(16270720, 8939590 )); pts.push_back(Point(17120780, 8939590 )); pts.push_back(Point(17120780, 9125000 )); pts.push_back(Point(17616200, 9125000 )); pts.push_back(Point(17616200,   75000 )); pts.push_back(Point(16024790,   75000 ));
2091 pts.push_back(Point(16021460,   80700 )); pts.push_back(Point(16016397,   88796 )); pts.push_back(Point(15972796,  132397 )); pts.push_back(Point(15915830,  155993 )); pts.push_back(Point(15908730,  157240 )); pts.push_back(Point(15905000,  157800 )); pts.push_back(Point(15516800,  546000 )); pts.push_back(Point(15905000,  934200 ));
2092 pts.push_back(Point(15908730,  934760 )); pts.push_back(Point(15915830,  936006 )); pts.push_back(Point(15972796,  959602 )); pts.push_back(Point(16016397, 1003203 )); pts.push_back(Point(16039993, 1060169 )); pts.push_back(Point(16039993, 1121830 )); pts.push_back(Point(16016397, 1178796 )); pts.push_back(Point(15972796, 1222397 ));
2093 pts.push_back(Point(15915830, 1245993 )); pts.push_back(Point(15854169, 1245993 )); pts.push_back(Point(15797203, 1222397 )); pts.push_back(Point(15753602, 1178796 )); pts.push_back(Point(15730006, 1121830 )); pts.push_back(Point(15728760, 1114730 )); pts.push_back(Point(15728200, 1111000 )); pts.push_back(Point(15363500,  746300 ));
2094 pts.push_back(Point(14602620,  746300 )); pts.push_back(Point(14598142,  746741 )); pts.push_back(Point(14589851,  750176 )); pts.push_back(Point(14583506,  756521 )); pts.push_back(Point(14580071,  764812 )); pts.push_back(Point(14580071,  773787 )); pts.push_back(Point(14583506,  782078 )); pts.push_back(Point(14586360,  785560 ));
2095 pts.push_back(Point(14586370,  785560 )); pts.push_back(Point(14735000,  934200 )); pts.push_back(Point(14738730,  934760 )); pts.push_back(Point(14745830,  936006 )); pts.push_back(Point(14802796,  959602 )); pts.push_back(Point(14846397, 1003203 )); pts.push_back(Point(14869993, 1060169 ));
2096 pts.push_back(Point(14870450, 1062550 )); pts.push_back(Point(14872170, 1071980 )); pts.push_back(Point(14972780, 1071980 )); pts.push_back(Point(15925000, 2024200 )); pts.push_back(Point(15928730, 2024760 )); pts.push_back(Point(15935830, 2026006 )); pts.push_back(Point(15992796, 2049602 ));
2097 pts.push_back(Point(16036397, 2093203 )); pts.push_back(Point(16059993, 2150169 )); pts.push_back(Point(16059993, 2211830 )); pts.push_back(Point(16036397, 2268796 )); pts.push_back(Point(15992796, 2312397 )); pts.push_back(Point(15935830, 2335993 )); pts.push_back(Point(15874169, 2335993 ));
2098 pts.push_back(Point(15817203, 2312397 )); pts.push_back(Point(15773602, 2268796 )); pts.push_back(Point(15750006, 2211830 )); pts.push_back(Point(15748760, 2204730 )); pts.push_back(Point(15748200, 2201000 )); pts.push_back(Point(14869220, 1322020 )); pts.push_back(Point(14088350, 1322020 ));
2099 pts.push_back(Point(14083862, 1322461 )); pts.push_back(Point(14075571, 1325896 )); pts.push_back(Point(14069226, 1332241 )); pts.push_back(Point(14065791, 1340532 )); pts.push_back(Point(14065791, 1349507 )); pts.push_back(Point(14069226, 1357798 )); pts.push_back(Point(14072080, 1361280 ));
2100 pts.push_back(Point(14072090, 1361280 )); pts.push_back(Point(14735000, 2024200 )); pts.push_back(Point(14738730, 2024760 )); pts.push_back(Point(14745830, 2026006 )); pts.push_back(Point(14802796, 2049602 )); pts.push_back(Point(14846397, 2093203 )); pts.push_back(Point(14869993, 2150169 ));
2101 pts.push_back(Point(14869993, 2211830 )); pts.push_back(Point(14846397, 2268796 )); pts.push_back(Point(14802796, 2312397 )); pts.push_back(Point(14745830, 2335993 )); pts.push_back(Point(14684169, 2335993 )); pts.push_back(Point(14627203, 2312397 )); pts.push_back(Point(14583602, 2268796 )); pts.push_back(Point(14560006, 2211830 ));
2102 pts.push_back(Point(14558760, 2204730 )); pts.push_back(Point(14558200, 2201000 )); pts.push_back(Point(13752220, 1395020 )); pts.push_back(Point(12991340, 1395020 )); pts.push_back(Point(12986862, 1395461 )); pts.push_back(Point(12978571, 1398896 )); pts.push_back(Point(12972226, 1405241 ));
2103 pts.push_back(Point(12968791, 1413532 )); pts.push_back(Point(12968791, 1422507 )); pts.push_back(Point(12972226, 1430798 )); pts.push_back(Point(12975080, 1434280 )); pts.push_back(Point(12975090, 1434280 )); pts.push_back(Point(13565000, 2024200 )); pts.push_back(Point(13568730, 2024760 )); pts.push_back(Point(13575830, 2026006 ));
2104 pts.push_back(Point(13632796, 2049602 )); pts.push_back(Point(13676397, 2093203 )); pts.push_back(Point(13699993, 2150169 )); pts.push_back(Point(13699993, 2211830 )); pts.push_back(Point(13676397, 2268796 )); pts.push_back(Point(13632796, 2312397 )); pts.push_back(Point(13575830, 2335993 ));
2105 pts.push_back(Point(13514169, 2335993 )); pts.push_back(Point(13457203, 2312397 )); pts.push_back(Point(13413602, 2268796 )); pts.push_back(Point(13390006, 2211830 )); pts.push_back(Point(13388760, 2204730 )); pts.push_back(Point(13388200, 2201000 )); pts.push_back(Point(12655220, 1468020 ));
2106 pts.push_back(Point(11894340, 1468020 )); pts.push_back(Point(11889862, 1468461 )); pts.push_back(Point(11881571, 1471896 )); pts.push_back(Point(11875226, 1478241 )); pts.push_back(Point(11871791, 1486532 )); pts.push_back(Point(11871791, 1495507 ));
2107 pts.push_back(Point(11875226, 1503798 )); pts.push_back(Point(11878090, 1507280 )); pts.push_back(Point(12395000, 2024200 )); pts.push_back(Point(12398730, 2024760 )); pts.push_back(Point(12405830, 2026006 )); pts.push_back(Point(12462796, 2049602 )); pts.push_back(Point(12506397, 2093203 ));
2108 pts.push_back(Point(12529993, 2150169 )); pts.push_back(Point(12529993, 2211830 )); pts.push_back(Point(12506397, 2268796 )); pts.push_back(Point(12462796, 2312397 )); pts.push_back(Point(12405830, 2335993 )); pts.push_back(Point(12344169, 2335993 ));
2109 pts.push_back(Point(12287203, 2312397 )); pts.push_back(Point(12243602, 2268796 )); pts.push_back(Point(12220006, 2211830 )); pts.push_back(Point(12218760, 2204730 )); pts.push_back(Point(12218200, 2201000 )); pts.push_back(Point(11558220, 1541020 ));
2110 pts.push_back(Point(10797340, 1541020 )); pts.push_back(Point(10792862, 1541461 )); pts.push_back(Point(10784571, 1544896 )); pts.push_back(Point(10778226, 1551241 )); pts.push_back(Point(10774791, 1559532 )); pts.push_back(Point(10774791, 1568507 )); pts.push_back(Point(10778226, 1576798 )); pts.push_back(Point(10781080, 1580280 ));
2111 pts.push_back(Point(10781090, 1580280 )); pts.push_back(Point(11225000, 2024200 )); pts.push_back(Point(11228730, 2024760 )); pts.push_back(Point(11235830, 2026006 )); pts.push_back(Point(11292796, 2049602 )); pts.push_back(Point(11336397, 2093203 )); pts.push_back(Point(11359993, 2150169 ));
2112 pts.push_back(Point(11359993, 2211830 )); pts.push_back(Point(11336397, 2268796 )); pts.push_back(Point(11292796, 2312397 )); pts.push_back(Point(11235830, 2335993 )); pts.push_back(Point(11174169, 2335993 )); pts.push_back(Point(11117203, 2312397 )); pts.push_back(Point(11073602, 2268796 )); pts.push_back(Point(11050006, 2211830 ));
2113 pts.push_back(Point(11048760, 2204730 )); pts.push_back(Point(11048200, 2201000 )); pts.push_back(Point(10461220, 1614020 )); pts.push_back(Point( 5647400, 1614020 )); pts.push_back(Point( 5642912, 1614461 ));
2114 pts.push_back(Point( 5634621, 1617896 )); pts.push_back(Point( 5628276, 1624241 )); pts.push_back(Point( 5624841, 1632532 )); pts.push_back(Point( 5624841, 1641507 )); pts.push_back(Point( 5628276, 1649798 )); pts.push_back(Point( 5631130, 1653280 ));
2115 pts.push_back(Point( 5688490, 1710640 )); pts.push_back(Point( 9722350, 1710640 )); pts.push_back(Point(10034620, 2022900 )); pts.push_back(Point(10039200, 2023030 )); pts.push_back(Point(10065830, 2026006 )); pts.push_back(Point(10122796, 2049602 ));
2116 pts.push_back(Point(10166397, 2093203 )); pts.push_back(Point(10189993, 2150169 )); pts.push_back(Point(10189993, 2211830 )); pts.push_back(Point(10166397, 2268796 )); pts.push_back(Point(10158620, 2279450 )); pts.push_back(Point(10158620, 2404900 )); pts.push_back(Point(10548950, 2795240 ));
2117 pts.push_back(Point(15586950, 2795240 )); pts.push_back(Point(15904620, 3112900 )); pts.push_back(Point(15909200, 3113030 )); pts.push_back(Point(15935830, 3116006 )); pts.push_back(Point(15992796, 3139602 )); pts.push_back(Point(16036397, 3183203 )); pts.push_back(Point(16059993, 3240169 )); pts.push_back(Point(16059993, 3301830 ));
2118 pts.push_back(Point(16036397, 3358796 )); pts.push_back(Point(15992796, 3402397 )); pts.push_back(Point(15935830, 3425993 )); pts.push_back(Point(15874169, 3425993 )); pts.push_back(Point(15817203, 3402397 )); pts.push_back(Point(15773602, 3358796 )); pts.push_back(Point(15750006, 3301830 )); pts.push_back(Point(15747030, 3275200 ));
2119 pts.push_back(Point(15746900, 3270620 )); pts.push_back(Point(15494570, 3018280 )); pts.push_back(Point(14675510, 3018280 )); pts.push_back(Point(14671032, 3018721 )); pts.push_back(Point(14662741, 3022156 )); pts.push_back(Point(14656396, 3028501 )); pts.push_back(Point(14652961, 3036792 ));
2120 pts.push_back(Point(14652961, 3045767 )); pts.push_back(Point(14656396, 3054058 )); pts.push_back(Point(14659260, 3057540 )); pts.push_back(Point(14714620, 3112900 )); pts.push_back(Point(14719200, 3113030 )); pts.push_back(Point(14745830, 3116006 )); pts.push_back(Point(14802796, 3139602 ));
2121 pts.push_back(Point(14846397, 3183203 )); pts.push_back(Point(14869993, 3240169 )); pts.push_back(Point(14869993, 3301830 )); pts.push_back(Point(14846397, 3358796 )); pts.push_back(Point(14802796, 3402397 )); pts.push_back(Point(14745830, 3425993 )); pts.push_back(Point(14684169, 3425993 )); pts.push_back(Point(14627203, 3402397 ));
2122 pts.push_back(Point(14583602, 3358796 )); pts.push_back(Point(14560006, 3301830 )); pts.push_back(Point(14557030, 3275200 )); pts.push_back(Point(14556900, 3270620 )); pts.push_back(Point(14370700, 3084410 )); pts.push_back(Point(13702830, 3084410 )); pts.push_back(Point(13702830, 3263160 ));
2123 pts.push_back(Point(13700003, 3302210 )); pts.push_back(Point(13676407, 3359176 )); pts.push_back(Point(13632806, 3402777 )); pts.push_back(Point(13575840, 3426373 )); pts.push_back(Point(13514179, 3426373 )); pts.push_back(Point(13457213, 3402777 )); pts.push_back(Point(13413612, 3359176 ));
2124 pts.push_back(Point(13390016, 3302210 )); pts.push_back(Point(13387030, 3275200 )); pts.push_back(Point(13386900, 3270620 )); pts.push_back(Point(13266840, 3150550 )); pts.push_back(Point(12532920, 3150550 )); pts.push_back(Point(12532920, 3264990 )); pts.push_back(Point(12529993, 3301820 ));
2125 pts.push_back(Point(12506397, 3358786 )); pts.push_back(Point(12462796, 3402387 )); pts.push_back(Point(12405830, 3425983 )); pts.push_back(Point(12344169, 3425983 )); pts.push_back(Point(12287203, 3402387 )); pts.push_back(Point(12243602, 3358786 )); pts.push_back(Point(12220006, 3301820 )); pts.push_back(Point(12217030, 3275200 ));
2126 pts.push_back(Point(12216900, 3270620 )); pts.push_back(Point(12157460, 3211170 )); pts.push_back(Point(11362030, 3211170 )); pts.push_back(Point(11360250, 3220520 )); pts.push_back(Point(11359993, 3221830 )); pts.push_back(Point(11336397, 3278796 ));
2127 pts.push_back(Point(11292796, 3322397 )); pts.push_back(Point(11235830, 3345993 )); pts.push_back(Point(11174169, 3345993 )); pts.push_back(Point(11117203, 3322397 )); pts.push_back(Point(11096270, 3305670 )); pts.push_back(Point(11092940, 3302520 )); pts.push_back(Point(10680760, 3302520 ));
2128 pts.push_back(Point(10676272, 3302961 )); pts.push_back(Point(10667981, 3306396 )); pts.push_back(Point(10661636, 3312741 )); pts.push_back(Point(10658201, 3321032 )); pts.push_back(Point(10658201, 3330007 )); pts.push_back(Point(10661636, 3338298 )); pts.push_back(Point(10664500, 3341780 ));
2129 pts.push_back(Point(11264260, 3941550 )); pts.push_back(Point(15643260, 3941550 )); pts.push_back(Point(15904620, 4202900 )); pts.push_back(Point(15909200, 4203030 )); pts.push_back(Point(15935830, 4206006 )); pts.push_back(Point(15992796, 4229602 ));
2130 pts.push_back(Point(16036397, 4273203 )); pts.push_back(Point(16059993, 4330169 )); pts.push_back(Point(16059993, 4391830 )); pts.push_back(Point(16036397, 4448796 )); pts.push_back(Point(15992796, 4492397 ));
2131 pts.push_back(Point(15935830, 4515993 )); pts.push_back(Point(15874169, 4515993 )); pts.push_back(Point(15817203, 4492397 )); pts.push_back(Point(15773602, 4448796 )); pts.push_back(Point(15750006, 4391830 )); pts.push_back(Point(15747030, 4365200 )); pts.push_back(Point(15746900, 4360620 ));
2132 pts.push_back(Point(15550880, 4164590 )); pts.push_back(Point(14825070, 4164590 )); pts.push_back(Point(14825070, 4247610 )); pts.push_back(Point(14846397, 4273213 )); pts.push_back(Point(14869993, 4330179 )); pts.push_back(Point(14869993, 4391840 )); pts.push_back(Point(14846397, 4448806 ));
2133 pts.push_back(Point(14802796, 4492407 )); pts.push_back(Point(14745830, 4516003 )); pts.push_back(Point(14684169, 4516003 )); pts.push_back(Point(14627203, 4492407 )); pts.push_back(Point(14583602, 4448806 )); pts.push_back(Point(14560006, 4391840 )); pts.push_back(Point(14557030, 4365200 ));
2134 pts.push_back(Point(14556900, 4360620 )); pts.push_back(Point(14432520, 4236230 )); pts.push_back(Point(13702830, 4236230 )); pts.push_back(Point(13702830, 4352930 )); pts.push_back(Point(13699993, 4391750 )); pts.push_back(Point(13676397, 4448716 ));
2135 pts.push_back(Point(13632796, 4492317 )); pts.push_back(Point(13575830, 4515913 )); pts.push_back(Point(13514169, 4515913 )); pts.push_back(Point(13457203, 4492317 )); pts.push_back(Point(13413602, 4448716 )); pts.push_back(Point(13390006, 4391750 )); pts.push_back(Point(13387030, 4365200 ));
2136 pts.push_back(Point(13386900, 4360620 )); pts.push_back(Point(13334170, 4307880 )); pts.push_back(Point(12532990, 4307880 )); pts.push_back(Point(12532990, 4357550 )); pts.push_back(Point(12529993, 4391760 )); pts.push_back(Point(12506397, 4448726 )); pts.push_back(Point(12462796, 4492327 ));
2137 pts.push_back(Point(12405830, 4515923 )); pts.push_back(Point(12344169, 4515923 )); pts.push_back(Point(12287203, 4492327 )); pts.push_back(Point(12243602, 4448726 )); pts.push_back(Point(12220006, 4391760 )); pts.push_back(Point(12217970, 4378710 )); pts.push_back(Point(12216810, 4368500 ));
2138 pts.push_back(Point(11363190, 4368500 )); pts.push_back(Point(11362030, 4378710 )); pts.push_back(Point(11359983, 4391828 )); pts.push_back(Point(11336388, 4448791 )); pts.push_back(Point(11292791, 4492388 )); pts.push_back(Point(11235828, 4515983 )); pts.push_back(Point(11174171, 4515983 )); pts.push_back(Point(11117208, 4492388 ));
2139 pts.push_back(Point(11096270, 4475670 )); pts.push_back(Point(11092940, 4472520 )); pts.push_back(Point(11057750, 4472520 )); pts.push_back(Point(11053272, 4472961 )); pts.push_back(Point(11044981, 4476396 )); pts.push_back(Point(11038636, 4482741 )); pts.push_back(Point(11035201, 4491032 ));
2140 pts.push_back(Point(11035201, 4500007 )); pts.push_back(Point(11038636, 4508298 )); pts.push_back(Point(11041490, 4511780 )); pts.push_back(Point(11573490, 5043780 )); pts.push_back(Point(15655490, 5043780 )); pts.push_back(Point(15904620, 5292900 ));
2141 pts.push_back(Point(15909200, 5293030 )); pts.push_back(Point(15935830, 5296006 )); pts.push_back(Point(15992796, 5319602 )); pts.push_back(Point(16036397, 5363203 )); pts.push_back(Point(16059993, 5420169 )); pts.push_back(Point(16059993, 5481830 )); pts.push_back(Point(16036397, 5538796 ));
2142 pts.push_back(Point(15992796, 5582397 )); pts.push_back(Point(15935830, 5605993 )); pts.push_back(Point(15874169, 5605993 )); pts.push_back(Point(15817203, 5582397 )); pts.push_back(Point(15773602, 5538796 )); pts.push_back(Point(15750006, 5481830 )); pts.push_back(Point(15747030, 5455200 ));
2143 pts.push_back(Point(15746900, 5450620 )); pts.push_back(Point(15563110, 5266820 )); pts.push_back(Point(14857380, 5266820 )); pts.push_back(Point(14857380, 5382430 )); pts.push_back(Point(14869993, 5420179 )); pts.push_back(Point(14869993, 5481840 )); pts.push_back(Point(14846397, 5538806 )); pts.push_back(Point(14802796, 5582407 ));
2144 pts.push_back(Point(14745830, 5606003 )); pts.push_back(Point(14684169, 5606003 )); pts.push_back(Point(14627203, 5582407 )); pts.push_back(Point(14583602, 5538806 )); pts.push_back(Point(14560006, 5481840 )); pts.push_back(Point(14557030, 5455200 )); pts.push_back(Point(14556900, 5450620 )); pts.push_back(Point(14444750, 5338460 ));
2145 pts.push_back(Point(13702890, 5338460 )); pts.push_back(Point(13702890, 5364400 )); pts.push_back(Point(13699993, 5401800 )); pts.push_back(Point(13676397, 5458766 )); pts.push_back(Point(13632796, 5502367 )); pts.push_back(Point(13575830, 5525963 )); pts.push_back(Point(13514169, 5525963 )); pts.push_back(Point(13457203, 5502367 ));
2146 pts.push_back(Point(13413602, 5458766 )); pts.push_back(Point(13390006, 5401800 )); pts.push_back(Point(13389230, 5397620 )); pts.push_back(Point(13387590, 5388060 )); pts.push_back(Point(12532960, 5388060 )); pts.push_back(Point(12532960, 5446220 )); pts.push_back(Point(12529993, 5481820 ));
2147 pts.push_back(Point(12506397, 5538786 )); pts.push_back(Point(12462796, 5582387 )); pts.push_back(Point(12405830, 5605983 )); pts.push_back(Point(12344169, 5605983 )); pts.push_back(Point(12287203, 5582387 )); pts.push_back(Point(12266270, 5565670 )); pts.push_back(Point(12262940, 5562520 )); pts.push_back(Point(11737750, 5562520 ));
2148 pts.push_back(Point(11733272, 5562961 )); pts.push_back(Point(11724981, 5566396 )); pts.push_back(Point(11718636, 5572741 )); pts.push_back(Point(11715201, 5581032 )); pts.push_back(Point(11715201, 5590007 )); pts.push_back(Point(11718636, 5598298 )); pts.push_back(Point(11721500, 5601780 ));
2149 pts.push_back(Point(12287760, 6168050 )); pts.push_back(Point(15689760, 6168050 )); pts.push_back(Point(15904620, 6382900 )); pts.push_back(Point(15909200, 6383030 )); pts.push_back(Point(15935830, 6386006 )); pts.push_back(Point(15992796, 6409602 ));
2150 pts.push_back(Point(16036397, 6453203 )); pts.push_back(Point(16059993, 6510169 )); pts.push_back(Point(16059993, 6571830 )); pts.push_back(Point(16036397, 6628796 )); pts.push_back(Point(15992796, 6672397 )); pts.push_back(Point(15935830, 6695993 )); pts.push_back(Point(15874169, 6695993 ));
2151 pts.push_back(Point(15817203, 6672397 )); pts.push_back(Point(15773602, 6628796 )); pts.push_back(Point(15750006, 6571830 )); pts.push_back(Point(15747030, 6545200 )); pts.push_back(Point(15746900, 6540620 )); pts.push_back(Point(15597380, 6391090 )); pts.push_back(Point(14858060, 6391090 ));
2152 pts.push_back(Point(14858060, 6473860 )); pts.push_back(Point(14869993, 6510179 )); pts.push_back(Point(14869993, 6571840 )); pts.push_back(Point(14846397, 6628806 )); pts.push_back(Point(14802796, 6672407 )); pts.push_back(Point(14745830, 6696003 )); pts.push_back(Point(14684169, 6696003 ));
2153 pts.push_back(Point(14627203, 6672407 )); pts.push_back(Point(14583602, 6628806 )); pts.push_back(Point(14560006, 6571840 )); pts.push_back(Point(14557030, 6545200 )); pts.push_back(Point(14556900, 6540620 )); pts.push_back(Point(14479020, 6462730 ));
2154 pts.push_back(Point(13702990, 6462730 )); pts.push_back(Point(13702990, 6537170 )); pts.push_back(Point(13700003, 6571840 )); pts.push_back(Point(13676407, 6628806 )); pts.push_back(Point(13632806, 6672407 )); pts.push_back(Point(13575840, 6696003 ));
2155 pts.push_back(Point(13514179, 6696003 )); pts.push_back(Point(13457213, 6672407 )); pts.push_back(Point(13413612, 6628806 )); pts.push_back(Point(13390016, 6571840 )); pts.push_back(Point(13387040, 6545550 )); pts.push_back(Point(13386710, 6534380 ));
2156 pts.push_back(Point(12533290, 6534380 )); pts.push_back(Point(12532960, 6545550 )); pts.push_back(Point(12529983, 6571828 )); pts.push_back(Point(12506388, 6628791 )); pts.push_back(Point(12462791, 6672388 )); pts.push_back(Point(12405828, 6695983 ));
2157 pts.push_back(Point(12344171, 6695983 )); pts.push_back(Point(12287208, 6672388 )); pts.push_back(Point(12266270, 6655670 ));
2158       poly.set(pts.begin(), pts.end());
2159       si.insert(poly, 444);
2160       result.clear();
2161       si.merge(result);
2162       si.verify1();
2163       print(stdcout, si.pmd) << "\n";
2164       if(!result.empty()) {
2165         psd = (*(result.begin())).second;
2166         stdcout << psd << "\n";
2167         std::vector<Point> outpts;
2168         for(typename polygon_set_data<Unit>::iterator_type itr = psd.begin();
2169             itr != psd.end(); ++itr) {
2170           outpts.push_back((*itr).first.first);
2171           outpts.push_back((*itr).first.second);
2172         }
2173         polygon_sort(outpts.begin(), outpts.end());
2174         for(std::size_t i = 0; i < outpts.size(); i+=2) {
2175           if(outpts[i] != outpts[i+1]) {
2176             stdcout << "Polygon set not a closed figure\n";
2177             stdcout << i << "\n";
2178             stdcout << outpts[i] << " " << outpts[i+1] << "\n";
2179             return 0;
2180           }
2181         }
2182         polys.clear();
2183         psd.get(polys);
2184         if(polys.size() == 0) {
2185           stdcout << "fail merge 10\n";
2186           return false;
2187         }
2188         stdcout << (polys[0]) << "\n";
2189       }
2190       for(unsigned int i = 0; i < 10; ++i) {
2191         stdcout << "random case # " << i << "\n";
2192         si.clear();
2193         pts.clear();
2194         pts.push_back(Point(rand()%9-4, rand()%9-4));
2195         pts.push_back(Point(rand()%9-4, rand()%9-4));
2196         pts.push_back(Point(rand()%9-4, rand()%9-4));
2197         polygon_data<Unit> poly1;
2198         poly1.set(pts.begin(), pts.end());
2199         stdcout << poly1 << "\n";
2200         si.insert(poly1, 444);
2201         pts.clear();
2202         pts.push_back(Point(rand()%9-4, rand()%9-4));
2203         pts.push_back(Point(rand()%9-4, rand()%9-4));
2204         pts.push_back(Point(rand()%9-4, rand()%9-4));
2205         polygon_data<Unit> poly2;
2206         poly2.set(pts.begin(), pts.end());
2207         stdcout << poly2 << "\n";
2208         si.insert(poly2, 444);
2209         result.clear();
2210         si.merge(result);
2211         print(stdcout, si.pmd) << "\n";
2212         if(!result.empty()) {
2213           psd = (*(result.begin())).second;
2214           stdcout << psd << "\n";
2215           polys.clear();
2216           psd.get(polys);
2217           if(polys.size() == 0) {
2218             si.clear();
2219             si.insert(poly1, 333);
2220             result.clear();
2221             si.merge(result);
2222             psd = (*(result.begin())).second;
2223             std::vector<polygon_data<Unit> > polys1;
2224             psd.get(polys1);
2225             si.clear();
2226             si.insert(poly2, 333);
2227             result.clear();
2228             si.merge(result);
2229             psd = (*(result.begin())).second;
2230             std::vector<polygon_data<Unit> > polys2;
2231             psd.get(polys2);
2232             if(!polys1.empty() || !polys2.empty()) {
2233               stdcout << "fail random merge " << i << "\n";
2234               return false;
2235             }
2236           }
2237         }
2238         if(!polys.empty())
2239           stdcout << polys.size() << ": " << (polys[0]) << "\n";
2240       }
2241       return true;
2242     }
2243 
2244     template <typename stream_type>
check_rectangle_trio(rectangle_data<Unit> rect1,rectangle_data<Unit> rect2,rectangle_data<Unit> rect3,stream_type & stdcout)2245     static inline bool check_rectangle_trio(rectangle_data<Unit> rect1, rectangle_data<Unit> rect2, rectangle_data<Unit> rect3, stream_type& stdcout) {
2246         property_merge si;
2247         std::map<std::set<property_type>, polygon_set_data<Unit> > result;
2248         std::vector<polygon_data<Unit> > polys;
2249         property_merge_90<property_type, Unit> si90;
2250         std::map<std::set<property_type>, polygon_90_set_data<Unit> > result90;
2251         std::vector<polygon_data<Unit> > polys90;
2252         si.insert(rect1, 111);
2253         si90.insert(rect1, 111);
2254         stdcout << rect1 << "\n";
2255         si.insert(rect2, 222);
2256         si90.insert(rect2, 222);
2257         stdcout << rect2 << "\n";
2258         si.insert(rect3, 333);
2259         si90.insert(rect3, 333);
2260         stdcout << rect3 << "\n";
2261         si.merge(result);
2262         si90.merge(result90);
2263         if(result.size() != result90.size()) {
2264           stdcout << "merge failed with size mismatch\n";
2265           return 0;
2266         }
2267         typename std::map<std::set<property_type>, polygon_90_set_data<Unit> >::iterator itr90 = result90.begin();
2268         for(typename std::map<std::set<property_type>, polygon_set_data<Unit> >::iterator itr = result.begin();
2269             itr != result.end(); ++itr) {
2270           for(typename std::set<property_type>::const_iterator set_itr = (*itr).first.begin();
2271               set_itr != (*itr).first.end(); ++set_itr) {
2272             stdcout << (*set_itr) << " ";
2273           } stdcout << ") \n";
2274           polygon_set_data<Unit> psd = (*itr).second;
2275           polygon_90_set_data<Unit> psd90 = (*itr90).second;
2276           polys.clear();
2277           polys90.clear();
2278           psd.get(polys);
2279           psd90.get(polys90);
2280           if(polys.size() != polys90.size()) {
2281             stdcout << "merge failed with polygon count mismatch\n";
2282             stdcout << psd << "\n";
2283             for(std::size_t j = 0; j < polys.size(); ++j) {
2284               stdcout << polys[j] << "\n";
2285             }
2286             stdcout << "reference\n";
2287             for(std::size_t j = 0; j < polys90.size(); ++j) {
2288               stdcout << polys90[j] << "\n";
2289             }
2290             return 0;
2291           }
2292           bool failed = false;
2293           for(std::size_t j = 0; j < polys.size(); ++j) {
2294             stdcout << polys[j] << "\n";
2295             stdcout << polys90[j] << "\n";
2296 #ifdef BOOST_POLYGON_ICC
2297 #pragma warning (push)
2298 #pragma warning (disable:1572)
2299 #endif
2300             if(area(polys[j]) != area(polys90[j])) {
2301 #ifdef BOOST_POLYGON_ICC
2302 #pragma warning (pop)
2303 #endif
2304               stdcout << "merge failed with area mismatch\n";
2305               failed = true;
2306             }
2307           }
2308           if(failed) return 0;
2309           ++itr90;
2310         }
2311         return true;
2312     }
2313 
2314     template <typename stream_type>
test_manhattan_intersection(stream_type & stdcout)2315     static inline bool test_manhattan_intersection(stream_type& stdcout) {
2316       rectangle_data<Unit> rect1, rect2, rect3;
2317       set_points(rect1, (Point(-1, 2)), (Point(1, 4)));
2318       set_points(rect2, (Point(-1, 2)), (Point(2, 3)));
2319       set_points(rect3, (Point(-3, 0)), (Point(4, 2)));
2320       if(!check_rectangle_trio(rect1, rect2, rect3, stdcout)) {
2321         return false;
2322       }
2323       for(unsigned int i = 0; i < 100; ++i) {
2324         property_merge si;
2325         std::map<std::set<property_type>, polygon_set_data<Unit> > result;
2326         std::vector<polygon_data<Unit> > polys;
2327         property_merge_90<property_type, Unit> si90;
2328         std::map<std::set<property_type>, polygon_90_set_data<Unit> > result90;
2329         std::vector<polygon_data<Unit> > polys90;
2330         stdcout << "random case # " << i << "\n";
2331         set_points(rect1, (Point(rand()%9-4, rand()%9-4)), (Point(rand()%9-4, rand()%9-4)));
2332         set_points(rect2, (Point(rand()%9-4, rand()%9-4)), (Point(rand()%9-4, rand()%9-4)));
2333         set_points(rect3, (Point(rand()%9-4, rand()%9-4)), (Point(rand()%9-4, rand()%9-4)));
2334         if(!check_rectangle_trio(rect1, rect2, rect3, stdcout)) {
2335           return false;
2336         }
2337       }
2338       return true;
2339     }
2340 
2341     template <typename stream_type>
test_intersection(stream_type & stdcout)2342     static inline bool test_intersection(stream_type& stdcout) {
2343       property_merge si;
2344       rectangle_data<Unit> rect;
2345       xl(rect, 0);
2346       yl(rect, 10);
2347       xh(rect, 30);
2348       yh(rect, 20);
2349       si.insert(rect, 333);
2350       xl(rect, 10);
2351       yl(rect, 0);
2352       xh(rect, 20);
2353       yh(rect, 30);
2354       si.insert(rect, 444);
2355       xl(rect, 15);
2356       yl(rect, 0);
2357       xh(rect, 25);
2358       yh(rect, 30);
2359       si.insert(rect, 555);
2360       std::map<std::set<property_type>, polygon_set_data<Unit> > result;
2361       si.merge(result);
2362       print(stdcout, si.pmd) << "\n";
2363       for(typename std::map<std::set<property_type>, polygon_set_data<Unit> >::iterator itr = result.begin();
2364           itr != result.end(); ++itr) {
2365         stdcout << "( ";
2366         for(typename std::set<property_type>::const_iterator set_itr = (*itr).first.begin();
2367             set_itr != (*itr).first.end(); ++set_itr) {
2368           stdcout << (*set_itr) << " ";
2369         } stdcout << ") \n";
2370         polygon_set_data<Unit> psd = (*itr).second;
2371         stdcout << psd << "\n";
2372         std::vector<polygon_data<Unit> > polys;
2373         psd.get(polys);
2374         for(std::size_t i = 0; i < polys.size(); ++i) {
2375           stdcout << polys[i] << "\n";
2376         }
2377       }
2378       std::vector<Point> pts;
2379       std::vector<polygon_data<Unit> > polys;
2380       for(unsigned int i = 0; i < 10; ++i) {
2381         property_merge si2;
2382         stdcout << "random case # " << i << "\n";
2383         si.clear();
2384         pts.clear();
2385         pts.push_back(Point(rand()%9-4, rand()%9-4));
2386         pts.push_back(Point(rand()%9-4, rand()%9-4));
2387         pts.push_back(Point(rand()%9-4, rand()%9-4));
2388         polygon_data<Unit> poly1;
2389         poly1.set(pts.begin(), pts.end());
2390         stdcout << poly1 << "\n";
2391         si.insert(poly1, 444);
2392         si2.insert(poly1, 333);
2393         pts.clear();
2394         pts.push_back(Point(rand()%9-4, rand()%9-4));
2395         pts.push_back(Point(rand()%9-4, rand()%9-4));
2396         pts.push_back(Point(rand()%9-4, rand()%9-4));
2397         polygon_data<Unit> poly2;
2398         poly2.set(pts.begin(), pts.end());
2399         stdcout << poly2 << "\n";
2400         si.insert(poly2, 444);
2401         si2.insert(poly2, 444);
2402         pts.clear();
2403         pts.push_back(Point(rand()%9-4, rand()%9-4));
2404         pts.push_back(Point(rand()%9-4, rand()%9-4));
2405         pts.push_back(Point(rand()%9-4, rand()%9-4));
2406         polygon_data<Unit> poly3;
2407         poly3.set(pts.begin(), pts.end());
2408         stdcout << poly3 << "\n";
2409         si.insert(poly3, 444);
2410         si2.insert(poly3, 555);
2411         result.clear();
2412         std::map<std::set<property_type>, polygon_set_data<Unit> > result2;
2413         si.merge(result);
2414         si2.merge(result2);
2415         stdcout << "merged result\n";
2416       for(typename std::map<std::set<property_type>, polygon_set_data<Unit> >::iterator itr = result.begin();
2417           itr != result.end(); ++itr) {
2418         stdcout << "( ";
2419         for(typename std::set<property_type>::const_iterator set_itr = (*itr).first.begin();
2420             set_itr != (*itr).first.end(); ++set_itr) {
2421           stdcout << (*set_itr) << " ";
2422         } stdcout << ") \n";
2423         polygon_set_data<Unit> psd = (*itr).second;
2424         stdcout << psd << "\n";
2425         std::vector<polygon_data<Unit> > polys2;
2426         psd.get(polys2);
2427         for(std::size_t ii = 0; ii < polys2.size(); ++ii) {
2428           stdcout << polys2[ii] << "\n";
2429         }
2430       }
2431       stdcout << "intersected pmd\n";
2432       print(stdcout, si2.pmd) << "\n";
2433       stdcout << "intersected result\n";
2434       for(typename std::map<std::set<property_type>, polygon_set_data<Unit> >::iterator itr = result2.begin();
2435           itr != result2.end(); ++itr) {
2436         stdcout << "( ";
2437         for(typename std::set<property_type>::const_iterator set_itr = (*itr).first.begin();
2438             set_itr != (*itr).first.end(); ++set_itr) {
2439           stdcout << (*set_itr) << " ";
2440         } stdcout << ") \n";
2441         polygon_set_data<Unit> psd = (*itr).second;
2442         stdcout << psd << "\n";
2443         std::vector<polygon_data<Unit> > polys2;
2444         psd.get(polys2);
2445         for(std::size_t ii = 0; ii < polys2.size(); ++ii) {
2446           stdcout << polys2[ii] << "\n";
2447         }
2448       }
2449         si.clear();
2450         for(typename std::map<std::set<property_type>, polygon_set_data<Unit> >::iterator itr = result2.begin();
2451             itr != result2.end(); ++itr) {
2452           polys.clear();
2453           (*itr).second.get(polys);
2454           for(std::size_t j = 0; j < polys.size(); ++j) {
2455             si.insert(polys[j], 444);
2456           }
2457         }
2458         result2.clear();
2459         si.merge(result2);
2460       stdcout << "remerged result\n";
2461       for(typename std::map<std::set<property_type>, polygon_set_data<Unit> >::iterator itr = result2.begin();
2462           itr != result2.end(); ++itr) {
2463         stdcout << "( ";
2464         for(typename std::set<property_type>::const_iterator set_itr = (*itr).first.begin();
2465             set_itr != (*itr).first.end(); ++set_itr) {
2466           stdcout << (*set_itr) << " ";
2467         } stdcout << ") \n";
2468         polygon_set_data<Unit> psd = (*itr).second;
2469         stdcout << psd << "\n";
2470         std::vector<polygon_data<Unit> > polys2;
2471         psd.get(polys2);
2472         for(std::size_t ii = 0; ii < polys2.size(); ++ii) {
2473           stdcout << polys2[ii] << "\n";
2474         }
2475       }
2476       std::vector<polygon_data<Unit> > polys2;
2477       polys.clear();
2478       (*(result.begin())).second.get(polys);
2479       (*(result2.begin())).second.get(polys2);
2480       if(!(polys == polys2)) {
2481           stdcout << "failed intersection check # " << i << "\n";
2482           return false;
2483         }
2484       }
2485       return true;
2486     }
2487   };
2488 
2489   template <typename Unit>
2490   class arbitrary_boolean_op : public scanline_base<Unit> {
2491   private:
2492 
2493     typedef int property_type;
2494     typedef typename scanline_base<Unit>::Point Point;
2495 
2496     //the first point is the vertex and and second point establishes the slope of an edge eminating from the vertex
2497     //typedef std::pair<Point, Point> half_edge;
2498     typedef typename scanline_base<Unit>::half_edge half_edge;
2499 
2500     //scanline comparator functor
2501     typedef typename scanline_base<Unit>::less_half_edge less_half_edge;
2502     typedef typename scanline_base<Unit>::less_point less_point;
2503 
2504     //this data structure assocates a property and count to a half edge
2505     typedef std::pair<half_edge, std::pair<property_type, int> > vertex_property;
2506     //this data type stores the combination of many half edges
2507     typedef std::vector<vertex_property> property_merge_data;
2508 
2509     //this is the data type used internally to store the combination of property counts at a given location
2510     typedef std::vector<std::pair<property_type, int> > property_map;
2511     //this data type is used internally to store the combined property data for a given half edge
2512     typedef std::pair<half_edge, property_map> vertex_data;
2513 
2514     property_merge_data pmd;
2515     typename scanline_base<Unit>::evalAtXforYPack evalAtXforYPack_;
2516 
2517     template<typename vertex_data_type>
2518     class less_vertex_data {
2519       typename scanline_base<Unit>::evalAtXforYPack* pack_;
2520     public:
less_vertex_data()2521       less_vertex_data() : pack_() {}
less_vertex_data(typename scanline_base<Unit>::evalAtXforYPack * pack)2522       less_vertex_data(typename scanline_base<Unit>::evalAtXforYPack* pack) : pack_(pack) {}
operator ()(const vertex_data_type & lvalue,const vertex_data_type & rvalue) const2523       bool operator()(const vertex_data_type& lvalue, const vertex_data_type& rvalue) const {
2524         less_point lp;
2525         if(lp(lvalue.first.first, rvalue.first.first)) return true;
2526         if(lp(rvalue.first.first, lvalue.first.first)) return false;
2527         Unit x = lvalue.first.first.get(HORIZONTAL);
2528         int just_before_ = 0;
2529         less_half_edge lhe(&x, &just_before_, pack_);
2530         return lhe(lvalue.first, rvalue.first);
2531       }
2532     };
2533 
2534     template <typename result_type, typename key_type, int op_type>
2535     class boolean_output_functor {
2536     public:
boolean_output_functor()2537       boolean_output_functor() {}
operator ()(result_type & result,const half_edge & edge,const key_type & left,const key_type & right)2538       void operator()(result_type& result, const half_edge& edge, const key_type& left, const key_type& right) {
2539         typename std::pair<half_edge, int> elem;
2540         elem.first = edge;
2541         elem.second = 1;
2542         if(edge.second < edge.first) elem.second *= -1;
2543         if(scanline_base<Unit>::is_vertical(edge)) elem.second *= -1;
2544 #ifdef BOOST_POLYGON_MSVC
2545 #pragma warning (push)
2546 #pragma warning (disable: 4127)
2547 #endif
2548         if(op_type == 0) { //OR
2549           if(!left.empty() && right.empty()) {
2550             result.insert_clean(elem);
2551           } else if(!right.empty() && left.empty()) {
2552             elem.second *= -1;
2553             result.insert_clean(elem);
2554           }
2555         } else if(op_type == 1) { //AND
2556           if(left.size() == 2 && right.size() != 2) {
2557             result.insert_clean(elem);
2558           } else if(right.size() == 2 && left.size() != 2) {
2559             elem.second *= -1;
2560             result.insert_clean(elem);
2561           }
2562         } else if(op_type == 2) { //XOR
2563           if(left.size() == 1 && right.size() != 1) {
2564             result.insert_clean(elem);
2565           } else if(right.size() == 1 && left.size() != 1) {
2566             elem.second *= -1;
2567             result.insert_clean(elem);
2568           }
2569         } else { //SUBTRACT
2570           if(left.size() == 1) {
2571             if((*(left.begin())) == 0) {
2572               result.insert_clean(elem);
2573             }
2574           }
2575 #ifdef BOOST_POLYGON_MSVC
2576 #pragma warning (pop)
2577 #endif
2578           if(right.size() == 1) {
2579             if((*(right.begin())) == 0) {
2580               elem.second *= -1;
2581               result.insert_clean(elem);
2582             }
2583           }
2584         }
2585       }
2586     };
2587 
sort_property_merge_data()2588     inline void sort_property_merge_data() {
2589       less_vertex_data<vertex_property> lvd(&evalAtXforYPack_);
2590       polygon_sort(pmd.begin(), pmd.end(), lvd);
2591     }
2592   public:
arbitrary_boolean_op()2593     inline arbitrary_boolean_op() : pmd(), evalAtXforYPack_() {}
arbitrary_boolean_op(const arbitrary_boolean_op & pm)2594     inline arbitrary_boolean_op(const arbitrary_boolean_op& pm) : pmd(pm.pmd), evalAtXforYPack_(pm.evalAtXforYPack_) {}
operator =(const arbitrary_boolean_op & pm)2595     inline arbitrary_boolean_op& operator=(const arbitrary_boolean_op& pm) { pmd = pm.pmd; return *this; }
2596 
2597     enum BOOLEAN_OP_TYPE {
2598       BOOLEAN_OR = 0,
2599       BOOLEAN_AND = 1,
2600       BOOLEAN_XOR = 2,
2601       BOOLEAN_NOT = 3
2602     };
2603     template <typename result_type, typename iT1, typename iT2>
execute(result_type & result,iT1 b1,iT1 e1,iT2 b2,iT2 e2,int op)2604     inline void execute(result_type& result, iT1 b1, iT1 e1, iT2 b2, iT2 e2, int op) {
2605       //intersect data
2606       insert(b1, e1, 0);
2607       insert(b2, e2, 1);
2608       property_merge_data tmp_pmd;
2609       //#define BOOST_POLYGON_DEBUG_FILE
2610 #ifdef BOOST_POLYGON_DEBUG_FILE
2611       std::fstream debug_file;
2612       debug_file.open("gtl_debug.txt", std::ios::out);
2613       property_merge<Unit, property_type, std::vector<property_type> >::print(debug_file, pmd);
2614       debug_file.close();
2615 #endif
2616       if(pmd.empty())
2617         return;
2618       line_intersection<Unit>::validate_scan(tmp_pmd, pmd.begin(), pmd.end());
2619       pmd.swap(tmp_pmd);
2620       sort_property_merge_data();
2621       scanline<Unit, property_type, std::vector<property_type> > sl;
2622       if(op == BOOLEAN_OR) {
2623         boolean_output_functor<result_type, std::vector<property_type>, 0> bof;
2624         sl.scan(result, bof, pmd.begin(), pmd.end());
2625       } else if(op == BOOLEAN_AND) {
2626         boolean_output_functor<result_type, std::vector<property_type>, 1> bof;
2627         sl.scan(result, bof, pmd.begin(), pmd.end());
2628       } else if(op == BOOLEAN_XOR) {
2629         boolean_output_functor<result_type, std::vector<property_type>, 2> bof;
2630         sl.scan(result, bof, pmd.begin(), pmd.end());
2631       } else if(op == BOOLEAN_NOT) {
2632         boolean_output_functor<result_type, std::vector<property_type>, 3> bof;
2633         sl.scan(result, bof, pmd.begin(), pmd.end());
2634       }
2635     }
2636 
clear()2637     inline void clear() {*this = arbitrary_boolean_op();}
2638 
2639   private:
2640     template <typename iT>
insert(iT b,iT e,int id)2641     void insert(iT b, iT e, int id) {
2642       for(;
2643           b != e; ++b) {
2644         pmd.push_back(vertex_property(half_edge((*b).first.first, (*b).first.second),
2645                                       std::pair<property_type, int>(id, (*b).second)));
2646       }
2647     }
2648 
2649   };
2650 
2651   template <typename Unit, typename stream_type>
test_arbitrary_boolean_op(stream_type & stdcout)2652   bool test_arbitrary_boolean_op(stream_type& stdcout) {
2653     polygon_set_data<Unit> psd;
2654     rectangle_data<Unit> rect;
2655     set_points(rect, point_data<Unit>(0, 0), point_data<Unit>(10, 10));
2656     psd.insert(rect);
2657     polygon_set_data<Unit> psd2;
2658     set_points(rect, point_data<Unit>(5, 5), point_data<Unit>(15, 15));
2659     psd2.insert(rect);
2660     std::vector<polygon_data<Unit> > pv;
2661     pv.clear();
2662     arbitrary_boolean_op<Unit> abo;
2663     polygon_set_data<Unit> psd3;
2664     abo.execute(psd3, psd.begin(), psd.end(), psd2.begin(), psd2.end(), arbitrary_boolean_op<Unit>::BOOLEAN_OR);
2665     psd3.get(pv);
2666     for(std::size_t i = 0; i < pv.size(); ++i) {
2667       stdcout << pv[i] << "\n";
2668     }
2669     pv.clear();
2670     abo.clear();
2671     psd3.clear();
2672     abo.execute(psd3, psd.begin(), psd.end(), psd2.begin(), psd2.end(), arbitrary_boolean_op<Unit>::BOOLEAN_AND);
2673     psd3.get(pv);
2674     for(std::size_t i = 0; i < pv.size(); ++i) {
2675       stdcout << pv[i] << "\n";
2676     }
2677     pv.clear();
2678     abo.clear();
2679     psd3.clear();
2680     abo.execute(psd3, psd.begin(), psd.end(), psd2.begin(), psd2.end(), arbitrary_boolean_op<Unit>::BOOLEAN_XOR);
2681     psd3.get(pv);
2682     for(std::size_t i = 0; i < pv.size(); ++i) {
2683       stdcout << pv[i] << "\n";
2684     }
2685     pv.clear();
2686     abo.clear();
2687     psd3.clear();
2688     abo.execute(psd3, psd.begin(), psd.end(), psd2.begin(), psd2.end(), arbitrary_boolean_op<Unit>::BOOLEAN_NOT);
2689     psd3.get(pv);
2690     for(std::size_t i = 0; i < pv.size(); ++i) {
2691       stdcout << pv[i] << "\n";
2692     }
2693     return true;
2694   }
2695 
2696 
2697 
2698 
2699 
2700 
2701 
2702 
2703 
2704 
2705 
2706 
2707 
2708 
2709   template <typename Unit, typename property_type>
2710   class arbitrary_connectivity_extraction : public scanline_base<Unit> {
2711   private:
2712 
2713     typedef typename scanline_base<Unit>::Point Point;
2714 
2715     //the first point is the vertex and and second point establishes the slope of an edge eminating from the vertex
2716     //typedef std::pair<Point, Point> half_edge;
2717     typedef typename scanline_base<Unit>::half_edge half_edge;
2718 
2719     //scanline comparator functor
2720     typedef typename scanline_base<Unit>::less_half_edge less_half_edge;
2721     typedef typename scanline_base<Unit>::less_point less_point;
2722 
2723     //this data structure assocates a property and count to a half edge
2724     typedef std::pair<half_edge, std::pair<property_type, int> > vertex_property;
2725     //this data type stores the combination of many half edges
2726     typedef std::vector<vertex_property> property_merge_data;
2727 
2728     //this is the data type used internally to store the combination of property counts at a given location
2729     typedef std::vector<std::pair<property_type, int> > property_map;
2730     //this data type is used internally to store the combined property data for a given half edge
2731     typedef std::pair<half_edge, property_map> vertex_data;
2732 
2733     property_merge_data pmd;
2734     typename scanline_base<Unit>::evalAtXforYPack evalAtXforYPack_;
2735 
2736     template<typename vertex_data_type>
2737     class less_vertex_data {
2738       typename scanline_base<Unit>::evalAtXforYPack* pack_;
2739     public:
less_vertex_data()2740       less_vertex_data() : pack_() {}
less_vertex_data(typename scanline_base<Unit>::evalAtXforYPack * pack)2741       less_vertex_data(typename scanline_base<Unit>::evalAtXforYPack* pack) : pack_(pack) {}
operator ()(const vertex_data_type & lvalue,const vertex_data_type & rvalue) const2742       bool operator()(const vertex_data_type& lvalue, const vertex_data_type& rvalue) const {
2743         less_point lp;
2744         if(lp(lvalue.first.first, rvalue.first.first)) return true;
2745         if(lp(rvalue.first.first, lvalue.first.first)) return false;
2746         Unit x = lvalue.first.first.get(HORIZONTAL);
2747         int just_before_ = 0;
2748         less_half_edge lhe(&x, &just_before_, pack_);
2749         return lhe(lvalue.first, rvalue.first);
2750       }
2751     };
2752 
2753 
2754     template <typename cT>
process_previous_x(cT & output)2755     static void process_previous_x(cT& output) {
2756       std::map<point_data<Unit>, std::set<property_type> >& y_prop_map = output.first.second;
2757       if(y_prop_map.empty()) return;
2758       Unit x = output.first.first;
2759       for(typename std::map<point_data<Unit>, std::set<property_type> >::iterator itr =
2760             y_prop_map.begin(); itr != y_prop_map.end(); ++itr) {
2761         if((*itr).first.x() < x) {
2762           y_prop_map.erase(y_prop_map.begin(), itr);
2763           continue;
2764         }
2765         for(typename std::set<property_type>::iterator inner_itr = itr->second.begin();
2766             inner_itr != itr->second.end(); ++inner_itr) {
2767           std::set<property_type>& output_edges = (*(output.second))[*inner_itr];
2768           typename std::set<property_type>::iterator inner_inner_itr = inner_itr;
2769           ++inner_inner_itr;
2770           for( ; inner_inner_itr != itr->second.end(); ++inner_inner_itr) {
2771             output_edges.insert(output_edges.end(), *inner_inner_itr);
2772             std::set<property_type>& output_edges_2 = (*(output.second))[*inner_inner_itr];
2773             output_edges_2.insert(output_edges_2.end(), *inner_itr);
2774           }
2775         }
2776       }
2777     }
2778 
2779     template <typename result_type, typename key_type>
2780     class connectivity_extraction_output_functor {
2781     public:
connectivity_extraction_output_functor()2782       connectivity_extraction_output_functor() {}
operator ()(result_type & result,const half_edge & edge,const key_type & left,const key_type & right)2783       void operator()(result_type& result, const half_edge& edge, const key_type& left, const key_type& right) {
2784         Unit& x = result.first.first;
2785         std::map<point_data<Unit>, std::set<property_type> >& y_prop_map = result.first.second;
2786         point_data<Unit> pt = edge.first;
2787         if(pt.x() != x) process_previous_x(result);
2788         x = pt.x();
2789         std::set<property_type>& output_set = y_prop_map[pt];
2790         {
2791           for(typename key_type::const_iterator itr1 =
2792                 left.begin(); itr1 != left.end(); ++itr1) {
2793             output_set.insert(output_set.end(), *itr1);
2794           }
2795           for(typename key_type::const_iterator itr2 =
2796                 right.begin(); itr2 != right.end(); ++itr2) {
2797             output_set.insert(output_set.end(), *itr2);
2798           }
2799         }
2800         std::set<property_type>& output_set2 = y_prop_map[edge.second];
2801         for(typename key_type::const_iterator itr1 =
2802               left.begin(); itr1 != left.end(); ++itr1) {
2803           output_set2.insert(output_set2.end(), *itr1);
2804         }
2805         for(typename key_type::const_iterator itr2 =
2806               right.begin(); itr2 != right.end(); ++itr2) {
2807           output_set2.insert(output_set2.end(), *itr2);
2808         }
2809       }
2810     };
2811 
sort_property_merge_data()2812     inline void sort_property_merge_data() {
2813       less_vertex_data<vertex_property> lvd(&evalAtXforYPack_);
2814       polygon_sort(pmd.begin(), pmd.end(), lvd);
2815     }
2816   public:
arbitrary_connectivity_extraction()2817     inline arbitrary_connectivity_extraction() : pmd(), evalAtXforYPack_() {}
arbitrary_connectivity_extraction(const arbitrary_connectivity_extraction & pm)2818     inline arbitrary_connectivity_extraction
2819     (const arbitrary_connectivity_extraction& pm) : pmd(pm.pmd), evalAtXforYPack_(pm.evalAtXforYPack_) {}
operator =(const arbitrary_connectivity_extraction & pm)2820     inline arbitrary_connectivity_extraction& operator=
2821       (const arbitrary_connectivity_extraction& pm) { pmd = pm.pmd; return *this; }
2822 
2823     template <typename result_type>
execute(result_type & result)2824     inline void execute(result_type& result) {
2825       //intersect data
2826       property_merge_data tmp_pmd;
2827       line_intersection<Unit>::validate_scan(tmp_pmd, pmd.begin(), pmd.end());
2828       pmd.swap(tmp_pmd);
2829       sort_property_merge_data();
2830       scanline<Unit, property_type, std::vector<property_type> > sl;
2831       std::pair<std::pair<Unit, std::map<point_data<Unit>, std::set<property_type> > >,
2832         result_type*> output
2833         (std::make_pair(std::make_pair((std::numeric_limits<Unit>::max)(),
2834                                        std::map<point_data<Unit>,
2835                                        std::set<property_type> >()), &result));
2836       connectivity_extraction_output_functor<std::pair<std::pair<Unit,
2837         std::map<point_data<Unit>, std::set<property_type> > >, result_type*>,
2838         std::vector<property_type> > ceof;
2839       sl.scan(output, ceof, pmd.begin(), pmd.end());
2840       process_previous_x(output);
2841     }
2842 
clear()2843     inline void clear() {*this = arbitrary_connectivity_extraction();}
2844 
2845     template <typename iT>
populateTouchSetData(iT begin,iT end,property_type property)2846     void populateTouchSetData(iT begin, iT end,
2847                                      property_type property) {
2848       for( ; begin != end; ++begin) {
2849         pmd.push_back(vertex_property(half_edge((*begin).first.first, (*begin).first.second),
2850                                       std::pair<property_type, int>(property, (*begin).second)));
2851       }
2852     }
2853 
2854   };
2855 
2856 }
2857 }
2858 #endif
2859