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_POLYGON_ARBITRARY_FORMATION_HPP
9 #define BOOST_POLYGON_POLYGON_ARBITRARY_FORMATION_HPP
10 namespace boost { namespace polygon{
11   template <typename T, typename T2>
12   struct PolyLineArbitraryByConcept {};
13 
14   template <typename T>
15   class poly_line_arbitrary_polygon_data;
16   template <typename T>
17   class poly_line_arbitrary_hole_data;
18 
19   template <typename Unit>
20   struct scanline_base {
21 
22     typedef point_data<Unit> Point;
23     typedef std::pair<Point, Point> half_edge;
24 
25     class less_point {
26     public:
27       typedef Point first_argument_type;
28       typedef Point second_argument_type;
29       typedef bool result_type;
less_point()30       inline less_point() {}
operator ()(const Point & pt1,const Point & pt2) const31       inline bool operator () (const Point& pt1, const Point& pt2) const {
32         if(pt1.get(HORIZONTAL) < pt2.get(HORIZONTAL)) return true;
33         if(pt1.get(HORIZONTAL) == pt2.get(HORIZONTAL)) {
34           if(pt1.get(VERTICAL) < pt2.get(VERTICAL)) return true;
35         }
36         return false;
37       }
38     };
39 
betweenboost::polygon::scanline_base40     static inline bool between(Point pt, Point pt1, Point pt2) {
41       less_point lp;
42       if(lp(pt1, pt2))
43         return lp(pt, pt2) && lp(pt1, pt);
44       return lp(pt, pt1) && lp(pt2, pt);
45     }
46 
47     template <typename area_type>
compute_interceptboost::polygon::scanline_base48     static inline Unit compute_intercept(const area_type& dy2,
49                                          const area_type& dx1,
50                                          const area_type& dx2) {
51       //intercept = dy2 * dx1 / dx2
52       //return (Unit)(((area_type)dy2 * (area_type)dx1) / (area_type)dx2);
53       area_type dx1_q = dx1 / dx2;
54       area_type dx1_r = dx1 % dx2;
55       return dx1_q * dy2 + (dy2 * dx1_r)/dx2;
56     }
57 
58     template <typename area_type>
equal_slopeboost::polygon::scanline_base59     static inline bool equal_slope(area_type dx1, area_type dy1, area_type dx2, area_type dy2) {
60       typedef typename coordinate_traits<Unit>::unsigned_area_type unsigned_product_type;
61       unsigned_product_type cross_1 = (unsigned_product_type)(dx2 < 0 ? -dx2 :dx2) * (unsigned_product_type)(dy1 < 0 ? -dy1 : dy1);
62       unsigned_product_type cross_2 = (unsigned_product_type)(dx1 < 0 ? -dx1 :dx1) * (unsigned_product_type)(dy2 < 0 ? -dy2 : dy2);
63       int dx1_sign = dx1 < 0 ? -1 : 1;
64       int dx2_sign = dx2 < 0 ? -1 : 1;
65       int dy1_sign = dy1 < 0 ? -1 : 1;
66       int dy2_sign = dy2 < 0 ? -1 : 1;
67       int cross_1_sign = dx2_sign * dy1_sign;
68       int cross_2_sign = dx1_sign * dy2_sign;
69       return cross_1 == cross_2 && (cross_1_sign == cross_2_sign || cross_1 == 0);
70     }
71 
72     template <typename T>
equal_slope_hpboost::polygon::scanline_base73     static inline bool equal_slope_hp(const T& dx1, const T& dy1, const T& dx2, const T& dy2) {
74       return dx1 * dy2 == dx2 * dy1;
75     }
76 
equal_slopeboost::polygon::scanline_base77     static inline bool equal_slope(const Unit& x, const Unit& y,
78                                    const Point& pt1, const Point& pt2) {
79       const Point* pts[2] = {&pt1, &pt2};
80       typedef typename coordinate_traits<Unit>::manhattan_area_type at;
81       at dy2 = (at)pts[1]->get(VERTICAL) - (at)y;
82       at dy1 = (at)pts[0]->get(VERTICAL) - (at)y;
83       at dx2 = (at)pts[1]->get(HORIZONTAL) - (at)x;
84       at dx1 = (at)pts[0]->get(HORIZONTAL) - (at)x;
85       return equal_slope(dx1, dy1, dx2, dy2);
86     }
87 
88     template <typename area_type>
less_slopeboost::polygon::scanline_base89     static inline bool less_slope(area_type dx1, area_type dy1, area_type dx2, area_type dy2) {
90       //reflext x and y slopes to right hand side half plane
91       if(dx1 < 0) {
92         dy1 *= -1;
93         dx1 *= -1;
94       } else if(dx1 == 0) {
95         //if the first slope is vertical the first cannot be less
96         return false;
97       }
98       if(dx2 < 0) {
99         dy2 *= -1;
100         dx2 *= -1;
101       } else if(dx2 == 0) {
102         //if the second slope is vertical the first is always less unless it is also vertical, in which case they are equal
103         return dx1 != 0;
104       }
105       typedef typename coordinate_traits<Unit>::unsigned_area_type unsigned_product_type;
106       unsigned_product_type cross_1 = (unsigned_product_type)(dx2 < 0 ? -dx2 :dx2) * (unsigned_product_type)(dy1 < 0 ? -dy1 : dy1);
107       unsigned_product_type cross_2 = (unsigned_product_type)(dx1 < 0 ? -dx1 :dx1) * (unsigned_product_type)(dy2 < 0 ? -dy2 : dy2);
108       int dx1_sign = dx1 < 0 ? -1 : 1;
109       int dx2_sign = dx2 < 0 ? -1 : 1;
110       int dy1_sign = dy1 < 0 ? -1 : 1;
111       int dy2_sign = dy2 < 0 ? -1 : 1;
112       int cross_1_sign = dx2_sign * dy1_sign;
113       int cross_2_sign = dx1_sign * dy2_sign;
114       if(cross_1_sign < cross_2_sign) return true;
115       if(cross_2_sign < cross_1_sign) return false;
116       if(cross_1_sign == -1) return cross_2 < cross_1;
117       return cross_1 < cross_2;
118     }
119 
less_slopeboost::polygon::scanline_base120     static inline bool less_slope(const Unit& x, const Unit& y,
121                                   const Point& pt1, const Point& pt2) {
122       const Point* pts[2] = {&pt1, &pt2};
123       //compute y value on edge from pt_ to pts[1] at the x value of pts[0]
124       typedef typename coordinate_traits<Unit>::manhattan_area_type at;
125       at dy2 = (at)pts[1]->get(VERTICAL) - (at)y;
126       at dy1 = (at)pts[0]->get(VERTICAL) - (at)y;
127       at dx2 = (at)pts[1]->get(HORIZONTAL) - (at)x;
128       at dx1 = (at)pts[0]->get(HORIZONTAL) - (at)x;
129       return less_slope(dx1, dy1, dx2, dy2);
130     }
131 
132     //return -1 below, 0 on and 1 above line
on_above_or_belowboost::polygon::scanline_base133     static inline int on_above_or_below(Point pt, const half_edge& he) {
134       if(pt == he.first || pt == he.second) return 0;
135       if(equal_slope(pt.get(HORIZONTAL), pt.get(VERTICAL), he.first, he.second)) return 0;
136       bool less_result = less_slope(pt.get(HORIZONTAL), pt.get(VERTICAL), he.first, he.second);
137       int retval = less_result ? -1 : 1;
138       less_point lp;
139       if(lp(he.second, he.first)) retval *= -1;
140       if(!between(pt, he.first, he.second)) retval *= -1;
141       return retval;
142     }
143 
144     //returns true is the segment intersects the integer grid square with lower
145     //left corner at point
intersects_gridboost::polygon::scanline_base146     static inline bool intersects_grid(Point pt, const half_edge& he) {
147       if(pt == he.second) return true;
148       if(pt == he.first) return true;
149       rectangle_data<Unit> rect1;
150       set_points(rect1, he.first, he.second);
151       if(contains(rect1, pt, true)) {
152         if(is_vertical(he) || is_horizontal(he)) return true;
153       } else {
154         return false; //can't intersect a grid not within bounding box
155       }
156       Unit x = pt.get(HORIZONTAL);
157       Unit y = pt.get(VERTICAL);
158       if(equal_slope(x, y, he.first, he.second) &&
159          between(pt, he.first, he.second)) return true;
160       Point pt01(pt.get(HORIZONTAL), pt.get(VERTICAL) + 1);
161       Point pt10(pt.get(HORIZONTAL) + 1, pt.get(VERTICAL));
162       Point pt11(pt.get(HORIZONTAL) + 1, pt.get(VERTICAL) + 1);
163 //       if(pt01 == he.first) return true;
164 //       if(pt10 == he.first) return true;
165 //       if(pt11 == he.first) return true;
166 //       if(pt01 == he.second) return true;
167 //       if(pt10 == he.second) return true;
168 //       if(pt11 == he.second) return true;
169       //check non-integer intersections
170       half_edge widget1(pt, pt11);
171       //intersects but not just at pt11
172       if(intersects(widget1, he) && on_above_or_below(pt11, he)) return true;
173       half_edge widget2(pt01, pt10);
174       //intersects but not just at pt01 or 10
175       if(intersects(widget2, he) && on_above_or_below(pt01, he) && on_above_or_below(pt10, he)) return true;
176       return false;
177     }
178 
evalAtXforYlazyboost::polygon::scanline_base179     static inline Unit evalAtXforYlazy(Unit xIn, Point pt, Point other_pt) {
180       long double
181         evalAtXforYret, evalAtXforYxIn, evalAtXforYx1, evalAtXforYy1, evalAtXforYdx1, evalAtXforYdx,
182         evalAtXforYdy, evalAtXforYx2, evalAtXforYy2, evalAtXforY0;
183       //y = (x - x1)dy/dx + y1
184       //y = (xIn - pt.x)*(other_pt.y-pt.y)/(other_pt.x-pt.x) + pt.y
185       //assert pt.x != other_pt.x
186       if(pt.y() == other_pt.y())
187         return pt.y();
188       evalAtXforYxIn = xIn;
189       evalAtXforYx1 = pt.get(HORIZONTAL);
190       evalAtXforYy1 = pt.get(VERTICAL);
191       evalAtXforYdx1 = evalAtXforYxIn - evalAtXforYx1;
192       evalAtXforY0 = 0;
193       if(evalAtXforYdx1 == evalAtXforY0) return (Unit)evalAtXforYy1;
194       evalAtXforYx2 = other_pt.get(HORIZONTAL);
195       evalAtXforYy2 = other_pt.get(VERTICAL);
196 
197       evalAtXforYdx = evalAtXforYx2 - evalAtXforYx1;
198       evalAtXforYdy = evalAtXforYy2 - evalAtXforYy1;
199       evalAtXforYret = ((evalAtXforYdx1) * evalAtXforYdy / evalAtXforYdx + evalAtXforYy1);
200       return (Unit)evalAtXforYret;
201     }
202 
evalAtXforYboost::polygon::scanline_base203     static inline typename high_precision_type<Unit>::type evalAtXforY(Unit xIn, Point pt, Point other_pt) {
204       typename high_precision_type<Unit>::type
205         evalAtXforYret, evalAtXforYxIn, evalAtXforYx1, evalAtXforYy1, evalAtXforYdx1, evalAtXforYdx,
206         evalAtXforYdy, evalAtXforYx2, evalAtXforYy2, evalAtXforY0;
207       //y = (x - x1)dy/dx + y1
208       //y = (xIn - pt.x)*(other_pt.y-pt.y)/(other_pt.x-pt.x) + pt.y
209       //assert pt.x != other_pt.x
210       typedef typename high_precision_type<Unit>::type high_precision;
211       if(pt.y() == other_pt.y())
212         return (high_precision)pt.y();
213       evalAtXforYxIn = (high_precision)xIn;
214       evalAtXforYx1 = pt.get(HORIZONTAL);
215       evalAtXforYy1 = pt.get(VERTICAL);
216       evalAtXforYdx1 = evalAtXforYxIn - evalAtXforYx1;
217       evalAtXforY0 = high_precision(0);
218       if(evalAtXforYdx1 == evalAtXforY0) return evalAtXforYret = evalAtXforYy1;
219       evalAtXforYx2 = (high_precision)other_pt.get(HORIZONTAL);
220       evalAtXforYy2 = (high_precision)other_pt.get(VERTICAL);
221 
222       evalAtXforYdx = evalAtXforYx2 - evalAtXforYx1;
223       evalAtXforYdy = evalAtXforYy2 - evalAtXforYy1;
224       evalAtXforYret = ((evalAtXforYdx1) * evalAtXforYdy / evalAtXforYdx + evalAtXforYy1);
225       return evalAtXforYret;
226     }
227 
228     struct evalAtXforYPack {
229     typename high_precision_type<Unit>::type
230     evalAtXforYret, evalAtXforYxIn, evalAtXforYx1, evalAtXforYy1, evalAtXforYdx1, evalAtXforYdx,
231                            evalAtXforYdy, evalAtXforYx2, evalAtXforYy2, evalAtXforY0;
evalAtXforYboost::polygon::scanline_base::evalAtXforYPack232       inline const typename high_precision_type<Unit>::type& evalAtXforY(Unit xIn, Point pt, Point other_pt) {
233         //y = (x - x1)dy/dx + y1
234         //y = (xIn - pt.x)*(other_pt.y-pt.y)/(other_pt.x-pt.x) + pt.y
235         //assert pt.x != other_pt.x
236         typedef typename high_precision_type<Unit>::type high_precision;
237         if(pt.y() == other_pt.y()) {
238           evalAtXforYret = (high_precision)pt.y();
239           return evalAtXforYret;
240         }
241         evalAtXforYxIn = (high_precision)xIn;
242         evalAtXforYx1 = pt.get(HORIZONTAL);
243         evalAtXforYy1 = pt.get(VERTICAL);
244         evalAtXforYdx1 = evalAtXforYxIn - evalAtXforYx1;
245         evalAtXforY0 = high_precision(0);
246         if(evalAtXforYdx1 == evalAtXforY0) return evalAtXforYret = evalAtXforYy1;
247         evalAtXforYx2 = (high_precision)other_pt.get(HORIZONTAL);
248         evalAtXforYy2 = (high_precision)other_pt.get(VERTICAL);
249 
250         evalAtXforYdx = evalAtXforYx2 - evalAtXforYx1;
251         evalAtXforYdy = evalAtXforYy2 - evalAtXforYy1;
252         evalAtXforYret = ((evalAtXforYdx1) * evalAtXforYdy / evalAtXforYdx + evalAtXforYy1);
253         return evalAtXforYret;
254       }
255     };
256 
is_verticalboost::polygon::scanline_base257     static inline bool is_vertical(const half_edge& he) {
258       return he.first.get(HORIZONTAL) == he.second.get(HORIZONTAL);
259     }
260 
is_horizontalboost::polygon::scanline_base261     static inline bool is_horizontal(const half_edge& he) {
262       return he.first.get(VERTICAL) == he.second.get(VERTICAL);
263     }
264 
is_45_degreeboost::polygon::scanline_base265     static inline bool is_45_degree(const half_edge& he) {
266       return euclidean_distance(he.first, he.second, HORIZONTAL) == euclidean_distance(he.first, he.second, VERTICAL);
267     }
268 
269     //scanline comparator functor
270     class less_half_edge {
271     private:
272       Unit *x_; //x value at which to apply comparison
273       int *justBefore_;
274       evalAtXforYPack * pack_;
275     public:
276       typedef half_edge first_argument_type;
277       typedef half_edge second_argument_type;
278       typedef bool result_type;
less_half_edge()279       inline less_half_edge() : x_(0), justBefore_(0), pack_(0) {}
less_half_edge(Unit * x,int * justBefore,evalAtXforYPack * packIn)280       inline less_half_edge(Unit *x, int *justBefore, evalAtXforYPack * packIn) : x_(x), justBefore_(justBefore), pack_(packIn) {}
less_half_edge(const less_half_edge & that)281       inline less_half_edge(const less_half_edge& that) : x_(that.x_), justBefore_(that.justBefore_),
282                                                           pack_(that.pack_){}
operator =(const less_half_edge & that)283       inline less_half_edge& operator=(const less_half_edge& that) {
284         x_ = that.x_;
285         justBefore_ = that.justBefore_;
286         pack_ = that.pack_;
287         return *this; }
operator ()(const half_edge & elm1,const half_edge & elm2) const288       inline bool operator () (const half_edge& elm1, const half_edge& elm2) const {
289         if((std::max)(elm1.first.y(), elm1.second.y()) < (std::min)(elm2.first.y(), elm2.second.y()))
290           return true;
291         if((std::min)(elm1.first.y(), elm1.second.y()) > (std::max)(elm2.first.y(), elm2.second.y()))
292           return false;
293 
294         //check if either x of elem1 is equal to x_
295         Unit localx = *x_;
296         Unit elm1y = 0;
297         bool elm1_at_x = false;
298         if(localx == elm1.first.get(HORIZONTAL)) {
299           elm1_at_x = true;
300           elm1y = elm1.first.get(VERTICAL);
301         } else if(localx == elm1.second.get(HORIZONTAL)) {
302           elm1_at_x = true;
303           elm1y = elm1.second.get(VERTICAL);
304         }
305         Unit elm2y = 0;
306         bool elm2_at_x = false;
307         if(localx == elm2.first.get(HORIZONTAL)) {
308           elm2_at_x = true;
309           elm2y = elm2.first.get(VERTICAL);
310         } else if(localx == elm2.second.get(HORIZONTAL)) {
311           elm2_at_x = true;
312           elm2y = elm2.second.get(VERTICAL);
313         }
314         bool retval = false;
315         if(!(elm1_at_x && elm2_at_x)) {
316           //at least one of the segments doesn't have an end point a the current x
317           //-1 below, 1 above
318           int pt1_oab = on_above_or_below(elm1.first, half_edge(elm2.first, elm2.second));
319           int pt2_oab = on_above_or_below(elm1.second, half_edge(elm2.first, elm2.second));
320           if(pt1_oab == pt2_oab) {
321             if(pt1_oab == -1)
322               retval = true; //pt1 is below elm2 so elm1 is below elm2
323           } else {
324             //the segments can't cross so elm2 is on whatever side of elm1 that one of its ends is
325             int pt3_oab = on_above_or_below(elm2.first, half_edge(elm1.first, elm1.second));
326             if(pt3_oab == 1)
327               retval = true; //elm1's point is above elm1
328           }
329         } else {
330           if(elm1y < elm2y) {
331             retval = true;
332           } else if(elm1y == elm2y) {
333             if(elm1 == elm2)
334               return false;
335             retval = less_slope(elm1.second.get(HORIZONTAL) - elm1.first.get(HORIZONTAL),
336                                      elm1.second.get(VERTICAL) - elm1.first.get(VERTICAL),
337                                      elm2.second.get(HORIZONTAL) - elm2.first.get(HORIZONTAL),
338                                      elm2.second.get(VERTICAL) - elm2.first.get(VERTICAL));
339             retval = ((*justBefore_) != 0) ^ retval;
340           }
341         }
342         return retval;
343       }
344     };
345 
346     template <typename unsigned_product_type>
unsigned_addboost::polygon::scanline_base347     static inline void unsigned_add(unsigned_product_type& result, int& result_sign, unsigned_product_type a, int a_sign, unsigned_product_type b, int b_sign) {
348       int switcher = 0;
349       if(a_sign < 0) switcher += 1;
350       if(b_sign < 0) switcher += 2;
351       if(a < b) switcher += 4;
352       switch (switcher) {
353       case 0: //both positive
354         result = a + b;
355         result_sign = 1;
356         break;
357       case 1: //a is negative
358         result = a - b;
359         result_sign = -1;
360         break;
361       case 2: //b is negative
362         result = a - b;
363         result_sign = 1;
364         break;
365       case 3: //both negative
366         result = a + b;
367         result_sign = -1;
368         break;
369       case 4: //both positive
370         result = a + b;
371         result_sign = 1;
372         break;
373       case 5: //a is negative
374         result = b - a;
375         result_sign = 1;
376         break;
377       case 6: //b is negative
378         result = b - a;
379         result_sign = -1;
380         break;
381       case 7: //both negative
382         result = b + a;
383         result_sign = -1;
384         break;
385       };
386     }
387 
388     struct compute_intersection_pack {
389       typedef typename high_precision_type<Unit>::type high_precision;
390       high_precision y_high, dx1, dy1, dx2, dy2, x11, x21, y11, y21, x_num, y_num, x_den, y_den, x, y;
compute_lazy_intersectionboost::polygon::scanline_base::compute_intersection_pack391       static inline bool compute_lazy_intersection(Point& intersection, const half_edge& he1, const half_edge& he2,
392                                                    bool projected = false, bool round_closest = false) {
393         long double y_high, dx1, dy1, dx2, dy2, x11, x21, y11, y21, x_num, y_num, x_den, y_den, x, y;
394         typedef rectangle_data<Unit> Rectangle;
395         Rectangle rect1, rect2;
396         set_points(rect1, he1.first, he1.second);
397         set_points(rect2, he2.first, he2.second);
398         if(!projected && !::boost::polygon::intersects(rect1, rect2, true)) return false;
399         if(is_vertical(he1)) {
400           if(is_vertical(he2)) return false;
401           y_high = evalAtXforYlazy(he1.first.get(HORIZONTAL), he2.first, he2.second);
402           Unit y_local = (Unit)y_high;
403           if(y_high < y_local) --y_local;
404           if(projected || contains(rect1.get(VERTICAL), y_local, true)) {
405             intersection = Point(he1.first.get(HORIZONTAL), y_local);
406             return true;
407           } else {
408             return false;
409           }
410         } else if(is_vertical(he2)) {
411           y_high = evalAtXforYlazy(he2.first.get(HORIZONTAL), he1.first, he1.second);
412           Unit y_local = (Unit)y_high;
413           if(y_high < y_local) --y_local;
414           if(projected || contains(rect2.get(VERTICAL), y_local, true)) {
415             intersection = Point(he2.first.get(HORIZONTAL), y_local);
416             return true;
417           } else {
418             return false;
419           }
420         }
421         //the bounding boxes of the two line segments intersect, so we check closer to find the intersection point
422         dy2 = (he2.second.get(VERTICAL)) -
423           (he2.first.get(VERTICAL));
424         dy1 = (he1.second.get(VERTICAL)) -
425           (he1.first.get(VERTICAL));
426         dx2 = (he2.second.get(HORIZONTAL)) -
427           (he2.first.get(HORIZONTAL));
428         dx1 = (he1.second.get(HORIZONTAL)) -
429           (he1.first.get(HORIZONTAL));
430         if(equal_slope_hp(dx1, dy1, dx2, dy2)) return false;
431         //the line segments have different slopes
432         //we can assume that the line segments are not vertical because such an intersection is handled elsewhere
433         x11 = (he1.first.get(HORIZONTAL));
434         x21 = (he2.first.get(HORIZONTAL));
435         y11 = (he1.first.get(VERTICAL));
436         y21 = (he2.first.get(VERTICAL));
437         //Unit exp_x = ((at)x11 * (at)dy1 * (at)dx2 - (at)x21 * (at)dy2 * (at)dx1 + (at)y21 * (at)dx1 * (at)dx2 - (at)y11 * (at)dx1 * (at)dx2) / ((at)dy1 * (at)dx2 - (at)dy2 * (at)dx1);
438         //Unit exp_y = ((at)y11 * (at)dx1 * (at)dy2 - (at)y21 * (at)dx2 * (at)dy1 + (at)x21 * (at)dy1 * (at)dy2 - (at)x11 * (at)dy1 * (at)dy2) / ((at)dx1 * (at)dy2 - (at)dx2 * (at)dy1);
439         x_num = (x11 * dy1 * dx2 - x21 * dy2 * dx1 + y21 * dx1 * dx2 - y11 * dx1 * dx2);
440         x_den = (dy1 * dx2 - dy2 * dx1);
441         y_num = (y11 * dx1 * dy2 - y21 * dx2 * dy1 + x21 * dy1 * dy2 - x11 * dy1 * dy2);
442         y_den = (dx1 * dy2 - dx2 * dy1);
443         x = x_num / x_den;
444         y = y_num / y_den;
445         //std::cout << "cross1 " << dy1 << " " << dx2 << " " << dy1 * dx2 << "\n";
446         //std::cout << "cross2 " << dy2 << " " << dx1 << " " << dy2 * dx1 << "\n";
447         //Unit exp_x = compute_x_intercept<at>(x11, x21, y11, y21, dy1, dy2, dx1, dx2);
448         //Unit exp_y = compute_x_intercept<at>(y11, y21, x11, x21, dx1, dx2, dy1, dy2);
449         if(round_closest) {
450           x = x + 0.5;
451           y = y + 0.5;
452         }
453         Unit x_unit = (Unit)(x);
454         Unit y_unit = (Unit)(y);
455         //truncate downward if it went up due to negative number
456         if(x < x_unit) --x_unit;
457         if(y < y_unit) --y_unit;
458         if(is_horizontal(he1))
459           y_unit = he1.first.y();
460         if(is_horizontal(he2))
461           y_unit = he2.first.y();
462         //if(x != exp_x || y != exp_y)
463         //  std::cout << exp_x << " " << exp_y << " " << x << " " << y << "\n";
464         //Unit y1 = evalAtXforY(exp_x, he1.first, he1.second);
465         //Unit y2 = evalAtXforY(exp_x, he2.first, he2.second);
466         //std::cout << exp_x << " " << exp_y << " " << y1 << " " << y2 << "\n";
467         Point result(x_unit, y_unit);
468         if(!projected && !contains(rect1, result, true)) return false;
469         if(!projected && !contains(rect2, result, true)) return false;
470         if(projected) {
471           rectangle_data<long double> inf_rect(-(long double)(std::numeric_limits<Unit>::max)(),
472                                                -(long double) (std::numeric_limits<Unit>::max)(),
473                                                (long double)(std::numeric_limits<Unit>::max)(),
474                                                (long double) (std::numeric_limits<Unit>::max)() );
475           if(contains(inf_rect, point_data<long double>(x, y), true)) {
476             intersection = result;
477             return true;
478           } else
479             return false;
480         }
481         intersection = result;
482         return true;
483       }
484 
compute_intersectionboost::polygon::scanline_base::compute_intersection_pack485       inline bool compute_intersection(Point& intersection, const half_edge& he1, const half_edge& he2,
486                                        bool projected = false, bool round_closest = false) {
487         if(!projected && !intersects(he1, he2))
488            return false;
489         bool lazy_success = compute_lazy_intersection(intersection, he1, he2, projected);
490         if(!projected) {
491           if(lazy_success) {
492             if(intersects_grid(intersection, he1) &&
493                intersects_grid(intersection, he2))
494               return true;
495           }
496         } else {
497           return lazy_success;
498         }
499         return compute_exact_intersection(intersection, he1, he2, projected, round_closest);
500       }
501 
compute_exact_intersectionboost::polygon::scanline_base::compute_intersection_pack502       inline bool compute_exact_intersection(Point& intersection, const half_edge& he1, const half_edge& he2,
503                                              bool projected = false, bool round_closest = false) {
504         if(!projected && !intersects(he1, he2))
505            return false;
506         typedef rectangle_data<Unit> Rectangle;
507         Rectangle rect1, rect2;
508         set_points(rect1, he1.first, he1.second);
509         set_points(rect2, he2.first, he2.second);
510         if(!::boost::polygon::intersects(rect1, rect2, true)) return false;
511         if(is_vertical(he1)) {
512           if(is_vertical(he2)) return false;
513           y_high = evalAtXforY(he1.first.get(HORIZONTAL), he2.first, he2.second);
514           Unit y = convert_high_precision_type<Unit>(y_high);
515           if(y_high < (high_precision)y) --y;
516           if(contains(rect1.get(VERTICAL), y, true)) {
517             intersection = Point(he1.first.get(HORIZONTAL), y);
518             return true;
519           } else {
520             return false;
521           }
522         } else if(is_vertical(he2)) {
523           y_high = evalAtXforY(he2.first.get(HORIZONTAL), he1.first, he1.second);
524           Unit y = convert_high_precision_type<Unit>(y_high);
525           if(y_high < (high_precision)y) --y;
526           if(contains(rect2.get(VERTICAL), y, true)) {
527             intersection = Point(he2.first.get(HORIZONTAL), y);
528             return true;
529           } else {
530             return false;
531           }
532         }
533         //the bounding boxes of the two line segments intersect, so we check closer to find the intersection point
534         dy2 = (high_precision)(he2.second.get(VERTICAL)) -
535           (high_precision)(he2.first.get(VERTICAL));
536         dy1 = (high_precision)(he1.second.get(VERTICAL)) -
537           (high_precision)(he1.first.get(VERTICAL));
538         dx2 = (high_precision)(he2.second.get(HORIZONTAL)) -
539           (high_precision)(he2.first.get(HORIZONTAL));
540         dx1 = (high_precision)(he1.second.get(HORIZONTAL)) -
541           (high_precision)(he1.first.get(HORIZONTAL));
542         if(equal_slope_hp(dx1, dy1, dx2, dy2)) return false;
543         //the line segments have different slopes
544         //we can assume that the line segments are not vertical because such an intersection is handled elsewhere
545         x11 = (high_precision)(he1.first.get(HORIZONTAL));
546         x21 = (high_precision)(he2.first.get(HORIZONTAL));
547         y11 = (high_precision)(he1.first.get(VERTICAL));
548         y21 = (high_precision)(he2.first.get(VERTICAL));
549         //Unit exp_x = ((at)x11 * (at)dy1 * (at)dx2 - (at)x21 * (at)dy2 * (at)dx1 + (at)y21 * (at)dx1 * (at)dx2 - (at)y11 * (at)dx1 * (at)dx2) / ((at)dy1 * (at)dx2 - (at)dy2 * (at)dx1);
550         //Unit exp_y = ((at)y11 * (at)dx1 * (at)dy2 - (at)y21 * (at)dx2 * (at)dy1 + (at)x21 * (at)dy1 * (at)dy2 - (at)x11 * (at)dy1 * (at)dy2) / ((at)dx1 * (at)dy2 - (at)dx2 * (at)dy1);
551         x_num = (x11 * dy1 * dx2 - x21 * dy2 * dx1 + y21 * dx1 * dx2 - y11 * dx1 * dx2);
552         x_den = (dy1 * dx2 - dy2 * dx1);
553         y_num = (y11 * dx1 * dy2 - y21 * dx2 * dy1 + x21 * dy1 * dy2 - x11 * dy1 * dy2);
554         y_den = (dx1 * dy2 - dx2 * dy1);
555         x = x_num / x_den;
556         y = y_num / y_den;
557         //std::cout << x << " " << y << "\n";
558         //std::cout << "cross1 " << dy1 << " " << dx2 << " " << dy1 * dx2 << "\n";
559         //std::cout << "cross2 " << dy2 << " " << dx1 << " " << dy2 * dx1 << "\n";
560         //Unit exp_x = compute_x_intercept<at>(x11, x21, y11, y21, dy1, dy2, dx1, dx2);
561         //Unit exp_y = compute_x_intercept<at>(y11, y21, x11, x21, dx1, dx2, dy1, dy2);
562         if(round_closest) {
563           x = x + (high_precision)0.5;
564           y = y + (high_precision)0.5;
565         }
566         Unit x_unit = convert_high_precision_type<Unit>(x);
567         Unit y_unit = convert_high_precision_type<Unit>(y);
568         //truncate downward if it went up due to negative number
569         if(x < (high_precision)x_unit) --x_unit;
570         if(y < (high_precision)y_unit) --y_unit;
571         if(is_horizontal(he1))
572           y_unit = he1.first.y();
573         if(is_horizontal(he2))
574           y_unit = he2.first.y();
575         //if(x != exp_x || y != exp_y)
576         //  std::cout << exp_x << " " << exp_y << " " << x << " " << y << "\n";
577         //Unit y1 = evalAtXforY(exp_x, he1.first, he1.second);
578         //Unit y2 = evalAtXforY(exp_x, he2.first, he2.second);
579         //std::cout << exp_x << " " << exp_y << " " << y1 << " " << y2 << "\n";
580         Point result(x_unit, y_unit);
581         if(!contains(rect1, result, true)) return false;
582         if(!contains(rect2, result, true)) return false;
583         if(projected) {
584           high_precision b1 = (high_precision) (std::numeric_limits<Unit>::min)();
585           high_precision b2 = (high_precision) (std::numeric_limits<Unit>::max)();
586           if(x > b2 || y > b2 || x < b1 || y < b1)
587             return false;
588         }
589         intersection = result;
590         return true;
591       }
592     };
593 
compute_intersectionboost::polygon::scanline_base594     static inline bool compute_intersection(Point& intersection, const half_edge& he1, const half_edge& he2) {
595       typedef typename high_precision_type<Unit>::type high_precision;
596       typedef rectangle_data<Unit> Rectangle;
597       Rectangle rect1, rect2;
598       set_points(rect1, he1.first, he1.second);
599       set_points(rect2, he2.first, he2.second);
600       if(!::boost::polygon::intersects(rect1, rect2, true)) return false;
601       if(is_vertical(he1)) {
602         if(is_vertical(he2)) return false;
603         high_precision y_high = evalAtXforY(he1.first.get(HORIZONTAL), he2.first, he2.second);
604         Unit y = convert_high_precision_type<Unit>(y_high);
605         if(y_high < (high_precision)y) --y;
606         if(contains(rect1.get(VERTICAL), y, true)) {
607           intersection = Point(he1.first.get(HORIZONTAL), y);
608           return true;
609         } else {
610           return false;
611         }
612       } else if(is_vertical(he2)) {
613         high_precision y_high = evalAtXforY(he2.first.get(HORIZONTAL), he1.first, he1.second);
614         Unit y = convert_high_precision_type<Unit>(y_high);
615         if(y_high < (high_precision)y) --y;
616         if(contains(rect2.get(VERTICAL), y, true)) {
617           intersection = Point(he2.first.get(HORIZONTAL), y);
618           return true;
619         } else {
620           return false;
621         }
622       }
623       //the bounding boxes of the two line segments intersect, so we check closer to find the intersection point
624       high_precision dy2 = (high_precision)(he2.second.get(VERTICAL)) -
625         (high_precision)(he2.first.get(VERTICAL));
626       high_precision dy1 = (high_precision)(he1.second.get(VERTICAL)) -
627         (high_precision)(he1.first.get(VERTICAL));
628       high_precision dx2 = (high_precision)(he2.second.get(HORIZONTAL)) -
629         (high_precision)(he2.first.get(HORIZONTAL));
630       high_precision dx1 = (high_precision)(he1.second.get(HORIZONTAL)) -
631         (high_precision)(he1.first.get(HORIZONTAL));
632       if(equal_slope_hp(dx1, dy1, dx2, dy2)) return false;
633       //the line segments have different slopes
634       //we can assume that the line segments are not vertical because such an intersection is handled elsewhere
635       high_precision x11 = (high_precision)(he1.first.get(HORIZONTAL));
636       high_precision x21 = (high_precision)(he2.first.get(HORIZONTAL));
637       high_precision y11 = (high_precision)(he1.first.get(VERTICAL));
638       high_precision y21 = (high_precision)(he2.first.get(VERTICAL));
639       //Unit exp_x = ((at)x11 * (at)dy1 * (at)dx2 - (at)x21 * (at)dy2 * (at)dx1 + (at)y21 * (at)dx1 * (at)dx2 - (at)y11 * (at)dx1 * (at)dx2) / ((at)dy1 * (at)dx2 - (at)dy2 * (at)dx1);
640       //Unit exp_y = ((at)y11 * (at)dx1 * (at)dy2 - (at)y21 * (at)dx2 * (at)dy1 + (at)x21 * (at)dy1 * (at)dy2 - (at)x11 * (at)dy1 * (at)dy2) / ((at)dx1 * (at)dy2 - (at)dx2 * (at)dy1);
641       high_precision x_num = (x11 * dy1 * dx2 - x21 * dy2 * dx1 + y21 * dx1 * dx2 - y11 * dx1 * dx2);
642       high_precision x_den = (dy1 * dx2 - dy2 * dx1);
643       high_precision y_num = (y11 * dx1 * dy2 - y21 * dx2 * dy1 + x21 * dy1 * dy2 - x11 * dy1 * dy2);
644       high_precision y_den = (dx1 * dy2 - dx2 * dy1);
645       high_precision x = x_num / x_den;
646       high_precision y = y_num / y_den;
647       //std::cout << "cross1 " << dy1 << " " << dx2 << " " << dy1 * dx2 << "\n";
648       //std::cout << "cross2 " << dy2 << " " << dx1 << " " << dy2 * dx1 << "\n";
649       //Unit exp_x = compute_x_intercept<at>(x11, x21, y11, y21, dy1, dy2, dx1, dx2);
650       //Unit exp_y = compute_x_intercept<at>(y11, y21, x11, x21, dx1, dx2, dy1, dy2);
651       Unit x_unit = convert_high_precision_type<Unit>(x);
652       Unit y_unit = convert_high_precision_type<Unit>(y);
653       //truncate downward if it went up due to negative number
654       if(x < (high_precision)x_unit) --x_unit;
655       if(y < (high_precision)y_unit) --y_unit;
656       if(is_horizontal(he1))
657         y_unit = he1.first.y();
658       if(is_horizontal(he2))
659         y_unit = he2.first.y();
660       //if(x != exp_x || y != exp_y)
661       //  std::cout << exp_x << " " << exp_y << " " << x << " " << y << "\n";
662       //Unit y1 = evalAtXforY(exp_x, he1.first, he1.second);
663       //Unit y2 = evalAtXforY(exp_x, he2.first, he2.second);
664       //std::cout << exp_x << " " << exp_y << " " << y1 << " " << y2 << "\n";
665       Point result(x_unit, y_unit);
666       if(!contains(rect1, result, true)) return false;
667       if(!contains(rect2, result, true)) return false;
668       intersection = result;
669       return true;
670     }
671 
intersectsboost::polygon::scanline_base672     static inline bool intersects(const half_edge& he1, const half_edge& he2) {
673       typedef rectangle_data<Unit> Rectangle;
674       Rectangle rect1, rect2;
675       set_points(rect1, he1.first, he1.second);
676       set_points(rect2, he2.first, he2.second);
677       if(::boost::polygon::intersects(rect1, rect2, false)) {
678         if(he1.first == he2.first) {
679           if(he1.second != he2.second && equal_slope(he1.first.get(HORIZONTAL), he1.first.get(VERTICAL),
680                                                      he1.second, he2.second)) {
681             return true;
682           } else {
683             return false;
684           }
685         }
686         if(he1.first == he2.second) {
687           if(he1.second != he2.first && equal_slope(he1.first.get(HORIZONTAL), he1.first.get(VERTICAL),
688                                                     he1.second, he2.first)) {
689             return true;
690           } else {
691             return false;
692           }
693         }
694         if(he1.second == he2.first) {
695           if(he1.first != he2.second && equal_slope(he1.second.get(HORIZONTAL), he1.second.get(VERTICAL),
696                                                     he1.first, he2.second)) {
697             return true;
698           } else {
699             return false;
700           }
701         }
702         if(he1.second == he2.second) {
703           if(he1.first != he2.first && equal_slope(he1.second.get(HORIZONTAL), he1.second.get(VERTICAL),
704                                                    he1.first, he2.first)) {
705             return true;
706           } else {
707             return false;
708           }
709         }
710         int oab1 = on_above_or_below(he1.first, he2);
711         if(oab1 == 0 && between(he1.first, he2.first, he2.second)) return true;
712         int oab2 = on_above_or_below(he1.second, he2);
713         if(oab2 == 0 && between(he1.second, he2.first, he2.second)) return true;
714         if(oab1 == oab2 && oab1 != 0) return false; //both points of he1 are on same side of he2
715         int oab3 = on_above_or_below(he2.first, he1);
716         if(oab3 == 0 && between(he2.first, he1.first, he1.second)) return true;
717         int oab4 = on_above_or_below(he2.second, he1);
718         if(oab4 == 0 && between(he2.second, he1.first, he1.second)) return true;
719         if(oab3 == oab4) return false; //both points of he2 are on same side of he1
720         return true; //they must cross
721       }
722       if(is_vertical(he1) && is_vertical(he2) && he1.first.get(HORIZONTAL) == he2.first.get(HORIZONTAL))
723         return ::boost::polygon::intersects(rect1.get(VERTICAL), rect2.get(VERTICAL), false) &&
724           rect1.get(VERTICAL) != rect2.get(VERTICAL);
725       if(is_horizontal(he1) && is_horizontal(he2) && he1.first.get(VERTICAL) == he2.first.get(VERTICAL))
726         return ::boost::polygon::intersects(rect1.get(HORIZONTAL), rect2.get(HORIZONTAL), false) &&
727           rect1.get(HORIZONTAL) != rect2.get(HORIZONTAL);
728       return false;
729     }
730 
731     class vertex_half_edge {
732     public:
733       typedef typename high_precision_type<Unit>::type high_precision;
734       Point pt;
735       Point other_pt; // 1, 0 or -1
736       int count; //dxdydTheta
vertex_half_edge()737       inline vertex_half_edge() : pt(), other_pt(), count() {}
vertex_half_edge(const Point & point,const Point & other_point,int countIn)738       inline vertex_half_edge(const Point& point, const Point& other_point, int countIn) : pt(point), other_pt(other_point), count(countIn) {}
vertex_half_edge(const vertex_half_edge & vertex)739       inline vertex_half_edge(const vertex_half_edge& vertex) : pt(vertex.pt), other_pt(vertex.other_pt), count(vertex.count) {}
operator =(const vertex_half_edge & vertex)740       inline vertex_half_edge& operator=(const vertex_half_edge& vertex){
741         pt = vertex.pt; other_pt = vertex.other_pt; count = vertex.count; return *this; }
operator ==(const vertex_half_edge & vertex) const742       inline bool operator==(const vertex_half_edge& vertex) const {
743         return pt == vertex.pt && other_pt == vertex.other_pt && count == vertex.count; }
operator !=(const vertex_half_edge & vertex) const744       inline bool operator!=(const vertex_half_edge& vertex) const { return !((*this) == vertex); }
operator <(const vertex_half_edge & vertex) const745       inline bool operator<(const vertex_half_edge& vertex) const {
746         if(pt.get(HORIZONTAL) < vertex.pt.get(HORIZONTAL)) return true;
747         if(pt.get(HORIZONTAL) == vertex.pt.get(HORIZONTAL)) {
748           if(pt.get(VERTICAL) < vertex.pt.get(VERTICAL)) return true;
749           if(pt.get(VERTICAL) == vertex.pt.get(VERTICAL)) { return less_slope(pt.get(HORIZONTAL), pt.get(VERTICAL),
750                                                                               other_pt, vertex.other_pt);
751           }
752         }
753         return false;
754       }
operator >(const vertex_half_edge & vertex) const755       inline bool operator>(const vertex_half_edge& vertex) const { return vertex < (*this); }
operator <=(const vertex_half_edge & vertex) const756       inline bool operator<=(const vertex_half_edge& vertex) const { return !((*this) > vertex); }
operator >=(const vertex_half_edge & vertex) const757       inline bool operator>=(const vertex_half_edge& vertex) const { return !((*this) < vertex); }
evalAtX(Unit xIn) const758       inline high_precision evalAtX(Unit xIn) const { return evalAtXforYlazy(xIn, pt, other_pt); }
is_vertical() const759       inline bool is_vertical() const {
760         return pt.get(HORIZONTAL) == other_pt.get(HORIZONTAL);
761       }
is_begin() const762       inline bool is_begin() const {
763         return pt.get(HORIZONTAL) < other_pt.get(HORIZONTAL) ||
764           (pt.get(HORIZONTAL) == other_pt.get(HORIZONTAL) &&
765            (pt.get(VERTICAL) < other_pt.get(VERTICAL)));
766       }
767     };
768 
769     //when scanning Vertex45 for polygon formation we need a scanline comparator functor
770     class less_vertex_half_edge {
771     private:
772       Unit *x_; //x value at which to apply comparison
773       int *justBefore_;
774     public:
775       typedef vertex_half_edge first_argument_type;
776       typedef vertex_half_edge second_argument_type;
777       typedef bool result_type;
less_vertex_half_edge()778       inline less_vertex_half_edge() : x_(0), justBefore_(0) {}
less_vertex_half_edge(Unit * x,int * justBefore)779       inline less_vertex_half_edge(Unit *x, int *justBefore) : x_(x), justBefore_(justBefore) {}
less_vertex_half_edge(const less_vertex_half_edge & that)780       inline less_vertex_half_edge(const less_vertex_half_edge& that) : x_(that.x_), justBefore_(that.justBefore_) {}
operator =(const less_vertex_half_edge & that)781       inline less_vertex_half_edge& operator=(const less_vertex_half_edge& that) { x_ = that.x_; justBefore_ = that.justBefore_; return *this; }
operator ()(const vertex_half_edge & elm1,const vertex_half_edge & elm2) const782       inline bool operator () (const vertex_half_edge& elm1, const vertex_half_edge& elm2) const {
783         if((std::max)(elm1.pt.y(), elm1.other_pt.y()) < (std::min)(elm2.pt.y(), elm2.other_pt.y()))
784           return true;
785         if((std::min)(elm1.pt.y(), elm1.other_pt.y()) > (std::max)(elm2.pt.y(), elm2.other_pt.y()))
786           return false;
787         //check if either x of elem1 is equal to x_
788         Unit localx = *x_;
789         Unit elm1y = 0;
790         bool elm1_at_x = false;
791         if(localx == elm1.pt.get(HORIZONTAL)) {
792           elm1_at_x = true;
793           elm1y = elm1.pt.get(VERTICAL);
794         } else if(localx == elm1.other_pt.get(HORIZONTAL)) {
795           elm1_at_x = true;
796           elm1y = elm1.other_pt.get(VERTICAL);
797         }
798         Unit elm2y = 0;
799         bool elm2_at_x = false;
800         if(localx == elm2.pt.get(HORIZONTAL)) {
801           elm2_at_x = true;
802           elm2y = elm2.pt.get(VERTICAL);
803         } else if(localx == elm2.other_pt.get(HORIZONTAL)) {
804           elm2_at_x = true;
805           elm2y = elm2.other_pt.get(VERTICAL);
806         }
807         bool retval = false;
808         if(!(elm1_at_x && elm2_at_x)) {
809           //at least one of the segments doesn't have an end point a the current x
810           //-1 below, 1 above
811           int pt1_oab = on_above_or_below(elm1.pt, half_edge(elm2.pt, elm2.other_pt));
812           int pt2_oab = on_above_or_below(elm1.other_pt, half_edge(elm2.pt, elm2.other_pt));
813           if(pt1_oab == pt2_oab) {
814             if(pt1_oab == -1)
815               retval = true; //pt1 is below elm2 so elm1 is below elm2
816           } else {
817             //the segments can't cross so elm2 is on whatever side of elm1 that one of its ends is
818             int pt3_oab = on_above_or_below(elm2.pt, half_edge(elm1.pt, elm1.other_pt));
819             if(pt3_oab == 1)
820               retval = true; //elm1's point is above elm1
821           }
822         } else {
823           if(elm1y < elm2y) {
824             retval = true;
825           } else if(elm1y == elm2y) {
826             if(elm1.pt == elm2.pt && elm1.other_pt == elm2.other_pt)
827               return false;
828             retval = less_slope(elm1.other_pt.get(HORIZONTAL) - elm1.pt.get(HORIZONTAL),
829                                      elm1.other_pt.get(VERTICAL) - elm1.pt.get(VERTICAL),
830                                      elm2.other_pt.get(HORIZONTAL) - elm2.pt.get(HORIZONTAL),
831                                      elm2.other_pt.get(VERTICAL) - elm2.pt.get(VERTICAL));
832             retval = ((*justBefore_) != 0) ^ retval;
833           }
834         }
835         return retval;
836       }
837     };
838 
839   };
840 
841   template <typename Unit>
842   class polygon_arbitrary_formation : public scanline_base<Unit> {
843   public:
844     typedef typename scanline_base<Unit>::Point Point;
845     typedef typename scanline_base<Unit>::half_edge half_edge;
846     typedef typename scanline_base<Unit>::vertex_half_edge vertex_half_edge;
847     typedef typename scanline_base<Unit>::less_vertex_half_edge less_vertex_half_edge;
848 
849     class poly_line_arbitrary {
850     public:
851       typedef typename std::list<Point>::const_iterator iterator;
852 
853       // default constructor of point does not initialize x and y
poly_line_arbitrary()854       inline poly_line_arbitrary() : points() {} //do nothing default constructor
855 
856       // initialize a polygon from x,y values, it is assumed that the first is an x
857       // and that the input is a well behaved polygon
858       template<class iT>
set(iT inputBegin,iT inputEnd)859       inline poly_line_arbitrary& set(iT inputBegin, iT inputEnd) {
860         points.clear();  //just in case there was some old data there
861         while(inputBegin != inputEnd) {
862           points.insert(points.end(), *inputBegin);
863           ++inputBegin;
864         }
865         return *this;
866       }
867 
868       // copy constructor (since we have dynamic memory)
poly_line_arbitrary(const poly_line_arbitrary & that)869       inline poly_line_arbitrary(const poly_line_arbitrary& that) : points(that.points) {}
870 
871       // assignment operator (since we have dynamic memory do a deep copy)
operator =(const poly_line_arbitrary & that)872       inline poly_line_arbitrary& operator=(const poly_line_arbitrary& that) {
873         points = that.points;
874         return *this;
875       }
876 
877       // get begin iterator, returns a pointer to a const Unit
begin() const878       inline iterator begin() const { return points.begin(); }
879 
880       // get end iterator, returns a pointer to a const Unit
end() const881       inline iterator end() const { return points.end(); }
882 
size() const883       inline std::size_t size() const { return points.size(); }
884 
885       //public data member
886       std::list<Point> points;
887     };
888 
889     class active_tail_arbitrary {
890     protected:
891       //data
892       poly_line_arbitrary* tailp_;
893       active_tail_arbitrary *otherTailp_;
894       std::list<active_tail_arbitrary*> holesList_;
895       bool head_;
896     public:
897 
898       /**
899        * @brief iterator over coordinates of the figure
900        */
901       typedef typename poly_line_arbitrary::iterator iterator;
902 
903       /**
904        * @brief iterator over holes contained within the figure
905        */
906       typedef typename std::list<active_tail_arbitrary*>::const_iterator iteratorHoles;
907 
908       //default constructor
active_tail_arbitrary()909       inline active_tail_arbitrary() : tailp_(), otherTailp_(), holesList_(), head_() {}
910 
911       //constructor
active_tail_arbitrary(const vertex_half_edge & vertex,active_tail_arbitrary * otherTailp=0)912       inline active_tail_arbitrary(const vertex_half_edge& vertex, active_tail_arbitrary* otherTailp = 0) : tailp_(), otherTailp_(), holesList_(), head_() {
913         tailp_ = new poly_line_arbitrary;
914         tailp_->points.push_back(vertex.pt);
915         //bool headArray[4] = {false, true, true, true};
916         bool inverted = vertex.count == -1;
917         head_ = (!vertex.is_vertical) ^ inverted;
918         otherTailp_ = otherTailp;
919       }
920 
active_tail_arbitrary(Point point,active_tail_arbitrary * otherTailp,bool head=true)921       inline active_tail_arbitrary(Point point, active_tail_arbitrary* otherTailp, bool head = true) :
922         tailp_(), otherTailp_(), holesList_(), head_() {
923         tailp_ = new poly_line_arbitrary;
924         tailp_->points.push_back(point);
925         head_ = head;
926         otherTailp_ = otherTailp;
927 
928       }
active_tail_arbitrary(active_tail_arbitrary * otherTailp)929       inline active_tail_arbitrary(active_tail_arbitrary* otherTailp) :
930         tailp_(), otherTailp_(), holesList_(), head_() {
931         tailp_ = otherTailp->tailp_;
932         otherTailp_ = otherTailp;
933       }
934 
935       //copy constructor
active_tail_arbitrary(const active_tail_arbitrary & that)936       inline active_tail_arbitrary(const active_tail_arbitrary& that) :
937         tailp_(), otherTailp_(), holesList_(), head_() { (*this) = that; }
938 
939       //destructor
~active_tail_arbitrary()940       inline ~active_tail_arbitrary() {
941         destroyContents();
942       }
943 
944       //assignment operator
operator =(const active_tail_arbitrary & that)945       inline active_tail_arbitrary& operator=(const active_tail_arbitrary& that) {
946         tailp_ = new poly_line_arbitrary(*(that.tailp_));
947         head_ = that.head_;
948         otherTailp_ = that.otherTailp_;
949         holesList_ = that.holesList_;
950         return *this;
951       }
952 
953       //equivalence operator
operator ==(const active_tail_arbitrary & b) const954       inline bool operator==(const active_tail_arbitrary& b) const {
955         return tailp_ == b.tailp_ && head_ == b.head_;
956       }
957 
958       /**
959        * @brief get the pointer to the polyline that this is an active tail of
960        */
getTail() const961       inline poly_line_arbitrary* getTail() const { return tailp_; }
962 
963       /**
964        * @brief get the pointer to the polyline at the other end of the chain
965        */
getOtherTail() const966       inline poly_line_arbitrary* getOtherTail() const { return otherTailp_->tailp_; }
967 
968       /**
969        * @brief get the pointer to the activetail at the other end of the chain
970        */
getOtherActiveTail() const971       inline active_tail_arbitrary* getOtherActiveTail() const { return otherTailp_; }
972 
973       /**
974        * @brief test if another active tail is the other end of the chain
975        */
isOtherTail(const active_tail_arbitrary & b) const976       inline bool isOtherTail(const active_tail_arbitrary& b) const { return &b == otherTailp_; }
977 
978       /**
979        * @brief update this end of chain pointer to new polyline
980        */
updateTail(poly_line_arbitrary * newTail)981       inline active_tail_arbitrary& updateTail(poly_line_arbitrary* newTail) { tailp_ = newTail; return *this; }
982 
join(active_tail_arbitrary * tail)983       inline bool join(active_tail_arbitrary* tail) {
984         if(tail == otherTailp_) {
985           //std::cout << "joining to other tail!\n";
986           return false;
987         }
988         if(tail->head_ == head_) {
989           //std::cout << "joining head to head!\n";
990           return false;
991         }
992         if(!tailp_) {
993           //std::cout << "joining empty tail!\n";
994           return false;
995         }
996         if(!(otherTailp_->head_)) {
997           otherTailp_->copyHoles(*tail);
998           otherTailp_->copyHoles(*this);
999         } else {
1000           tail->otherTailp_->copyHoles(*this);
1001           tail->otherTailp_->copyHoles(*tail);
1002         }
1003         poly_line_arbitrary* tail1 = tailp_;
1004         poly_line_arbitrary* tail2 = tail->tailp_;
1005         if(head_) std::swap(tail1, tail2);
1006         typename std::list<point_data<Unit> >::reverse_iterator riter = tail1->points.rbegin();
1007         typename std::list<point_data<Unit> >::iterator iter = tail2->points.begin();
1008         if(*riter == *iter) {
1009           tail1->points.pop_back(); //remove duplicate point
1010         }
1011         tail1->points.splice(tail1->points.end(), tail2->points);
1012         delete tail2;
1013         otherTailp_->tailp_ = tail1;
1014         tail->otherTailp_->tailp_ = tail1;
1015         otherTailp_->otherTailp_ = tail->otherTailp_;
1016         tail->otherTailp_->otherTailp_ = otherTailp_;
1017         tailp_ = 0;
1018         tail->tailp_ = 0;
1019         tail->otherTailp_ = 0;
1020         otherTailp_ = 0;
1021         return true;
1022       }
1023 
1024       /**
1025        * @brief associate a hole to this active tail by the specified policy
1026        */
addHole(active_tail_arbitrary * hole)1027       inline active_tail_arbitrary* addHole(active_tail_arbitrary* hole) {
1028         holesList_.push_back(hole);
1029         copyHoles(*hole);
1030         copyHoles(*(hole->otherTailp_));
1031         return this;
1032       }
1033 
1034       /**
1035        * @brief get the list of holes
1036        */
getHoles() const1037       inline const std::list<active_tail_arbitrary*>& getHoles() const { return holesList_; }
1038 
1039       /**
1040        * @brief copy holes from that to this
1041        */
copyHoles(active_tail_arbitrary & that)1042       inline void copyHoles(active_tail_arbitrary& that) { holesList_.splice(holesList_.end(), that.holesList_); }
1043 
1044       /**
1045        * @brief find out if solid to right
1046        */
solidToRight() const1047       inline bool solidToRight() const { return !head_; }
solidToLeft() const1048       inline bool solidToLeft() const { return head_; }
1049 
1050       /**
1051        * @brief get vertex
1052        */
getPoint() const1053       inline Point getPoint() const {
1054         if(head_) return tailp_->points.front();
1055         return tailp_->points.back();
1056       }
1057 
1058       /**
1059        * @brief add a coordinate to the polygon at this active tail end, properly handle degenerate edges by removing redundant coordinate
1060        */
pushPoint(Point point)1061       inline void pushPoint(Point point) {
1062         if(head_) {
1063           //if(tailp_->points.size() < 2) {
1064           //   tailp_->points.push_front(point);
1065           //   return;
1066           //}
1067           typename std::list<Point>::iterator iter = tailp_->points.begin();
1068           if(iter == tailp_->points.end()) {
1069             tailp_->points.push_front(point);
1070             return;
1071           }
1072           ++iter;
1073           if(iter == tailp_->points.end()) {
1074             tailp_->points.push_front(point);
1075             return;
1076           }
1077           --iter;
1078           if(*iter != point) {
1079             tailp_->points.push_front(point);
1080           }
1081           return;
1082         }
1083         //if(tailp_->points.size() < 2) {
1084         //   tailp_->points.push_back(point);
1085         //   return;
1086         //}
1087         typename std::list<Point>::reverse_iterator iter = tailp_->points.rbegin();
1088         if(iter == tailp_->points.rend()) {
1089           tailp_->points.push_back(point);
1090           return;
1091         }
1092         ++iter;
1093         if(iter == tailp_->points.rend()) {
1094           tailp_->points.push_back(point);
1095           return;
1096         }
1097         --iter;
1098         if(*iter != point) {
1099           tailp_->points.push_back(point);
1100         }
1101       }
1102 
1103       /**
1104        * @brief joins the two chains that the two active tail tails are ends of
1105        * checks for closure of figure and writes out polygons appropriately
1106        * returns a handle to a hole if one is closed
1107        */
1108       template <class cT>
joinChains(Point point,active_tail_arbitrary * at1,active_tail_arbitrary * at2,bool solid,cT & output)1109       static inline active_tail_arbitrary* joinChains(Point point, active_tail_arbitrary* at1, active_tail_arbitrary* at2, bool solid,
1110                                                       cT& output) {
1111         if(at1->otherTailp_ == at2) {
1112           //if(at2->otherTailp_ != at1) std::cout << "half closed error\n";
1113           //we are closing a figure
1114           at1->pushPoint(point);
1115           at2->pushPoint(point);
1116           if(solid) {
1117             //we are closing a solid figure, write to output
1118             //std::cout << "test1\n";
1119             at1->copyHoles(*(at1->otherTailp_));
1120             typename PolyLineArbitraryByConcept<Unit, typename geometry_concept<typename cT::value_type>::type>::type polyData(at1);
1121             //poly_line_arbitrary_polygon_data polyData(at1);
1122             //std::cout << "test2\n";
1123             //std::cout << poly << "\n";
1124             //std::cout << "test3\n";
1125             typedef typename cT::value_type result_type;
1126             output.push_back(result_type());
1127             assign(output.back(), polyData);
1128             //std::cout << "test4\n";
1129             //std::cout << "delete " << at1->otherTailp_ << "\n";
1130             //at1->print();
1131             //at1->otherTailp_->print();
1132             delete at1->otherTailp_;
1133             //at1->print();
1134             //at1->otherTailp_->print();
1135             //std::cout << "test5\n";
1136             //std::cout << "delete " << at1 << "\n";
1137             delete at1;
1138             //std::cout << "test6\n";
1139             return 0;
1140           } else {
1141             //we are closing a hole, return the tail end active tail of the figure
1142             return at1;
1143           }
1144         }
1145         //we are not closing a figure
1146         at1->pushPoint(point);
1147         at1->join(at2);
1148         delete at1;
1149         delete at2;
1150         return 0;
1151       }
1152 
destroyContents()1153       inline void destroyContents() {
1154         if(otherTailp_) {
1155           //std::cout << "delete p " << tailp_ << "\n";
1156           if(tailp_) delete tailp_;
1157           tailp_ = 0;
1158           otherTailp_->otherTailp_ = 0;
1159           otherTailp_->tailp_ = 0;
1160           otherTailp_ = 0;
1161         }
1162         for(typename std::list<active_tail_arbitrary*>::iterator itr = holesList_.begin(); itr != holesList_.end(); ++itr) {
1163           //std::cout << "delete p " << (*itr) << "\n";
1164           if(*itr) {
1165             if((*itr)->otherTailp_) {
1166               delete (*itr)->otherTailp_;
1167               (*itr)->otherTailp_ = 0;
1168             }
1169             delete (*itr);
1170           }
1171           (*itr) = 0;
1172         }
1173         holesList_.clear();
1174       }
1175 
print()1176       inline void print() {
1177         //std::cout << this << " " << tailp_ << " " << otherTailp_ << " " << holesList_.size() << " " << head_ << "\n";
1178       }
1179 
createActiveTailsAsPair(Point point,bool solid,active_tail_arbitrary * phole,bool fractureHoles)1180       static inline std::pair<active_tail_arbitrary*, active_tail_arbitrary*> createActiveTailsAsPair(Point point, bool solid,
1181                                                                                                       active_tail_arbitrary* phole, bool fractureHoles) {
1182         active_tail_arbitrary* at1 = 0;
1183         active_tail_arbitrary* at2 = 0;
1184         if(phole && fractureHoles) {
1185           //std::cout << "adding hole\n";
1186           at1 = phole;
1187           //assert solid == false, we should be creating a corner with solid below and to the left if there was a hole
1188           at2 = at1->getOtherActiveTail();
1189           at2->pushPoint(point);
1190           at1->pushPoint(point);
1191         } else {
1192           at1 = new active_tail_arbitrary(point, at2, solid);
1193           at2 = new active_tail_arbitrary(at1);
1194           at1->otherTailp_ = at2;
1195           at2->head_ = !solid;
1196           if(phole)
1197             at2->addHole(phole); //assert fractureHoles == false
1198         }
1199         return std::pair<active_tail_arbitrary*, active_tail_arbitrary*>(at1, at2);
1200       }
1201 
1202     };
1203 
1204 
1205     typedef std::vector<std::pair<Point, int> > vertex_arbitrary_count;
1206 
1207     class less_half_edge_count {
1208     private:
1209       Point pt_;
1210     public:
1211       typedef vertex_half_edge first_argument_type;
1212       typedef vertex_half_edge second_argument_type;
1213       typedef bool result_type;
less_half_edge_count()1214       inline less_half_edge_count() : pt_() {}
less_half_edge_count(Point point)1215       inline less_half_edge_count(Point point) : pt_(point) {}
operator ()(const std::pair<Point,int> & elm1,const std::pair<Point,int> & elm2) const1216       inline bool operator () (const std::pair<Point, int>& elm1, const std::pair<Point, int>& elm2) const {
1217         return scanline_base<Unit>::less_slope(pt_.get(HORIZONTAL), pt_.get(VERTICAL), elm1.first, elm2.first);
1218       }
1219     };
1220 
sort_vertex_arbitrary_count(vertex_arbitrary_count & count,const Point & pt)1221     static inline void sort_vertex_arbitrary_count(vertex_arbitrary_count& count, const Point& pt) {
1222       less_half_edge_count lfec(pt);
1223       polygon_sort(count.begin(), count.end(), lfec);
1224     }
1225 
1226     typedef std::vector<std::pair<std::pair<std::pair<Point, Point>, int>, active_tail_arbitrary*> > incoming_count;
1227 
1228     class less_incoming_count {
1229     private:
1230       Point pt_;
1231     public:
1232       typedef std::pair<std::pair<std::pair<Point, Point>, int>, active_tail_arbitrary*> first_argument_type;
1233       typedef std::pair<std::pair<std::pair<Point, Point>, int>, active_tail_arbitrary*> second_argument_type;
1234       typedef bool result_type;
less_incoming_count()1235       inline less_incoming_count() : pt_() {}
less_incoming_count(Point point)1236       inline less_incoming_count(Point point) : pt_(point) {}
operator ()(const std::pair<std::pair<std::pair<Point,Point>,int>,active_tail_arbitrary * > & elm1,const std::pair<std::pair<std::pair<Point,Point>,int>,active_tail_arbitrary * > & elm2) const1237       inline bool operator () (const std::pair<std::pair<std::pair<Point, Point>, int>, active_tail_arbitrary*>& elm1,
1238                                const std::pair<std::pair<std::pair<Point, Point>, int>, active_tail_arbitrary*>& elm2) const {
1239         Unit dx1 = elm1.first.first.first.get(HORIZONTAL) - elm1.first.first.second.get(HORIZONTAL);
1240         Unit dx2 = elm2.first.first.first.get(HORIZONTAL) - elm2.first.first.second.get(HORIZONTAL);
1241         Unit dy1 = elm1.first.first.first.get(VERTICAL) - elm1.first.first.second.get(VERTICAL);
1242         Unit dy2 = elm2.first.first.first.get(VERTICAL) - elm2.first.first.second.get(VERTICAL);
1243         return scanline_base<Unit>::less_slope(dx1, dy1, dx2, dy2);
1244       }
1245     };
1246 
sort_incoming_count(incoming_count & count,const Point & pt)1247     static inline void sort_incoming_count(incoming_count& count, const Point& pt) {
1248       less_incoming_count lfec(pt);
1249       polygon_sort(count.begin(), count.end(), lfec);
1250     }
1251 
compact_vertex_arbitrary_count(const Point & pt,vertex_arbitrary_count & count)1252     static inline void compact_vertex_arbitrary_count(const Point& pt, vertex_arbitrary_count &count) {
1253       if(count.empty()) return;
1254       vertex_arbitrary_count tmp;
1255       tmp.reserve(count.size());
1256       tmp.push_back(count[0]);
1257       //merge duplicates
1258       for(std::size_t i = 1; i < count.size(); ++i) {
1259         if(!equal_slope(pt.get(HORIZONTAL), pt.get(VERTICAL), tmp[i-1].first, count[i].first)) {
1260           tmp.push_back(count[i]);
1261         } else {
1262           tmp.back().second += count[i].second;
1263         }
1264       }
1265       count.clear();
1266       count.swap(tmp);
1267     }
1268 
1269     // inline std::ostream& operator<< (std::ostream& o, const vertex_arbitrary_count& c) {
1270 //       for(unsinged int i = 0; i < c.size(); ++i) {
1271 //         o << c[i].first << " " << c[i].second << " ";
1272 //       }
1273 //       return o;
1274 //     }
1275 
1276     class vertex_arbitrary_compact {
1277     public:
1278       Point pt;
1279       vertex_arbitrary_count count;
vertex_arbitrary_compact()1280       inline vertex_arbitrary_compact() : pt(), count() {}
vertex_arbitrary_compact(const Point & point,const Point & other_point,int countIn)1281       inline vertex_arbitrary_compact(const Point& point, const Point& other_point, int countIn) : pt(point), count() {
1282         count.push_back(std::pair<Point, int>(other_point, countIn));
1283       }
vertex_arbitrary_compact(const vertex_half_edge & vertex)1284       inline vertex_arbitrary_compact(const vertex_half_edge& vertex) : pt(vertex.pt), count() {
1285         count.push_back(std::pair<Point, int>(vertex.other_pt, vertex.count));
1286       }
vertex_arbitrary_compact(const vertex_arbitrary_compact & vertex)1287       inline vertex_arbitrary_compact(const vertex_arbitrary_compact& vertex) : pt(vertex.pt), count(vertex.count) {}
operator =(const vertex_arbitrary_compact & vertex)1288       inline vertex_arbitrary_compact& operator=(const vertex_arbitrary_compact& vertex){
1289         pt = vertex.pt; count = vertex.count; return *this; }
operator ==(const vertex_arbitrary_compact & vertex) const1290       inline bool operator==(const vertex_arbitrary_compact& vertex) const {
1291         return pt == vertex.pt && count == vertex.count; }
operator !=(const vertex_arbitrary_compact & vertex) const1292       inline bool operator!=(const vertex_arbitrary_compact& vertex) const { return !((*this) == vertex); }
operator <(const vertex_arbitrary_compact & vertex) const1293       inline bool operator<(const vertex_arbitrary_compact& vertex) const {
1294         if(pt.get(HORIZONTAL) < vertex.pt.get(HORIZONTAL)) return true;
1295         if(pt.get(HORIZONTAL) == vertex.pt.get(HORIZONTAL)) {
1296           return pt.get(VERTICAL) < vertex.pt.get(VERTICAL);
1297         }
1298         return false;
1299       }
operator >(const vertex_arbitrary_compact & vertex) const1300       inline bool operator>(const vertex_arbitrary_compact& vertex) const { return vertex < (*this); }
operator <=(const vertex_arbitrary_compact & vertex) const1301       inline bool operator<=(const vertex_arbitrary_compact& vertex) const { return !((*this) > vertex); }
operator >=(const vertex_arbitrary_compact & vertex) const1302       inline bool operator>=(const vertex_arbitrary_compact& vertex) const { return !((*this) < vertex); }
have_vertex_half_edge(int index) const1303       inline bool have_vertex_half_edge(int index) const { return count[index]; }
operator [](int index) const1304       inline vertex_half_edge operator[](int index) const { return vertex_half_edge(pt, count[index]); }
1305       };
1306 
1307 //     inline std::ostream& operator<< (std::ostream& o, const vertex_arbitrary_compact& c) {
1308 //       o << c.pt << ", " << c.count;
1309 //       return o;
1310 //     }
1311 
1312   protected:
1313     //definitions
1314     typedef std::map<vertex_half_edge, active_tail_arbitrary*, less_vertex_half_edge> scanline_data;
1315     typedef typename scanline_data::iterator iterator;
1316     typedef typename scanline_data::const_iterator const_iterator;
1317 
1318     //data
1319     scanline_data scanData_;
1320     Unit x_;
1321     int justBefore_;
1322     int fractureHoles_;
1323   public:
polygon_arbitrary_formation()1324     inline polygon_arbitrary_formation() :
1325       scanData_(), x_((std::numeric_limits<Unit>::min)()), justBefore_(false), fractureHoles_(0) {
1326       less_vertex_half_edge lessElm(&x_, &justBefore_);
1327       scanData_ = scanline_data(lessElm);
1328     }
polygon_arbitrary_formation(bool fractureHoles)1329     inline polygon_arbitrary_formation(bool fractureHoles) :
1330       scanData_(), x_((std::numeric_limits<Unit>::min)()), justBefore_(false), fractureHoles_(fractureHoles) {
1331       less_vertex_half_edge lessElm(&x_, &justBefore_);
1332       scanData_ = scanline_data(lessElm);
1333     }
polygon_arbitrary_formation(const polygon_arbitrary_formation & that)1334     inline polygon_arbitrary_formation(const polygon_arbitrary_formation& that) :
1335       scanData_(), x_((std::numeric_limits<Unit>::min)()), justBefore_(false), fractureHoles_(0) { (*this) = that; }
operator =(const polygon_arbitrary_formation & that)1336     inline polygon_arbitrary_formation& operator=(const polygon_arbitrary_formation& that) {
1337       x_ = that.x_;
1338       justBefore_ = that.justBefore_;
1339       fractureHoles_ = that.fractureHoles_;
1340       less_vertex_half_edge lessElm(&x_, &justBefore_);
1341       scanData_ = scanline_data(lessElm);
1342       for(const_iterator itr = that.scanData_.begin(); itr != that.scanData_.end(); ++itr){
1343         scanData_.insert(scanData_.end(), *itr);
1344       }
1345       return *this;
1346     }
1347 
1348     //cT is an output container of Polygon45 or Polygon45WithHoles
1349     //iT is an iterator over vertex_half_edge elements
1350     //inputBegin - inputEnd is a range of sorted iT that represents
1351     //one or more scanline stops worth of data
1352     template <class cT, class iT>
scan(cT & output,iT inputBegin,iT inputEnd)1353     void scan(cT& output, iT inputBegin, iT inputEnd) {
1354       //std::cout << "1\n";
1355       while(inputBegin != inputEnd) {
1356         //std::cout << "2\n";
1357         x_ = (*inputBegin).pt.get(HORIZONTAL);
1358         //std::cout << "SCAN FORMATION " << x_ << "\n";
1359         //std::cout << "x_ = " << x_ << "\n";
1360         //std::cout << "scan line size: " << scanData_.size() << "\n";
1361         inputBegin = processEvent_(output, inputBegin, inputEnd);
1362       }
1363       //std::cout << "scan line size: " << scanData_.size() << "\n";
1364     }
1365 
1366   protected:
1367     //functions
1368     template <class cT, class cT2>
processPoint_(cT & output,cT2 & elements,Point point,incoming_count & counts_from_scanline,vertex_arbitrary_count & incoming_count)1369     inline std::pair<std::pair<Point, int>, active_tail_arbitrary*> processPoint_(cT& output, cT2& elements, Point point,
1370                                                                                   incoming_count& counts_from_scanline, vertex_arbitrary_count& incoming_count) {
1371       //std::cout << "\nAT POINT: " <<  point << "\n";
1372       //join any closing solid corners
1373       std::vector<int> counts;
1374       std::vector<int> incoming;
1375       std::vector<active_tail_arbitrary*> tails;
1376       counts.reserve(counts_from_scanline.size());
1377       tails.reserve(counts_from_scanline.size());
1378       incoming.reserve(incoming_count.size());
1379       for(std::size_t i = 0; i < counts_from_scanline.size(); ++i) {
1380         counts.push_back(counts_from_scanline[i].first.second);
1381         tails.push_back(counts_from_scanline[i].second);
1382       }
1383       for(std::size_t i = 0; i < incoming_count.size(); ++i) {
1384         incoming.push_back(incoming_count[i].second);
1385         if(incoming_count[i].first < point) {
1386           incoming.back() = 0;
1387         }
1388       }
1389 
1390       active_tail_arbitrary* returnValue = 0;
1391       std::pair<Point, int> returnCount(Point(0, 0), 0);
1392       int i_size_less_1 = (int)(incoming.size()) -1;
1393       int c_size_less_1 = (int)(counts.size()) -1;
1394       int i_size = incoming.size();
1395       int c_size = counts.size();
1396 
1397       bool have_vertical_tail_from_below = false;
1398       if(c_size &&
1399          scanline_base<Unit>::is_vertical(counts_from_scanline.back().first.first)) {
1400         have_vertical_tail_from_below = true;
1401       }
1402       //assert size = size_less_1 + 1
1403       //std::cout << tails.size() << " " << incoming.size() << " " << counts_from_scanline.size() << " " << incoming_count.size() << "\n";
1404       //         for(std::size_t i = 0; i < counts.size(); ++i) {
1405       //           std::cout << counts_from_scanline[i].first.first.first.get(HORIZONTAL) << ",";
1406       //           std::cout << counts_from_scanline[i].first.first.first.get(VERTICAL) << " ";
1407       //           std::cout << counts_from_scanline[i].first.first.second.get(HORIZONTAL) << ",";
1408       //           std::cout << counts_from_scanline[i].first.first.second.get(VERTICAL) << ":";
1409       //           std::cout << counts_from_scanline[i].first.second << " ";
1410       //         } std::cout << "\n";
1411       //         print(incoming_count);
1412       {
1413         for(int i = 0; i < c_size_less_1; ++i) {
1414           //std::cout << i << "\n";
1415           if(counts[i] == -1) {
1416             //std::cout << "fixed i\n";
1417             for(int j = i + 1; j < c_size; ++j) {
1418               //std::cout << j << "\n";
1419               if(counts[j]) {
1420                 if(counts[j] == 1) {
1421                   //std::cout << "case1: " << i << " " << j << "\n";
1422                   //if a figure is closed it will be written out by this function to output
1423                   active_tail_arbitrary::joinChains(point, tails[i], tails[j], true, output);
1424                   counts[i] = 0;
1425                   counts[j] = 0;
1426                   tails[i] = 0;
1427                   tails[j] = 0;
1428                 }
1429                 break;
1430               }
1431             }
1432           }
1433         }
1434       }
1435       //find any pairs of incoming edges that need to create pair for leading solid
1436       //std::cout << "checking case2\n";
1437       {
1438         for(int i = 0; i < i_size_less_1; ++i) {
1439           //std::cout << i << "\n";
1440           if(incoming[i] == 1) {
1441             //std::cout << "fixed i\n";
1442             for(int j = i + 1; j < i_size; ++j) {
1443               //std::cout << j << "\n";
1444               if(incoming[j]) {
1445                 //std::cout << incoming[j] << "\n";
1446                 if(incoming[j] == -1) {
1447                   //std::cout << "case2: " << i << " " << j << "\n";
1448                   //std::cout << "creating active tail pair\n";
1449                   std::pair<active_tail_arbitrary*, active_tail_arbitrary*> tailPair =
1450                     active_tail_arbitrary::createActiveTailsAsPair(point, true, 0, fractureHoles_ != 0);
1451                   //tailPair.first->print();
1452                   //tailPair.second->print();
1453                   if(j == i_size_less_1 && incoming_count[j].first.get(HORIZONTAL) == point.get(HORIZONTAL)) {
1454                     //vertical active tail becomes return value
1455                     returnValue = tailPair.first;
1456                     returnCount.first = point;
1457                     returnCount.second = 1;
1458                   } else {
1459                     //std::cout << "new element " << j-1 << " " << -1 << "\n";
1460                     //std::cout << point << " " <<  incoming_count[j].first << "\n";
1461                     elements.push_back(std::pair<vertex_half_edge,
1462                                        active_tail_arbitrary*>(vertex_half_edge(point,
1463                                                                                 incoming_count[j].first, -1), tailPair.first));
1464                   }
1465                   //std::cout << "new element " << i-1 << " " << 1 << "\n";
1466                   //std::cout << point << " " <<  incoming_count[i].first << "\n";
1467                   elements.push_back(std::pair<vertex_half_edge,
1468                                      active_tail_arbitrary*>(vertex_half_edge(point,
1469                                                                               incoming_count[i].first, 1), tailPair.second));
1470                   incoming[i] = 0;
1471                   incoming[j] = 0;
1472                 }
1473                 break;
1474               }
1475             }
1476           }
1477         }
1478       }
1479       //find any active tail that needs to pass through to an incoming edge
1480       //we expect to find no more than two pass through
1481 
1482       //find pass through with solid on top
1483       {
1484         //std::cout << "checking case 3\n";
1485         for(int i = 0; i < c_size; ++i) {
1486           //std::cout << i << "\n";
1487           if(counts[i] != 0) {
1488             if(counts[i] == 1) {
1489               //std::cout << "fixed i\n";
1490               for(int j = i_size_less_1; j >= 0; --j) {
1491                 if(incoming[j] != 0) {
1492                   if(incoming[j] == 1) {
1493                     //std::cout << "case3: " << i << " " << j << "\n";
1494                     //tails[i]->print();
1495                     //pass through solid on top
1496                     tails[i]->pushPoint(point);
1497                     //std::cout << "after push\n";
1498                     if(j == i_size_less_1 && incoming_count[j].first.get(HORIZONTAL) == point.get(HORIZONTAL)) {
1499                       returnValue = tails[i];
1500                       returnCount.first = point;
1501                       returnCount.second = -1;
1502                     } else {
1503                       elements.push_back(std::pair<vertex_half_edge,
1504                                          active_tail_arbitrary*>(vertex_half_edge(point,
1505                                                                                   incoming_count[j].first, incoming[j]), tails[i]));
1506                     }
1507                     tails[i] = 0;
1508                     counts[i] = 0;
1509                     incoming[j] = 0;
1510                   }
1511                   break;
1512                 }
1513               }
1514             }
1515             break;
1516           }
1517         }
1518       }
1519       //std::cout << "checking case 4\n";
1520       //find pass through with solid on bottom
1521       {
1522         for(int i = c_size_less_1; i >= 0; --i) {
1523           //std::cout << "i = " << i << " with count " << counts[i] << "\n";
1524           if(counts[i] != 0) {
1525             if(counts[i] == -1) {
1526               for(int j = 0; j < i_size; ++j) {
1527                 if(incoming[j] != 0) {
1528                   if(incoming[j] == -1) {
1529                     //std::cout << "case4: " << i << " " << j << "\n";
1530                     //pass through solid on bottom
1531                     tails[i]->pushPoint(point);
1532                     if(j == i_size_less_1 && incoming_count[j].first.get(HORIZONTAL) == point.get(HORIZONTAL)) {
1533                       returnValue = tails[i];
1534                       returnCount.first = point;
1535                       returnCount.second = 1;
1536                     } else {
1537                       //std::cout << "new element " << j-1 << " " << incoming[j] << "\n";
1538                       //std::cout << point << " " <<  incoming_count[j].first << "\n";
1539                       elements.push_back(std::pair<vertex_half_edge,
1540                                          active_tail_arbitrary*>(vertex_half_edge(point,
1541                                                                                   incoming_count[j].first, incoming[j]), tails[i]));
1542                     }
1543                     tails[i] = 0;
1544                     counts[i] = 0;
1545                     incoming[j] = 0;
1546                   }
1547                   break;
1548                 }
1549               }
1550             }
1551             break;
1552           }
1553         }
1554       }
1555       //find the end of a hole or the beginning of a hole
1556 
1557       //find end of a hole
1558       {
1559         for(int i = 0; i < c_size_less_1; ++i) {
1560           if(counts[i] != 0) {
1561             for(int j = i+1; j < c_size; ++j) {
1562               if(counts[j] != 0) {
1563                 //std::cout << "case5: " << i << " " << j << "\n";
1564                 //we are ending a hole and may potentially close a figure and have to handle the hole
1565                 returnValue = active_tail_arbitrary::joinChains(point, tails[i], tails[j], false, output);
1566                 if(returnValue) returnCount.first = point;
1567                 //std::cout << returnValue << "\n";
1568                 tails[i] = 0;
1569                 tails[j] = 0;
1570                 counts[i] = 0;
1571                 counts[j] = 0;
1572                 break;
1573               }
1574             }
1575             break;
1576           }
1577         }
1578       }
1579       //find beginning of a hole
1580       {
1581         for(int i = 0; i < i_size_less_1; ++i) {
1582           if(incoming[i] != 0) {
1583             for(int j = i+1; j < i_size; ++j) {
1584               if(incoming[j] != 0) {
1585                 //std::cout << "case6: " << i << " " << j << "\n";
1586                 //we are beginning a empty space
1587                 active_tail_arbitrary* holep = 0;
1588                 //if(c_size && counts[c_size_less_1] == 0 &&
1589                 //   counts_from_scanline[c_size_less_1].first.first.first.get(HORIZONTAL) == point.get(HORIZONTAL))
1590                 if(have_vertical_tail_from_below) {
1591                   holep = tails[c_size_less_1];
1592                   tails[c_size_less_1] = 0;
1593                   have_vertical_tail_from_below = false;
1594                 }
1595                 std::pair<active_tail_arbitrary*, active_tail_arbitrary*> tailPair =
1596                   active_tail_arbitrary::createActiveTailsAsPair(point, false, holep, fractureHoles_ != 0);
1597                 if(j == i_size_less_1 && incoming_count[j].first.get(HORIZONTAL) == point.get(HORIZONTAL)) {
1598                   //std::cout << "vertical element " << point << "\n";
1599                   returnValue = tailPair.first;
1600                   returnCount.first = point;
1601                   //returnCount = incoming_count[j];
1602                   returnCount.second = -1;
1603                 } else {
1604                   //std::cout << "new element " << j-1 << " " << incoming[j] << "\n";
1605                   //std::cout << point << " " <<  incoming_count[j].first << "\n";
1606                   elements.push_back(std::pair<vertex_half_edge,
1607                                      active_tail_arbitrary*>(vertex_half_edge(point,
1608                                                                               incoming_count[j].first, incoming[j]), tailPair.first));
1609                 }
1610                 //std::cout << "new element " << i-1 << " " << incoming[i] << "\n";
1611                 //std::cout << point << " " <<  incoming_count[i].first << "\n";
1612                 elements.push_back(std::pair<vertex_half_edge,
1613                                    active_tail_arbitrary*>(vertex_half_edge(point,
1614                                                                             incoming_count[i].first, incoming[i]), tailPair.second));
1615                 incoming[i] = 0;
1616                 incoming[j] = 0;
1617                 break;
1618               }
1619             }
1620             break;
1621           }
1622         }
1623       }
1624       if(have_vertical_tail_from_below) {
1625         if(tails.back()) {
1626           tails.back()->pushPoint(point);
1627           returnValue = tails.back();
1628           returnCount.first = point;
1629           returnCount.second = counts.back();
1630         }
1631       }
1632       //assert that tails, counts and incoming are all null
1633       return std::pair<std::pair<Point, int>, active_tail_arbitrary*>(returnCount, returnValue);
1634     }
1635 
print(const vertex_arbitrary_count & count)1636     static inline void print(const vertex_arbitrary_count& count) {
1637       for(unsigned i = 0; i < count.size(); ++i) {
1638         //std::cout << count[i].first.get(HORIZONTAL) << ",";
1639         //std::cout << count[i].first.get(VERTICAL) << ":";
1640         //std::cout << count[i].second << " ";
1641       } //std::cout << "\n";
1642     }
1643 
print(const scanline_data & data)1644     static inline void print(const scanline_data& data) {
1645       for(typename scanline_data::const_iterator itr = data.begin(); itr != data.end(); ++itr){
1646         //std::cout << itr->first.pt << ", " << itr->first.other_pt << "; ";
1647       } //std::cout << "\n";
1648     }
1649 
1650     template <class cT, class iT>
processEvent_(cT & output,iT inputBegin,iT inputEnd)1651     inline iT processEvent_(cT& output, iT inputBegin, iT inputEnd) {
1652       typedef typename high_precision_type<Unit>::type high_precision;
1653       //std::cout << "processEvent_\n";
1654       justBefore_ = true;
1655       //collect up all elements from the tree that are at the y
1656       //values of events in the input queue
1657       //create vector of new elements to add into tree
1658       active_tail_arbitrary* verticalTail = 0;
1659       std::pair<Point, int> verticalCount(Point(0, 0), 0);
1660       iT currentIter = inputBegin;
1661       std::vector<iterator> elementIters;
1662       std::vector<std::pair<vertex_half_edge, active_tail_arbitrary*> > elements;
1663       while(currentIter != inputEnd && currentIter->pt.get(HORIZONTAL) == x_) {
1664         //std::cout << "loop\n";
1665         Unit currentY = (*currentIter).pt.get(VERTICAL);
1666         //std::cout << "current Y " << currentY << "\n";
1667         //std::cout << "scanline size " << scanData_.size() << "\n";
1668         //print(scanData_);
1669         iterator iter = lookUp_(currentY);
1670         //std::cout << "found element in scanline " << (iter != scanData_.end()) << "\n";
1671         //int counts[4] = {0, 0, 0, 0};
1672         incoming_count counts_from_scanline;
1673         //std::cout << "finding elements in tree\n";
1674         //if(iter != scanData_.end())
1675         //  std::cout << "first iter y is " << iter->first.evalAtX(x_) << "\n";
1676         while(iter != scanData_.end() &&
1677               ((iter->first.pt.x() == x_ && iter->first.pt.y() == currentY) ||
1678                (iter->first.other_pt.x() == x_ && iter->first.other_pt.y() == currentY))) {
1679                 //iter->first.evalAtX(x_) == (high_precision)currentY) {
1680           //std::cout << "loop2\n";
1681           elementIters.push_back(iter);
1682           counts_from_scanline.push_back(std::pair<std::pair<std::pair<Point, Point>, int>, active_tail_arbitrary*>
1683                                          (std::pair<std::pair<Point, Point>, int>(std::pair<Point, Point>(iter->first.pt,
1684                                                                                                           iter->first.other_pt),
1685                                                                                   iter->first.count),
1686                                           iter->second));
1687           ++iter;
1688         }
1689         Point currentPoint(x_, currentY);
1690         //std::cout << "counts_from_scanline size " << counts_from_scanline.size() << "\n";
1691         sort_incoming_count(counts_from_scanline, currentPoint);
1692 
1693         vertex_arbitrary_count incoming;
1694         //std::cout << "aggregating\n";
1695         do {
1696           //std::cout << "loop3\n";
1697           const vertex_half_edge& elem = *currentIter;
1698           incoming.push_back(std::pair<Point, int>(elem.other_pt, elem.count));
1699           ++currentIter;
1700         } while(currentIter != inputEnd && currentIter->pt.get(VERTICAL) == currentY &&
1701                 currentIter->pt.get(HORIZONTAL) == x_);
1702         //print(incoming);
1703         sort_vertex_arbitrary_count(incoming, currentPoint);
1704         //std::cout << currentPoint.get(HORIZONTAL) << "," << currentPoint.get(VERTICAL) << "\n";
1705         //print(incoming);
1706         //std::cout << "incoming counts from input size " << incoming.size() << "\n";
1707         //compact_vertex_arbitrary_count(currentPoint, incoming);
1708         vertex_arbitrary_count tmp;
1709         tmp.reserve(incoming.size());
1710         for(std::size_t i = 0; i < incoming.size(); ++i) {
1711           if(currentPoint < incoming[i].first) {
1712             tmp.push_back(incoming[i]);
1713           }
1714         }
1715         incoming.swap(tmp);
1716         //std::cout << "incoming counts from input size " << incoming.size() << "\n";
1717         //now counts_from_scanline has the data from the left and
1718         //incoming has the data from the right at this point
1719         //cancel out any end points
1720         if(verticalTail) {
1721           //std::cout << "adding vertical tail to counts from scanline\n";
1722           //std::cout << -verticalCount.second << "\n";
1723           counts_from_scanline.push_back(std::pair<std::pair<std::pair<Point, Point>, int>, active_tail_arbitrary*>
1724                                          (std::pair<std::pair<Point, Point>, int>(std::pair<Point, Point>(verticalCount.first,
1725                                                                                                           currentPoint),
1726                                                                                   -verticalCount.second),
1727                                           verticalTail));
1728         }
1729         if(!incoming.empty() && incoming.back().first.get(HORIZONTAL) == x_) {
1730           //std::cout << "inverted vertical event\n";
1731           incoming.back().second *= -1;
1732         }
1733         //std::cout << "calling processPoint_\n";
1734         std::pair<std::pair<Point, int>, active_tail_arbitrary*> result = processPoint_(output, elements, Point(x_, currentY), counts_from_scanline, incoming);
1735         verticalCount = result.first;
1736         verticalTail = result.second;
1737         //if(verticalTail) {
1738         //  std::cout << "have vertical tail\n";
1739         //  std::cout << verticalCount.second << "\n";
1740         //}
1741         if(verticalTail && !(verticalCount.second)) {
1742           //we got a hole out of the point we just processed
1743           //iter is still at the next y element above the current y value in the tree
1744           //std::cout << "checking whether ot handle hole\n";
1745           if(currentIter == inputEnd ||
1746              currentIter->pt.get(HORIZONTAL) != x_ ||
1747              scanline_base<Unit>::on_above_or_below(currentIter->pt, half_edge(iter->first.pt, iter->first.other_pt)) != -1) {
1748             //(high_precision)(currentIter->pt.get(VERTICAL)) >= iter->first.evalAtX(x_)) {
1749 
1750             //std::cout << "handle hole here\n";
1751             if(fractureHoles_) {
1752               //std::cout << "fracture hole here\n";
1753               //we need to handle the hole now and not at the next input vertex
1754               active_tail_arbitrary* at = iter->second;
1755               high_precision precise_y = iter->first.evalAtX(x_);
1756               Unit fracture_y = convert_high_precision_type<Unit>(precise_y);
1757               if(precise_y < fracture_y) --fracture_y;
1758               Point point(x_, fracture_y);
1759               verticalTail->getOtherActiveTail()->pushPoint(point);
1760               iter->second = verticalTail->getOtherActiveTail();
1761               at->pushPoint(point);
1762               verticalTail->join(at);
1763               delete at;
1764               delete verticalTail;
1765               verticalTail = 0;
1766             } else {
1767               //std::cout << "push hole onto list\n";
1768               iter->second->addHole(verticalTail);
1769               verticalTail = 0;
1770             }
1771           }
1772         }
1773       }
1774       //std::cout << "erasing\n";
1775       //erase all elements from the tree
1776       for(typename std::vector<iterator>::iterator iter = elementIters.begin();
1777           iter != elementIters.end(); ++iter) {
1778         //std::cout << "erasing loop\n";
1779         scanData_.erase(*iter);
1780       }
1781       //switch comparison tie breaking policy
1782       justBefore_ = false;
1783       //add new elements into tree
1784       //std::cout << "inserting\n";
1785       for(typename std::vector<std::pair<vertex_half_edge, active_tail_arbitrary*> >::iterator iter = elements.begin();
1786           iter != elements.end(); ++iter) {
1787         //std::cout << "inserting loop\n";
1788         scanData_.insert(scanData_.end(), *iter);
1789       }
1790       //std::cout << "end processEvent\n";
1791       return currentIter;
1792     }
1793 
lookUp_(Unit y)1794     inline iterator lookUp_(Unit y){
1795       //if just before then we need to look from 1 not -1
1796       //std::cout << "just before " << justBefore_ << "\n";
1797       return scanData_.lower_bound(vertex_half_edge(Point(x_, y), Point(x_, y+1), 0));
1798     }
1799 
1800   public: //test functions
1801 
1802     template <typename stream_type>
testPolygonArbitraryFormationRect(stream_type & stdcout)1803     static inline bool testPolygonArbitraryFormationRect(stream_type& stdcout) {
1804       stdcout << "testing polygon formation\n";
1805       polygon_arbitrary_formation pf(true);
1806       std::vector<polygon_data<Unit> > polys;
1807       std::vector<vertex_half_edge> data;
1808       data.push_back(vertex_half_edge(Point(0, 0), Point(10, 0), 1));
1809       data.push_back(vertex_half_edge(Point(0, 0), Point(0, 10), 1));
1810       data.push_back(vertex_half_edge(Point(0, 10), Point(0, 0), -1));
1811       data.push_back(vertex_half_edge(Point(0, 10), Point(10, 10), -1));
1812       data.push_back(vertex_half_edge(Point(10, 0), Point(0, 0), -1));
1813       data.push_back(vertex_half_edge(Point(10, 0), Point(10, 10), -1));
1814       data.push_back(vertex_half_edge(Point(10, 10), Point(10, 0), 1));
1815       data.push_back(vertex_half_edge(Point(10, 10), Point(0, 10), 1));
1816       polygon_sort(data.begin(), data.end());
1817       pf.scan(polys, data.begin(), data.end());
1818       stdcout << "result size: " << polys.size() << "\n";
1819       for(std::size_t i = 0; i < polys.size(); ++i) {
1820         stdcout << polys[i] << "\n";
1821       }
1822       stdcout << "done testing polygon formation\n";
1823       return true;
1824     }
1825 
1826     template <typename stream_type>
testPolygonArbitraryFormationP1(stream_type & stdcout)1827     static inline bool testPolygonArbitraryFormationP1(stream_type& stdcout) {
1828       stdcout << "testing polygon formation P1\n";
1829       polygon_arbitrary_formation pf(true);
1830       std::vector<polygon_data<Unit> > polys;
1831       std::vector<vertex_half_edge> data;
1832       data.push_back(vertex_half_edge(Point(0, 0), Point(10, 10), 1));
1833       data.push_back(vertex_half_edge(Point(0, 0), Point(0, 10), 1));
1834       data.push_back(vertex_half_edge(Point(0, 10), Point(0, 0), -1));
1835       data.push_back(vertex_half_edge(Point(0, 10), Point(10, 20), -1));
1836       data.push_back(vertex_half_edge(Point(10, 10), Point(0, 0), -1));
1837       data.push_back(vertex_half_edge(Point(10, 10), Point(10, 20), -1));
1838       data.push_back(vertex_half_edge(Point(10, 20), Point(10, 10), 1));
1839       data.push_back(vertex_half_edge(Point(10, 20), Point(0, 10), 1));
1840       polygon_sort(data.begin(), data.end());
1841       pf.scan(polys, data.begin(), data.end());
1842       stdcout << "result size: " << polys.size() << "\n";
1843       for(std::size_t i = 0; i < polys.size(); ++i) {
1844         stdcout << polys[i] << "\n";
1845       }
1846       stdcout << "done testing polygon formation\n";
1847       return true;
1848     }
1849 
1850     template <typename stream_type>
testPolygonArbitraryFormationP2(stream_type & stdcout)1851     static inline bool testPolygonArbitraryFormationP2(stream_type& stdcout) {
1852       stdcout << "testing polygon formation P2\n";
1853       polygon_arbitrary_formation pf(true);
1854       std::vector<polygon_data<Unit> > polys;
1855       std::vector<vertex_half_edge> data;
1856       data.push_back(vertex_half_edge(Point(-3, 1), Point(2, -4), 1));
1857       data.push_back(vertex_half_edge(Point(-3, 1), Point(-2, 2), -1));
1858       data.push_back(vertex_half_edge(Point(-2, 2), Point(2, 4), -1));
1859       data.push_back(vertex_half_edge(Point(-2, 2), Point(-3, 1), 1));
1860       data.push_back(vertex_half_edge(Point(2, -4), Point(-3, 1), -1));
1861       data.push_back(vertex_half_edge(Point(2, -4), Point(2, 4), -1));
1862       data.push_back(vertex_half_edge(Point(2, 4), Point(-2, 2), 1));
1863       data.push_back(vertex_half_edge(Point(2, 4), Point(2, -4), 1));
1864       polygon_sort(data.begin(), data.end());
1865       pf.scan(polys, data.begin(), data.end());
1866       stdcout << "result size: " << polys.size() << "\n";
1867       for(std::size_t i = 0; i < polys.size(); ++i) {
1868         stdcout << polys[i] << "\n";
1869       }
1870       stdcout << "done testing polygon formation\n";
1871       return true;
1872     }
1873 
1874 
1875     template <typename stream_type>
testPolygonArbitraryFormationPolys(stream_type & stdcout)1876     static inline bool testPolygonArbitraryFormationPolys(stream_type& stdcout) {
1877       stdcout << "testing polygon formation polys\n";
1878       polygon_arbitrary_formation pf(false);
1879       std::vector<polygon_with_holes_data<Unit> > polys;
1880       polygon_arbitrary_formation pf2(true);
1881       std::vector<polygon_with_holes_data<Unit> > polys2;
1882       std::vector<vertex_half_edge> data;
1883       data.push_back(vertex_half_edge(Point(0, 0), Point(100, 1), 1));
1884       data.push_back(vertex_half_edge(Point(0, 0), Point(1, 100), -1));
1885       data.push_back(vertex_half_edge(Point(1, 100), Point(0, 0), 1));
1886       data.push_back(vertex_half_edge(Point(1, 100), Point(101, 101), -1));
1887       data.push_back(vertex_half_edge(Point(100, 1), Point(0, 0), -1));
1888       data.push_back(vertex_half_edge(Point(100, 1), Point(101, 101), 1));
1889       data.push_back(vertex_half_edge(Point(101, 101), Point(100, 1), -1));
1890       data.push_back(vertex_half_edge(Point(101, 101), Point(1, 100), 1));
1891 
1892       data.push_back(vertex_half_edge(Point(2, 2), Point(10, 2), -1));
1893       data.push_back(vertex_half_edge(Point(2, 2), Point(2, 10), -1));
1894       data.push_back(vertex_half_edge(Point(2, 10), Point(2, 2), 1));
1895       data.push_back(vertex_half_edge(Point(2, 10), Point(10, 10), 1));
1896       data.push_back(vertex_half_edge(Point(10, 2), Point(2, 2), 1));
1897       data.push_back(vertex_half_edge(Point(10, 2), Point(10, 10), 1));
1898       data.push_back(vertex_half_edge(Point(10, 10), Point(10, 2), -1));
1899       data.push_back(vertex_half_edge(Point(10, 10), Point(2, 10), -1));
1900 
1901       data.push_back(vertex_half_edge(Point(2, 12), Point(10, 12), -1));
1902       data.push_back(vertex_half_edge(Point(2, 12), Point(2, 22), -1));
1903       data.push_back(vertex_half_edge(Point(2, 22), Point(2, 12), 1));
1904       data.push_back(vertex_half_edge(Point(2, 22), Point(10, 22), 1));
1905       data.push_back(vertex_half_edge(Point(10, 12), Point(2, 12), 1));
1906       data.push_back(vertex_half_edge(Point(10, 12), Point(10, 22), 1));
1907       data.push_back(vertex_half_edge(Point(10, 22), Point(10, 12), -1));
1908       data.push_back(vertex_half_edge(Point(10, 22), Point(2, 22), -1));
1909 
1910       polygon_sort(data.begin(), data.end());
1911       pf.scan(polys, data.begin(), data.end());
1912       stdcout << "result size: " << polys.size() << "\n";
1913       for(std::size_t i = 0; i < polys.size(); ++i) {
1914         stdcout << polys[i] << "\n";
1915       }
1916       pf2.scan(polys2, data.begin(), data.end());
1917       stdcout << "result size: " << polys2.size() << "\n";
1918       for(std::size_t i = 0; i < polys2.size(); ++i) {
1919         stdcout << polys2[i] << "\n";
1920       }
1921       stdcout << "done testing polygon formation\n";
1922       return true;
1923     }
1924 
1925     template <typename stream_type>
testPolygonArbitraryFormationSelfTouch1(stream_type & stdcout)1926     static inline bool testPolygonArbitraryFormationSelfTouch1(stream_type& stdcout) {
1927       stdcout << "testing polygon formation self touch 1\n";
1928       polygon_arbitrary_formation pf(true);
1929       std::vector<polygon_data<Unit> > polys;
1930       std::vector<vertex_half_edge> data;
1931       data.push_back(vertex_half_edge(Point(0, 0), Point(10, 0), 1));
1932       data.push_back(vertex_half_edge(Point(0, 0), Point(0, 10), 1));
1933 
1934       data.push_back(vertex_half_edge(Point(0, 10), Point(0, 0), -1));
1935       data.push_back(vertex_half_edge(Point(0, 10), Point(5, 10), -1));
1936 
1937       data.push_back(vertex_half_edge(Point(10, 0), Point(0, 0), -1));
1938       data.push_back(vertex_half_edge(Point(10, 0), Point(10, 5), -1));
1939 
1940       data.push_back(vertex_half_edge(Point(10, 5), Point(10, 0), 1));
1941       data.push_back(vertex_half_edge(Point(10, 5), Point(5, 5), 1));
1942 
1943       data.push_back(vertex_half_edge(Point(5, 10), Point(5, 5), 1));
1944       data.push_back(vertex_half_edge(Point(5, 10), Point(0, 10), 1));
1945 
1946       data.push_back(vertex_half_edge(Point(5, 2), Point(5, 5), -1));
1947       data.push_back(vertex_half_edge(Point(5, 2), Point(7, 2), -1));
1948 
1949       data.push_back(vertex_half_edge(Point(5, 5), Point(5, 10), -1));
1950       data.push_back(vertex_half_edge(Point(5, 5), Point(5, 2), 1));
1951       data.push_back(vertex_half_edge(Point(5, 5), Point(10, 5), -1));
1952       data.push_back(vertex_half_edge(Point(5, 5), Point(7, 2), 1));
1953 
1954       data.push_back(vertex_half_edge(Point(7, 2), Point(5, 5), -1));
1955       data.push_back(vertex_half_edge(Point(7, 2), Point(5, 2), 1));
1956 
1957       polygon_sort(data.begin(), data.end());
1958       pf.scan(polys, data.begin(), data.end());
1959       stdcout << "result size: " << polys.size() << "\n";
1960       for(std::size_t i = 0; i < polys.size(); ++i) {
1961         stdcout << polys[i] << "\n";
1962       }
1963       stdcout << "done testing polygon formation\n";
1964       return true;
1965     }
1966 
1967     template <typename stream_type>
testPolygonArbitraryFormationSelfTouch2(stream_type & stdcout)1968     static inline bool testPolygonArbitraryFormationSelfTouch2(stream_type& stdcout) {
1969       stdcout << "testing polygon formation self touch 2\n";
1970       polygon_arbitrary_formation pf(true);
1971       std::vector<polygon_data<Unit> > polys;
1972       std::vector<vertex_half_edge> data;
1973       data.push_back(vertex_half_edge(Point(0, 0), Point(10, 0), 1));
1974       data.push_back(vertex_half_edge(Point(0, 0), Point(0, 10), 1));
1975 
1976       data.push_back(vertex_half_edge(Point(0, 10), Point(0, 0), -1));
1977       data.push_back(vertex_half_edge(Point(0, 10), Point(5, 10), -1));
1978 
1979       data.push_back(vertex_half_edge(Point(10, 0), Point(0, 0), -1));
1980       data.push_back(vertex_half_edge(Point(10, 0), Point(10, 5), -1));
1981 
1982       data.push_back(vertex_half_edge(Point(10, 5), Point(10, 0), 1));
1983       data.push_back(vertex_half_edge(Point(10, 5), Point(5, 5), 1));
1984 
1985       data.push_back(vertex_half_edge(Point(5, 10), Point(4, 1), -1));
1986       data.push_back(vertex_half_edge(Point(5, 10), Point(0, 10), 1));
1987 
1988       data.push_back(vertex_half_edge(Point(4, 1), Point(5, 10), 1));
1989       data.push_back(vertex_half_edge(Point(4, 1), Point(7, 2), -1));
1990 
1991       data.push_back(vertex_half_edge(Point(5, 5), Point(10, 5), -1));
1992       data.push_back(vertex_half_edge(Point(5, 5), Point(7, 2), 1));
1993 
1994       data.push_back(vertex_half_edge(Point(7, 2), Point(5, 5), -1));
1995       data.push_back(vertex_half_edge(Point(7, 2), Point(4, 1), 1));
1996 
1997       polygon_sort(data.begin(), data.end());
1998       pf.scan(polys, data.begin(), data.end());
1999       stdcout << "result size: " << polys.size() << "\n";
2000       for(std::size_t i = 0; i < polys.size(); ++i) {
2001         stdcout << polys[i] << "\n";
2002       }
2003       stdcout << "done testing polygon formation\n";
2004       return true;
2005     }
2006 
2007     template <typename stream_type>
testPolygonArbitraryFormationSelfTouch3(stream_type & stdcout)2008     static inline bool testPolygonArbitraryFormationSelfTouch3(stream_type& stdcout) {
2009       stdcout << "testing polygon formation self touch 3\n";
2010       polygon_arbitrary_formation pf(true);
2011       std::vector<polygon_data<Unit> > polys;
2012       std::vector<vertex_half_edge> data;
2013       data.push_back(vertex_half_edge(Point(0, 0), Point(10, 0), 1));
2014       data.push_back(vertex_half_edge(Point(0, 0), Point(0, 10), 1));
2015 
2016       data.push_back(vertex_half_edge(Point(0, 10), Point(0, 0), -1));
2017       data.push_back(vertex_half_edge(Point(0, 10), Point(6, 10), -1));
2018 
2019       data.push_back(vertex_half_edge(Point(10, 0), Point(0, 0), -1));
2020       data.push_back(vertex_half_edge(Point(10, 0), Point(10, 5), -1));
2021 
2022       data.push_back(vertex_half_edge(Point(10, 5), Point(10, 0), 1));
2023       data.push_back(vertex_half_edge(Point(10, 5), Point(5, 5), 1));
2024 
2025       data.push_back(vertex_half_edge(Point(6, 10), Point(4, 1), -1));
2026       data.push_back(vertex_half_edge(Point(6, 10), Point(0, 10), 1));
2027 
2028       data.push_back(vertex_half_edge(Point(4, 1), Point(6, 10), 1));
2029       data.push_back(vertex_half_edge(Point(4, 1), Point(7, 2), -1));
2030 
2031       data.push_back(vertex_half_edge(Point(5, 5), Point(10, 5), -1));
2032       data.push_back(vertex_half_edge(Point(5, 5), Point(7, 2), 1));
2033 
2034       data.push_back(vertex_half_edge(Point(7, 2), Point(5, 5), -1));
2035       data.push_back(vertex_half_edge(Point(7, 2), Point(4, 1), 1));
2036 
2037       polygon_sort(data.begin(), data.end());
2038       pf.scan(polys, data.begin(), data.end());
2039       stdcout << "result size: " << polys.size() << "\n";
2040       for(std::size_t i = 0; i < polys.size(); ++i) {
2041         stdcout << polys[i] << "\n";
2042       }
2043       stdcout << "done testing polygon formation\n";
2044       return true;
2045     }
2046 
2047     template <typename stream_type>
testPolygonArbitraryFormationColinear(stream_type & stdcout)2048     static inline bool testPolygonArbitraryFormationColinear(stream_type& stdcout) {
2049       stdcout << "testing polygon formation colinear 3\n";
2050       stdcout << "Polygon Set Data { <-3 2, -2 2>:1 <-3 2, -1 4>:-1 <-2 2, 0 2>:1 <-1 4, 0 2>:-1 } \n";
2051       polygon_arbitrary_formation pf(true);
2052       std::vector<polygon_data<Unit> > polys;
2053       std::vector<vertex_half_edge> data;
2054       data.push_back(vertex_half_edge(Point(-3, 2), Point(-2, 2), 1));
2055       data.push_back(vertex_half_edge(Point(-2, 2), Point(-3, 2), -1));
2056 
2057       data.push_back(vertex_half_edge(Point(-3, 2), Point(-1, 4), -1));
2058       data.push_back(vertex_half_edge(Point(-1, 4), Point(-3, 2), 1));
2059 
2060       data.push_back(vertex_half_edge(Point(-2, 2), Point(0, 2), 1));
2061       data.push_back(vertex_half_edge(Point(0, 2), Point(-2, 2), -1));
2062 
2063       data.push_back(vertex_half_edge(Point(-1, 4), Point(0, 2), -1));
2064       data.push_back(vertex_half_edge(Point(0, 2), Point(-1, 4), 1));
2065       polygon_sort(data.begin(), data.end());
2066       pf.scan(polys, data.begin(), data.end());
2067       stdcout << "result size: " << polys.size() << "\n";
2068       for(std::size_t i = 0; i < polys.size(); ++i) {
2069         stdcout << polys[i] << "\n";
2070       }
2071       stdcout << "done testing polygon formation\n";
2072       return true;
2073     }
2074 
2075     template <typename stream_type>
testSegmentIntersection(stream_type & stdcout)2076     static inline bool testSegmentIntersection(stream_type& stdcout) {
2077       stdcout << "testing segment intersection\n";
2078       half_edge he1, he2;
2079       he1.first = Point(0, 0);
2080       he1.second = Point(10, 10);
2081       he2.first = Point(0, 0);
2082       he2.second = Point(10, 20);
2083       Point result;
2084       bool b = scanline_base<Unit>::compute_intersection(result, he1, he2);
2085       if(!b || result != Point(0, 0)) return false;
2086       he1.first = Point(0, 10);
2087       b = scanline_base<Unit>::compute_intersection(result, he1, he2);
2088       if(!b || result != Point(5, 10)) return false;
2089       he1.first = Point(0, 11);
2090       b = scanline_base<Unit>::compute_intersection(result, he1, he2);
2091       if(!b || result != Point(5, 10)) return false;
2092       he1.first = Point(0, 0);
2093       he1.second = Point(1, 9);
2094       he2.first = Point(0, 9);
2095       he2.second = Point(1, 0);
2096       b = scanline_base<Unit>::compute_intersection(result, he1, he2);
2097       if(!b || result != Point(0, 4)) return false;
2098 
2099       he1.first = Point(0, -10);
2100       he1.second = Point(1, -1);
2101       he2.first = Point(0, -1);
2102       he2.second = Point(1, -10);
2103       b = scanline_base<Unit>::compute_intersection(result, he1, he2);
2104       if(!b || result != Point(0, -5)) return false;
2105       he1.first = Point((std::numeric_limits<int>::max)(), (std::numeric_limits<int>::max)()-1);
2106       he1.second = Point((std::numeric_limits<int>::min)(), (std::numeric_limits<int>::max)());
2107       //he1.second = Point(0, (std::numeric_limits<int>::max)());
2108       he2.first = Point((std::numeric_limits<int>::max)()-1, (std::numeric_limits<int>::max)());
2109       he2.second = Point((std::numeric_limits<int>::max)(), (std::numeric_limits<int>::min)());
2110       //he2.second = Point((std::numeric_limits<int>::max)(), 0);
2111       b = scanline_base<Unit>::compute_intersection(result, he1, he2);
2112       //b is false because of overflow error
2113       he1.first = Point(1000, 2000);
2114       he1.second = Point(1010, 2010);
2115       he2.first = Point(1000, 2000);
2116       he2.second = Point(1010, 2020);
2117       b = scanline_base<Unit>::compute_intersection(result, he1, he2);
2118       if(!b || result != Point(1000, 2000)) return false;
2119 
2120       return b;
2121     }
2122 
2123   };
2124 
2125   template <typename Unit>
2126   class poly_line_arbitrary_hole_data {
2127   private:
2128     typedef typename polygon_arbitrary_formation<Unit>::active_tail_arbitrary active_tail_arbitrary;
2129     active_tail_arbitrary* p_;
2130   public:
2131     typedef point_data<Unit> Point;
2132     typedef Point point_type;
2133     typedef Unit coordinate_type;
2134     typedef typename active_tail_arbitrary::iterator iterator_type;
2135     //typedef iterator_points_to_compact<iterator_type, Point> compact_iterator_type;
2136 
2137     typedef iterator_type iterator;
poly_line_arbitrary_hole_data()2138     inline poly_line_arbitrary_hole_data() : p_(0) {}
poly_line_arbitrary_hole_data(active_tail_arbitrary * p)2139     inline poly_line_arbitrary_hole_data(active_tail_arbitrary* p) : p_(p) {}
2140     //use default copy and assign
begin() const2141     inline iterator begin() const { return p_->getTail()->begin(); }
end() const2142     inline iterator end() const { return p_->getTail()->end(); }
size() const2143     inline std::size_t size() const { return 0; }
2144   };
2145 
2146   template <typename Unit>
2147   class poly_line_arbitrary_polygon_data {
2148   private:
2149     typedef typename polygon_arbitrary_formation<Unit>::active_tail_arbitrary active_tail_arbitrary;
2150     active_tail_arbitrary* p_;
2151   public:
2152     typedef point_data<Unit> Point;
2153     typedef Point point_type;
2154     typedef Unit coordinate_type;
2155     typedef typename active_tail_arbitrary::iterator iterator_type;
2156     //typedef iterator_points_to_compact<iterator_type, Point> compact_iterator_type;
2157     typedef typename coordinate_traits<Unit>::coordinate_distance area_type;
2158 
2159     class iterator_holes_type {
2160     private:
2161       typedef poly_line_arbitrary_hole_data<Unit> holeType;
2162       mutable holeType hole_;
2163       typename active_tail_arbitrary::iteratorHoles itr_;
2164 
2165     public:
2166       typedef std::forward_iterator_tag iterator_category;
2167       typedef holeType value_type;
2168       typedef std::ptrdiff_t difference_type;
2169       typedef const holeType* pointer; //immutable
2170       typedef const holeType& reference; //immutable
iterator_holes_type()2171       inline iterator_holes_type() : hole_(), itr_() {}
iterator_holes_type(typename active_tail_arbitrary::iteratorHoles itr)2172       inline iterator_holes_type(typename active_tail_arbitrary::iteratorHoles itr) : hole_(), itr_(itr) {}
iterator_holes_type(const iterator_holes_type & that)2173       inline iterator_holes_type(const iterator_holes_type& that) : hole_(that.hole_), itr_(that.itr_) {}
operator =(const iterator_holes_type & that)2174       inline iterator_holes_type& operator=(const iterator_holes_type& that) {
2175         itr_ = that.itr_;
2176         return *this;
2177       }
operator ==(const iterator_holes_type & that)2178       inline bool operator==(const iterator_holes_type& that) { return itr_ == that.itr_; }
operator !=(const iterator_holes_type & that)2179       inline bool operator!=(const iterator_holes_type& that) { return itr_ != that.itr_; }
operator ++()2180       inline iterator_holes_type& operator++() {
2181         ++itr_;
2182         return *this;
2183       }
operator ++(int)2184       inline const iterator_holes_type operator++(int) {
2185         iterator_holes_type tmp = *this;
2186         ++(*this);
2187         return tmp;
2188       }
operator *()2189       inline reference operator*() {
2190         hole_ = holeType(*itr_);
2191         return hole_;
2192       }
2193     };
2194 
2195     typedef poly_line_arbitrary_hole_data<Unit> hole_type;
2196 
poly_line_arbitrary_polygon_data()2197     inline poly_line_arbitrary_polygon_data() : p_(0) {}
poly_line_arbitrary_polygon_data(active_tail_arbitrary * p)2198     inline poly_line_arbitrary_polygon_data(active_tail_arbitrary* p) : p_(p) {}
2199     //use default copy and assign
begin() const2200     inline iterator_type begin() const { return p_->getTail()->begin(); }
end() const2201     inline iterator_type end() const { return p_->getTail()->end(); }
2202     //inline compact_iterator_type begin_compact() const { return p_->getTail()->begin(); }
2203     //inline compact_iterator_type end_compact() const { return p_->getTail()->end(); }
begin_holes() const2204     inline iterator_holes_type begin_holes() const { return iterator_holes_type(p_->getHoles().begin()); }
end_holes() const2205     inline iterator_holes_type end_holes() const { return iterator_holes_type(p_->getHoles().end()); }
yield()2206     inline active_tail_arbitrary* yield() { return p_; }
2207     //stub out these four required functions that will not be used but are needed for the interface
size_holes() const2208     inline std::size_t size_holes() const { return 0; }
size() const2209     inline std::size_t size() const { return 0; }
2210   };
2211 
2212   template <typename Unit>
2213   class trapezoid_arbitrary_formation : public polygon_arbitrary_formation<Unit> {
2214   private:
2215     typedef typename scanline_base<Unit>::Point Point;
2216     typedef typename scanline_base<Unit>::half_edge half_edge;
2217     typedef typename scanline_base<Unit>::vertex_half_edge vertex_half_edge;
2218     typedef typename scanline_base<Unit>::less_vertex_half_edge less_vertex_half_edge;
2219 
2220     typedef typename polygon_arbitrary_formation<Unit>::poly_line_arbitrary poly_line_arbitrary;
2221 
2222     typedef typename polygon_arbitrary_formation<Unit>::active_tail_arbitrary active_tail_arbitrary;
2223 
2224     typedef std::vector<std::pair<Point, int> > vertex_arbitrary_count;
2225 
2226     typedef typename polygon_arbitrary_formation<Unit>::less_half_edge_count less_half_edge_count;
2227 
2228     typedef std::vector<std::pair<std::pair<std::pair<Point, Point>, int>, active_tail_arbitrary*> > incoming_count;
2229 
2230     typedef typename polygon_arbitrary_formation<Unit>::less_incoming_count less_incoming_count;
2231 
2232     typedef typename polygon_arbitrary_formation<Unit>::vertex_arbitrary_compact vertex_arbitrary_compact;
2233 
2234   private:
2235     //definitions
2236     typedef std::map<vertex_half_edge, active_tail_arbitrary*, less_vertex_half_edge> scanline_data;
2237     typedef typename scanline_data::iterator iterator;
2238     typedef typename scanline_data::const_iterator const_iterator;
2239 
2240     //data
2241   public:
trapezoid_arbitrary_formation()2242     inline trapezoid_arbitrary_formation() : polygon_arbitrary_formation<Unit>() {}
trapezoid_arbitrary_formation(const trapezoid_arbitrary_formation & that)2243     inline trapezoid_arbitrary_formation(const trapezoid_arbitrary_formation& that) : polygon_arbitrary_formation<Unit>(that) {}
operator =(const trapezoid_arbitrary_formation & that)2244     inline trapezoid_arbitrary_formation& operator=(const trapezoid_arbitrary_formation& that) {
2245       * static_cast<polygon_arbitrary_formation<Unit>*>(this) = * static_cast<polygon_arbitrary_formation<Unit>*>(&that);
2246       return *this;
2247     }
2248 
2249     //cT is an output container of Polygon45 or Polygon45WithHoles
2250     //iT is an iterator over vertex_half_edge elements
2251     //inputBegin - inputEnd is a range of sorted iT that represents
2252     //one or more scanline stops worth of data
2253     template <class cT, class iT>
scan(cT & output,iT inputBegin,iT inputEnd)2254     void scan(cT& output, iT inputBegin, iT inputEnd) {
2255       //std::cout << "1\n";
2256       while(inputBegin != inputEnd) {
2257         //std::cout << "2\n";
2258         polygon_arbitrary_formation<Unit>::x_ = (*inputBegin).pt.get(HORIZONTAL);
2259         //std::cout << "SCAN FORMATION " << x_ << "\n";
2260         //std::cout << "x_ = " << x_ << "\n";
2261         //std::cout << "scan line size: " << scanData_.size() << "\n";
2262         inputBegin = processEvent_(output, inputBegin, inputEnd);
2263       }
2264       //std::cout << "scan line size: " << scanData_.size() << "\n";
2265     }
2266 
2267   private:
2268     //functions
getVerticalPair_(std::pair<active_tail_arbitrary *,active_tail_arbitrary * > & verticalPair,iterator previter)2269     inline void getVerticalPair_(std::pair<active_tail_arbitrary*,
2270                                  active_tail_arbitrary*>& verticalPair,
2271                                  iterator previter) {
2272       active_tail_arbitrary* iterTail = (*previter).second;
2273       Point prevPoint(polygon_arbitrary_formation<Unit>::x_,
2274                       convert_high_precision_type<Unit>(previter->first.evalAtX(polygon_arbitrary_formation<Unit>::x_)));
2275       iterTail->pushPoint(prevPoint);
2276       std::pair<active_tail_arbitrary*, active_tail_arbitrary*> tailPair =
2277         active_tail_arbitrary::createActiveTailsAsPair(prevPoint, true, 0, false);
2278       verticalPair.first = iterTail;
2279       verticalPair.second = tailPair.first;
2280       (*previter).second = tailPair.second;
2281     }
2282 
2283     template <class cT, class cT2>
2284     inline std::pair<std::pair<Point, int>, active_tail_arbitrary*>
processPoint_(cT & output,cT2 & elements,std::pair<active_tail_arbitrary *,active_tail_arbitrary * > & verticalPair,iterator previter,Point point,incoming_count & counts_from_scanline,vertex_arbitrary_count & incoming_count)2285     processPoint_(cT& output, cT2& elements,
2286                   std::pair<active_tail_arbitrary*, active_tail_arbitrary*>& verticalPair,
2287                   iterator previter, Point point, incoming_count& counts_from_scanline,
2288                   vertex_arbitrary_count& incoming_count) {
2289       //std::cout << "\nAT POINT: " <<  point << "\n";
2290       //join any closing solid corners
2291       std::vector<int> counts;
2292       std::vector<int> incoming;
2293       std::vector<active_tail_arbitrary*> tails;
2294       counts.reserve(counts_from_scanline.size());
2295       tails.reserve(counts_from_scanline.size());
2296       incoming.reserve(incoming_count.size());
2297       for(std::size_t i = 0; i < counts_from_scanline.size(); ++i) {
2298         counts.push_back(counts_from_scanline[i].first.second);
2299         tails.push_back(counts_from_scanline[i].second);
2300       }
2301       for(std::size_t i = 0; i < incoming_count.size(); ++i) {
2302         incoming.push_back(incoming_count[i].second);
2303         if(incoming_count[i].first < point) {
2304           incoming.back() = 0;
2305         }
2306       }
2307 
2308       active_tail_arbitrary* returnValue = 0;
2309       std::pair<active_tail_arbitrary*, active_tail_arbitrary*> verticalPairOut;
2310       verticalPairOut.first = 0;
2311       verticalPairOut.second = 0;
2312       std::pair<Point, int> returnCount(Point(0, 0), 0);
2313       int i_size_less_1 = (int)(incoming.size()) -1;
2314       int c_size_less_1 = (int)(counts.size()) -1;
2315       int i_size = incoming.size();
2316       int c_size = counts.size();
2317 
2318       bool have_vertical_tail_from_below = false;
2319       if(c_size &&
2320          scanline_base<Unit>::is_vertical(counts_from_scanline.back().first.first)) {
2321         have_vertical_tail_from_below = true;
2322       }
2323       //assert size = size_less_1 + 1
2324       //std::cout << tails.size() << " " << incoming.size() << " " << counts_from_scanline.size() << " " << incoming_count.size() << "\n";
2325       //         for(std::size_t i = 0; i < counts.size(); ++i) {
2326       //           std::cout << counts_from_scanline[i].first.first.first.get(HORIZONTAL) << ",";
2327       //           std::cout << counts_from_scanline[i].first.first.first.get(VERTICAL) << " ";
2328       //           std::cout << counts_from_scanline[i].first.first.second.get(HORIZONTAL) << ",";
2329       //           std::cout << counts_from_scanline[i].first.first.second.get(VERTICAL) << ":";
2330       //           std::cout << counts_from_scanline[i].first.second << " ";
2331       //         } std::cout << "\n";
2332       //         print(incoming_count);
2333       {
2334         for(int i = 0; i < c_size_less_1; ++i) {
2335           //std::cout << i << "\n";
2336           if(counts[i] == -1) {
2337             //std::cout << "fixed i\n";
2338             for(int j = i + 1; j < c_size; ++j) {
2339               //std::cout << j << "\n";
2340               if(counts[j]) {
2341                 if(counts[j] == 1) {
2342                   //std::cout << "case1: " << i << " " << j << "\n";
2343                   //if a figure is closed it will be written out by this function to output
2344                   active_tail_arbitrary::joinChains(point, tails[i], tails[j], true, output);
2345                   counts[i] = 0;
2346                   counts[j] = 0;
2347                   tails[i] = 0;
2348                   tails[j] = 0;
2349                 }
2350                 break;
2351               }
2352             }
2353           }
2354         }
2355       }
2356       //find any pairs of incoming edges that need to create pair for leading solid
2357       //std::cout << "checking case2\n";
2358       {
2359         for(int i = 0; i < i_size_less_1; ++i) {
2360           //std::cout << i << "\n";
2361           if(incoming[i] == 1) {
2362             //std::cout << "fixed i\n";
2363             for(int j = i + 1; j < i_size; ++j) {
2364               //std::cout << j << "\n";
2365               if(incoming[j]) {
2366                 //std::cout << incoming[j] << "\n";
2367                 if(incoming[j] == -1) {
2368                   //std::cout << "case2: " << i << " " << j << "\n";
2369                   //std::cout << "creating active tail pair\n";
2370                   std::pair<active_tail_arbitrary*, active_tail_arbitrary*> tailPair =
2371                     active_tail_arbitrary::createActiveTailsAsPair(point, true, 0, polygon_arbitrary_formation<Unit>::fractureHoles_ != 0);
2372                   //tailPair.first->print();
2373                   //tailPair.second->print();
2374                   if(j == i_size_less_1 && incoming_count[j].first.get(HORIZONTAL) == point.get(HORIZONTAL)) {
2375                     //vertical active tail becomes return value
2376                     returnValue = tailPair.first;
2377                     returnCount.first = point;
2378                     returnCount.second = 1;
2379                   } else {
2380                     //std::cout << "new element " << j-1 << " " << -1 << "\n";
2381                     //std::cout << point << " " <<  incoming_count[j].first << "\n";
2382                     elements.push_back(std::pair<vertex_half_edge,
2383                                        active_tail_arbitrary*>(vertex_half_edge(point,
2384                                                                                 incoming_count[j].first, -1), tailPair.first));
2385                   }
2386                   //std::cout << "new element " << i-1 << " " << 1 << "\n";
2387                   //std::cout << point << " " <<  incoming_count[i].first << "\n";
2388                   elements.push_back(std::pair<vertex_half_edge,
2389                                      active_tail_arbitrary*>(vertex_half_edge(point,
2390                                                                               incoming_count[i].first, 1), tailPair.second));
2391                   incoming[i] = 0;
2392                   incoming[j] = 0;
2393                 }
2394                 break;
2395               }
2396             }
2397           }
2398         }
2399       }
2400       //find any active tail that needs to pass through to an incoming edge
2401       //we expect to find no more than two pass through
2402 
2403       //find pass through with solid on top
2404       {
2405         //std::cout << "checking case 3\n";
2406         for(int i = 0; i < c_size; ++i) {
2407           //std::cout << i << "\n";
2408           if(counts[i] != 0) {
2409             if(counts[i] == 1) {
2410               //std::cout << "fixed i\n";
2411               for(int j = i_size_less_1; j >= 0; --j) {
2412                 if(incoming[j] != 0) {
2413                   if(incoming[j] == 1) {
2414                     //std::cout << "case3: " << i << " " << j << "\n";
2415                     //tails[i]->print();
2416                     //pass through solid on top
2417                     tails[i]->pushPoint(point);
2418                     //std::cout << "after push\n";
2419                     if(j == i_size_less_1 && incoming_count[j].first.get(HORIZONTAL) == point.get(HORIZONTAL)) {
2420                       returnValue = tails[i];
2421                       returnCount.first = point;
2422                       returnCount.second = -1;
2423                     } else {
2424                       std::pair<active_tail_arbitrary*, active_tail_arbitrary*> tailPair =
2425                         active_tail_arbitrary::createActiveTailsAsPair(point, true, 0, false);
2426                       verticalPairOut.first = tails[i];
2427                       verticalPairOut.second = tailPair.first;
2428                       elements.push_back(std::pair<vertex_half_edge,
2429                                          active_tail_arbitrary*>(vertex_half_edge(point,
2430                                                                                   incoming_count[j].first, incoming[j]), tailPair.second));
2431                     }
2432                     tails[i] = 0;
2433                     counts[i] = 0;
2434                     incoming[j] = 0;
2435                   }
2436                   break;
2437                 }
2438               }
2439             }
2440             break;
2441           }
2442         }
2443       }
2444       //std::cout << "checking case 4\n";
2445       //find pass through with solid on bottom
2446       {
2447         for(int i = c_size_less_1; i >= 0; --i) {
2448           //std::cout << "i = " << i << " with count " << counts[i] << "\n";
2449           if(counts[i] != 0) {
2450             if(counts[i] == -1) {
2451               for(int j = 0; j < i_size; ++j) {
2452                 if(incoming[j] != 0) {
2453                   if(incoming[j] == -1) {
2454                     //std::cout << "case4: " << i << " " << j << "\n";
2455                     //pass through solid on bottom
2456 
2457                     //if count from scanline is vertical
2458                     if(i == c_size_less_1 &&
2459                        counts_from_scanline[i].first.first.first.get(HORIZONTAL) ==
2460                        point.get(HORIZONTAL)) {
2461                        //if incoming count is vertical
2462                        if(j == i_size_less_1 &&
2463                           incoming_count[j].first.get(HORIZONTAL) == point.get(HORIZONTAL)) {
2464                          returnValue = tails[i];
2465                          returnCount.first = point;
2466                          returnCount.second = 1;
2467                        } else {
2468                          tails[i]->pushPoint(point);
2469                          elements.push_back(std::pair<vertex_half_edge,
2470                                          active_tail_arbitrary*>(vertex_half_edge(point,
2471                                                                                   incoming_count[j].first, incoming[j]), tails[i]));
2472                        }
2473                     } else if(j == i_size_less_1 &&
2474                               incoming_count[j].first.get(HORIZONTAL) ==
2475                               point.get(HORIZONTAL)) {
2476                       if(verticalPair.first == 0) {
2477                         getVerticalPair_(verticalPair, previter);
2478                       }
2479                       active_tail_arbitrary::joinChains(point, tails[i], verticalPair.first, true, output);
2480                       returnValue = verticalPair.second;
2481                       returnCount.first = point;
2482                       returnCount.second = 1;
2483                     } else {
2484                       //neither is vertical
2485                       if(verticalPair.first == 0) {
2486                         getVerticalPair_(verticalPair, previter);
2487                       }
2488                       active_tail_arbitrary::joinChains(point, tails[i], verticalPair.first, true, output);
2489                       verticalPair.second->pushPoint(point);
2490                       elements.push_back(std::pair<vertex_half_edge,
2491                                          active_tail_arbitrary*>(vertex_half_edge(point,
2492                                                                                   incoming_count[j].first, incoming[j]), verticalPair.second));
2493                     }
2494                     tails[i] = 0;
2495                     counts[i] = 0;
2496                     incoming[j] = 0;
2497                   }
2498                   break;
2499                 }
2500               }
2501             }
2502             break;
2503           }
2504         }
2505       }
2506       //find the end of a hole or the beginning of a hole
2507 
2508       //find end of a hole
2509       {
2510         for(int i = 0; i < c_size_less_1; ++i) {
2511           if(counts[i] != 0) {
2512             for(int j = i+1; j < c_size; ++j) {
2513               if(counts[j] != 0) {
2514                 //std::cout << "case5: " << i << " " << j << "\n";
2515                 //we are ending a hole and may potentially close a figure and have to handle the hole
2516                 tails[i]->pushPoint(point);
2517                 verticalPairOut.first = tails[i];
2518                 if(j == c_size_less_1 &&
2519                    counts_from_scanline[j].first.first.first.get(HORIZONTAL) ==
2520                    point.get(HORIZONTAL)) {
2521                   verticalPairOut.second = tails[j];
2522                 } else {
2523                   //need to close a trapezoid below
2524                   if(verticalPair.first == 0) {
2525                     getVerticalPair_(verticalPair, previter);
2526                   }
2527                   active_tail_arbitrary::joinChains(point, tails[j], verticalPair.first, true, output);
2528                   verticalPairOut.second = verticalPair.second;
2529                 }
2530                 tails[i] = 0;
2531                 tails[j] = 0;
2532                 counts[i] = 0;
2533                 counts[j] = 0;
2534                 break;
2535               }
2536             }
2537             break;
2538           }
2539         }
2540       }
2541       //find beginning of a hole
2542       {
2543         for(int i = 0; i < i_size_less_1; ++i) {
2544           if(incoming[i] != 0) {
2545             for(int j = i+1; j < i_size; ++j) {
2546               if(incoming[j] != 0) {
2547                 //std::cout << "case6: " << i << " " << j << "\n";
2548                 //we are beginning a empty space
2549                 if(verticalPair.first == 0) {
2550                   getVerticalPair_(verticalPair, previter);
2551                 }
2552                 verticalPair.second->pushPoint(point);
2553                 if(j == i_size_less_1 &&
2554                    incoming_count[j].first.get(HORIZONTAL) == point.get(HORIZONTAL)) {
2555                   returnValue = verticalPair.first;
2556                   returnCount.first = point;
2557                   returnCount.second = -1;
2558                 } else {
2559                   std::pair<active_tail_arbitrary*, active_tail_arbitrary*> tailPair =
2560                   active_tail_arbitrary::createActiveTailsAsPair(point, false, 0, false);
2561                   elements.push_back(std::pair<vertex_half_edge,
2562                                      active_tail_arbitrary*>(vertex_half_edge(point,
2563                                                                               incoming_count[j].first, incoming[j]), tailPair.second));
2564                   verticalPairOut.second = tailPair.first;
2565                   verticalPairOut.first = verticalPair.first;
2566                 }
2567                 elements.push_back(std::pair<vertex_half_edge,
2568                                    active_tail_arbitrary*>(vertex_half_edge(point,
2569                                                                             incoming_count[i].first, incoming[i]), verticalPair.second));
2570                 incoming[i] = 0;
2571                 incoming[j] = 0;
2572                 break;
2573               }
2574             }
2575             break;
2576           }
2577         }
2578       }
2579       if(have_vertical_tail_from_below) {
2580         if(tails.back()) {
2581           tails.back()->pushPoint(point);
2582           returnValue = tails.back();
2583           returnCount.first = point;
2584           returnCount.second = counts.back();
2585         }
2586       }
2587       verticalPair = verticalPairOut;
2588       //assert that tails, counts and incoming are all null
2589       return std::pair<std::pair<Point, int>, active_tail_arbitrary*>(returnCount, returnValue);
2590     }
2591 
print(const vertex_arbitrary_count & count)2592     static inline void print(const vertex_arbitrary_count& count) {
2593       for(unsigned i = 0; i < count.size(); ++i) {
2594         //std::cout << count[i].first.get(HORIZONTAL) << ",";
2595         //std::cout << count[i].first.get(VERTICAL) << ":";
2596         //std::cout << count[i].second << " ";
2597       } //std::cout << "\n";
2598     }
2599 
print(const scanline_data & data)2600     static inline void print(const scanline_data& data) {
2601       for(typename scanline_data::const_iterator itr = data.begin(); itr != data.end(); ++itr){
2602         //std::cout << itr->first.pt << ", " << itr->first.other_pt << "; ";
2603       } //std::cout << "\n";
2604     }
2605 
2606     template <class cT, class iT>
processEvent_(cT & output,iT inputBegin,iT inputEnd)2607     inline iT processEvent_(cT& output, iT inputBegin, iT inputEnd) {
2608       //typedef typename high_precision_type<Unit>::type high_precision;
2609       //std::cout << "processEvent_\n";
2610       polygon_arbitrary_formation<Unit>::justBefore_ = true;
2611       //collect up all elements from the tree that are at the y
2612       //values of events in the input queue
2613       //create vector of new elements to add into tree
2614       active_tail_arbitrary* verticalTail = 0;
2615       std::pair<active_tail_arbitrary*, active_tail_arbitrary*> verticalPair;
2616       std::pair<Point, int> verticalCount(Point(0, 0), 0);
2617       iT currentIter = inputBegin;
2618       std::vector<iterator> elementIters;
2619       std::vector<std::pair<vertex_half_edge, active_tail_arbitrary*> > elements;
2620       while(currentIter != inputEnd && currentIter->pt.get(HORIZONTAL) == polygon_arbitrary_formation<Unit>::x_) {
2621         //std::cout << "loop\n";
2622         Unit currentY = (*currentIter).pt.get(VERTICAL);
2623         //std::cout << "current Y " << currentY << "\n";
2624         //std::cout << "scanline size " << scanData_.size() << "\n";
2625         //print(scanData_);
2626         iterator iter = this->lookUp_(currentY);
2627         //std::cout << "found element in scanline " << (iter != scanData_.end()) << "\n";
2628         //int counts[4] = {0, 0, 0, 0};
2629         incoming_count counts_from_scanline;
2630         //std::cout << "finding elements in tree\n";
2631         //if(iter != scanData_.end())
2632         //  std::cout << "first iter y is " << iter->first.evalAtX(x_) << "\n";
2633         iterator previter = iter;
2634         if(previter != polygon_arbitrary_formation<Unit>::scanData_.end() &&
2635              previter->first.evalAtX(polygon_arbitrary_formation<Unit>::x_) >= currentY &&
2636              previter != polygon_arbitrary_formation<Unit>::scanData_.begin())
2637            --previter;
2638         while(iter != polygon_arbitrary_formation<Unit>::scanData_.end() &&
2639               ((iter->first.pt.x() == polygon_arbitrary_formation<Unit>::x_ && iter->first.pt.y() == currentY) ||
2640                (iter->first.other_pt.x() == polygon_arbitrary_formation<Unit>::x_ && iter->first.other_pt.y() == currentY))) {
2641                //iter->first.evalAtX(polygon_arbitrary_formation<Unit>::x_) == (high_precision)currentY) {
2642           //std::cout << "loop2\n";
2643           elementIters.push_back(iter);
2644           counts_from_scanline.push_back(std::pair<std::pair<std::pair<Point, Point>, int>, active_tail_arbitrary*>
2645                                          (std::pair<std::pair<Point, Point>, int>(std::pair<Point, Point>(iter->first.pt,
2646                                                                                                           iter->first.other_pt),
2647                                                                                   iter->first.count),
2648                                           iter->second));
2649           ++iter;
2650         }
2651         Point currentPoint(polygon_arbitrary_formation<Unit>::x_, currentY);
2652         //std::cout << "counts_from_scanline size " << counts_from_scanline.size() << "\n";
2653         this->sort_incoming_count(counts_from_scanline, currentPoint);
2654 
2655         vertex_arbitrary_count incoming;
2656         //std::cout << "aggregating\n";
2657         do {
2658           //std::cout << "loop3\n";
2659           const vertex_half_edge& elem = *currentIter;
2660           incoming.push_back(std::pair<Point, int>(elem.other_pt, elem.count));
2661           ++currentIter;
2662         } while(currentIter != inputEnd && currentIter->pt.get(VERTICAL) == currentY &&
2663                 currentIter->pt.get(HORIZONTAL) == polygon_arbitrary_formation<Unit>::x_);
2664         //print(incoming);
2665         this->sort_vertex_arbitrary_count(incoming, currentPoint);
2666         //std::cout << currentPoint.get(HORIZONTAL) << "," << currentPoint.get(VERTICAL) << "\n";
2667         //print(incoming);
2668         //std::cout << "incoming counts from input size " << incoming.size() << "\n";
2669         //compact_vertex_arbitrary_count(currentPoint, incoming);
2670         vertex_arbitrary_count tmp;
2671         tmp.reserve(incoming.size());
2672         for(std::size_t i = 0; i < incoming.size(); ++i) {
2673           if(currentPoint < incoming[i].first) {
2674             tmp.push_back(incoming[i]);
2675           }
2676         }
2677         incoming.swap(tmp);
2678         //std::cout << "incoming counts from input size " << incoming.size() << "\n";
2679         //now counts_from_scanline has the data from the left and
2680         //incoming has the data from the right at this point
2681         //cancel out any end points
2682         if(verticalTail) {
2683           //std::cout << "adding vertical tail to counts from scanline\n";
2684           //std::cout << -verticalCount.second << "\n";
2685           counts_from_scanline.push_back(std::pair<std::pair<std::pair<Point, Point>, int>, active_tail_arbitrary*>
2686                                          (std::pair<std::pair<Point, Point>, int>(std::pair<Point, Point>(verticalCount.first,
2687                                                                                                           currentPoint),
2688                                                                                   -verticalCount.second),
2689                                           verticalTail));
2690         }
2691         if(!incoming.empty() && incoming.back().first.get(HORIZONTAL) == polygon_arbitrary_formation<Unit>::x_) {
2692           //std::cout << "inverted vertical event\n";
2693           incoming.back().second *= -1;
2694         }
2695         //std::cout << "calling processPoint_\n";
2696            std::pair<std::pair<Point, int>, active_tail_arbitrary*> result = processPoint_(output, elements, verticalPair, previter, Point(polygon_arbitrary_formation<Unit>::x_, currentY), counts_from_scanline, incoming);
2697         verticalCount = result.first;
2698         verticalTail = result.second;
2699         if(verticalPair.first != 0 && iter != polygon_arbitrary_formation<Unit>::scanData_.end() &&
2700            (currentIter == inputEnd || currentIter->pt.x() != polygon_arbitrary_formation<Unit>::x_ ||
2701             currentIter->pt.y() > (*iter).first.evalAtX(polygon_arbitrary_formation<Unit>::x_))) {
2702           //splice vertical pair into edge above
2703           active_tail_arbitrary* tailabove = (*iter).second;
2704           Point point(polygon_arbitrary_formation<Unit>::x_,
2705                       convert_high_precision_type<Unit>((*iter).first.evalAtX(polygon_arbitrary_formation<Unit>::x_)));
2706           verticalPair.second->pushPoint(point);
2707           active_tail_arbitrary::joinChains(point, tailabove, verticalPair.first, true, output);
2708           (*iter).second = verticalPair.second;
2709           verticalPair.first = 0;
2710           verticalPair.second = 0;
2711         }
2712       }
2713       //std::cout << "erasing\n";
2714       //erase all elements from the tree
2715       for(typename std::vector<iterator>::iterator iter = elementIters.begin();
2716           iter != elementIters.end(); ++iter) {
2717         //std::cout << "erasing loop\n";
2718         polygon_arbitrary_formation<Unit>::scanData_.erase(*iter);
2719       }
2720       //switch comparison tie breaking policy
2721       polygon_arbitrary_formation<Unit>::justBefore_ = false;
2722       //add new elements into tree
2723       //std::cout << "inserting\n";
2724       for(typename std::vector<std::pair<vertex_half_edge, active_tail_arbitrary*> >::iterator iter = elements.begin();
2725           iter != elements.end(); ++iter) {
2726         //std::cout << "inserting loop\n";
2727         polygon_arbitrary_formation<Unit>::scanData_.insert(polygon_arbitrary_formation<Unit>::scanData_.end(), *iter);
2728       }
2729       //std::cout << "end processEvent\n";
2730       return currentIter;
2731     }
2732   public:
2733     template <typename stream_type>
testTrapezoidArbitraryFormationRect(stream_type & stdcout)2734     static inline bool testTrapezoidArbitraryFormationRect(stream_type& stdcout) {
2735       stdcout << "testing trapezoid formation\n";
2736       trapezoid_arbitrary_formation pf;
2737       std::vector<polygon_data<Unit> > polys;
2738       std::vector<vertex_half_edge> data;
2739       data.push_back(vertex_half_edge(Point(0, 0), Point(10, 0), 1));
2740       data.push_back(vertex_half_edge(Point(0, 0), Point(0, 10), 1));
2741       data.push_back(vertex_half_edge(Point(0, 10), Point(0, 0), -1));
2742       data.push_back(vertex_half_edge(Point(0, 10), Point(10, 10), -1));
2743       data.push_back(vertex_half_edge(Point(10, 0), Point(0, 0), -1));
2744       data.push_back(vertex_half_edge(Point(10, 0), Point(10, 10), -1));
2745       data.push_back(vertex_half_edge(Point(10, 10), Point(10, 0), 1));
2746       data.push_back(vertex_half_edge(Point(10, 10), Point(0, 10), 1));
2747       polygon_sort(data.begin(), data.end());
2748       pf.scan(polys, data.begin(), data.end());
2749       stdcout << "result size: " << polys.size() << "\n";
2750       for(std::size_t i = 0; i < polys.size(); ++i) {
2751         stdcout << polys[i] << "\n";
2752       }
2753       stdcout << "done testing trapezoid formation\n";
2754       return true;
2755     }
2756     template <typename stream_type>
testTrapezoidArbitraryFormationP1(stream_type & stdcout)2757     static inline bool testTrapezoidArbitraryFormationP1(stream_type& stdcout) {
2758       stdcout << "testing trapezoid formation P1\n";
2759       trapezoid_arbitrary_formation pf;
2760       std::vector<polygon_data<Unit> > polys;
2761       std::vector<vertex_half_edge> data;
2762       data.push_back(vertex_half_edge(Point(0, 0), Point(10, 10), 1));
2763       data.push_back(vertex_half_edge(Point(0, 0), Point(0, 10), 1));
2764       data.push_back(vertex_half_edge(Point(0, 10), Point(0, 0), -1));
2765       data.push_back(vertex_half_edge(Point(0, 10), Point(10, 20), -1));
2766       data.push_back(vertex_half_edge(Point(10, 10), Point(0, 0), -1));
2767       data.push_back(vertex_half_edge(Point(10, 10), Point(10, 20), -1));
2768       data.push_back(vertex_half_edge(Point(10, 20), Point(10, 10), 1));
2769       data.push_back(vertex_half_edge(Point(10, 20), Point(0, 10), 1));
2770       polygon_sort(data.begin(), data.end());
2771       pf.scan(polys, data.begin(), data.end());
2772       stdcout << "result size: " << polys.size() << "\n";
2773       for(std::size_t i = 0; i < polys.size(); ++i) {
2774         stdcout << polys[i] << "\n";
2775       }
2776       stdcout << "done testing trapezoid formation\n";
2777       return true;
2778     }
2779     template <typename stream_type>
testTrapezoidArbitraryFormationP2(stream_type & stdcout)2780     static inline bool testTrapezoidArbitraryFormationP2(stream_type& stdcout) {
2781       stdcout << "testing trapezoid formation P2\n";
2782       trapezoid_arbitrary_formation pf;
2783       std::vector<polygon_data<Unit> > polys;
2784       std::vector<vertex_half_edge> data;
2785       data.push_back(vertex_half_edge(Point(-3, 1), Point(2, -4), 1));
2786       data.push_back(vertex_half_edge(Point(-3, 1), Point(-2, 2), -1));
2787       data.push_back(vertex_half_edge(Point(-2, 2), Point(2, 4), -1));
2788       data.push_back(vertex_half_edge(Point(-2, 2), Point(-3, 1), 1));
2789       data.push_back(vertex_half_edge(Point(2, -4), Point(-3, 1), -1));
2790       data.push_back(vertex_half_edge(Point(2, -4), Point(2, 4), -1));
2791       data.push_back(vertex_half_edge(Point(2, 4), Point(-2, 2), 1));
2792       data.push_back(vertex_half_edge(Point(2, 4), Point(2, -4), 1));
2793       polygon_sort(data.begin(), data.end());
2794       pf.scan(polys, data.begin(), data.end());
2795       stdcout << "result size: " << polys.size() << "\n";
2796       for(std::size_t i = 0; i < polys.size(); ++i) {
2797         stdcout << polys[i] << "\n";
2798       }
2799       stdcout << "done testing trapezoid formation\n";
2800       return true;
2801     }
2802 
2803     template <typename stream_type>
testTrapezoidArbitraryFormationPolys(stream_type & stdcout)2804     static inline bool testTrapezoidArbitraryFormationPolys(stream_type& stdcout) {
2805       stdcout << "testing trapezoid formation polys\n";
2806       trapezoid_arbitrary_formation pf;
2807       std::vector<polygon_with_holes_data<Unit> > polys;
2808       //trapezoid_arbitrary_formation pf2(true);
2809       //std::vector<polygon_with_holes_data<Unit> > polys2;
2810       std::vector<vertex_half_edge> data;
2811       data.push_back(vertex_half_edge(Point(0, 0), Point(100, 1), 1));
2812       data.push_back(vertex_half_edge(Point(0, 0), Point(1, 100), -1));
2813       data.push_back(vertex_half_edge(Point(1, 100), Point(0, 0), 1));
2814       data.push_back(vertex_half_edge(Point(1, 100), Point(101, 101), -1));
2815       data.push_back(vertex_half_edge(Point(100, 1), Point(0, 0), -1));
2816       data.push_back(vertex_half_edge(Point(100, 1), Point(101, 101), 1));
2817       data.push_back(vertex_half_edge(Point(101, 101), Point(100, 1), -1));
2818       data.push_back(vertex_half_edge(Point(101, 101), Point(1, 100), 1));
2819 
2820       data.push_back(vertex_half_edge(Point(2, 2), Point(10, 2), -1));
2821       data.push_back(vertex_half_edge(Point(2, 2), Point(2, 10), -1));
2822       data.push_back(vertex_half_edge(Point(2, 10), Point(2, 2), 1));
2823       data.push_back(vertex_half_edge(Point(2, 10), Point(10, 10), 1));
2824       data.push_back(vertex_half_edge(Point(10, 2), Point(2, 2), 1));
2825       data.push_back(vertex_half_edge(Point(10, 2), Point(10, 10), 1));
2826       data.push_back(vertex_half_edge(Point(10, 10), Point(10, 2), -1));
2827       data.push_back(vertex_half_edge(Point(10, 10), Point(2, 10), -1));
2828 
2829       data.push_back(vertex_half_edge(Point(2, 12), Point(10, 12), -1));
2830       data.push_back(vertex_half_edge(Point(2, 12), Point(2, 22), -1));
2831       data.push_back(vertex_half_edge(Point(2, 22), Point(2, 12), 1));
2832       data.push_back(vertex_half_edge(Point(2, 22), Point(10, 22), 1));
2833       data.push_back(vertex_half_edge(Point(10, 12), Point(2, 12), 1));
2834       data.push_back(vertex_half_edge(Point(10, 12), Point(10, 22), 1));
2835       data.push_back(vertex_half_edge(Point(10, 22), Point(10, 12), -1));
2836       data.push_back(vertex_half_edge(Point(10, 22), Point(2, 22), -1));
2837 
2838       polygon_sort(data.begin(), data.end());
2839       pf.scan(polys, data.begin(), data.end());
2840       stdcout << "result size: " << polys.size() << "\n";
2841       for(std::size_t i = 0; i < polys.size(); ++i) {
2842         stdcout << polys[i] << "\n";
2843       }
2844       //pf2.scan(polys2, data.begin(), data.end());
2845       //stdcout << "result size: " << polys2.size() << "\n";
2846       //for(std::size_t i = 0; i < polys2.size(); ++i) {
2847       //  stdcout << polys2[i] << "\n";
2848       //}
2849       stdcout << "done testing trapezoid formation\n";
2850       return true;
2851     }
2852 
2853     template <typename stream_type>
testTrapezoidArbitraryFormationSelfTouch1(stream_type & stdcout)2854     static inline bool testTrapezoidArbitraryFormationSelfTouch1(stream_type& stdcout) {
2855       stdcout << "testing trapezoid formation self touch 1\n";
2856       trapezoid_arbitrary_formation pf;
2857       std::vector<polygon_data<Unit> > polys;
2858       std::vector<vertex_half_edge> data;
2859       data.push_back(vertex_half_edge(Point(0, 0), Point(10, 0), 1));
2860       data.push_back(vertex_half_edge(Point(0, 0), Point(0, 10), 1));
2861 
2862       data.push_back(vertex_half_edge(Point(0, 10), Point(0, 0), -1));
2863       data.push_back(vertex_half_edge(Point(0, 10), Point(5, 10), -1));
2864 
2865       data.push_back(vertex_half_edge(Point(10, 0), Point(0, 0), -1));
2866       data.push_back(vertex_half_edge(Point(10, 0), Point(10, 5), -1));
2867 
2868       data.push_back(vertex_half_edge(Point(10, 5), Point(10, 0), 1));
2869       data.push_back(vertex_half_edge(Point(10, 5), Point(5, 5), 1));
2870 
2871       data.push_back(vertex_half_edge(Point(5, 10), Point(5, 5), 1));
2872       data.push_back(vertex_half_edge(Point(5, 10), Point(0, 10), 1));
2873 
2874       data.push_back(vertex_half_edge(Point(5, 2), Point(5, 5), -1));
2875       data.push_back(vertex_half_edge(Point(5, 2), Point(7, 2), -1));
2876 
2877       data.push_back(vertex_half_edge(Point(5, 5), Point(5, 10), -1));
2878       data.push_back(vertex_half_edge(Point(5, 5), Point(5, 2), 1));
2879       data.push_back(vertex_half_edge(Point(5, 5), Point(10, 5), -1));
2880       data.push_back(vertex_half_edge(Point(5, 5), Point(7, 2), 1));
2881 
2882       data.push_back(vertex_half_edge(Point(7, 2), Point(5, 5), -1));
2883       data.push_back(vertex_half_edge(Point(7, 2), Point(5, 2), 1));
2884 
2885       polygon_sort(data.begin(), data.end());
2886       pf.scan(polys, data.begin(), data.end());
2887       stdcout << "result size: " << polys.size() << "\n";
2888       for(std::size_t i = 0; i < polys.size(); ++i) {
2889         stdcout << polys[i] << "\n";
2890       }
2891       stdcout << "done testing trapezoid formation\n";
2892       return true;
2893     }
2894   };
2895 
2896   template <typename T>
2897   struct PolyLineArbitraryByConcept<T, polygon_with_holes_concept> { typedef poly_line_arbitrary_polygon_data<T> type; };
2898   template <typename T>
2899   struct PolyLineArbitraryByConcept<T, polygon_concept> { typedef poly_line_arbitrary_hole_data<T> type; };
2900 
2901   template <typename T>
2902   struct geometry_concept<poly_line_arbitrary_polygon_data<T> > { typedef polygon_45_with_holes_concept type; };
2903   template <typename T>
2904   struct geometry_concept<poly_line_arbitrary_hole_data<T> > { typedef polygon_45_concept type; };
2905 }
2906 }
2907 #endif
2908