polygon/detail/polygon_arbitrary_formation.hpp

0001 /*
0002     Copyright 2008 Intel Corporation
0003
0004     Use, modification and distribution are subject to the Boost Software License,
0005     Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0006     http://www.boost.org/LICENSE_1_0.txt).
0007 */
0008 #ifndef BOOST_POLYGON_POLYGON_ARBITRARY_FORMATION_HPP
0009 #define BOOST_POLYGON_POLYGON_ARBITRARY_FORMATION_HPP
0010 namespace boost { namespace polygon{
0011   template <typename T, typename T2>
0012   struct PolyLineArbitraryByConcept {};
0013
0014   template <typename T>
0015   class poly_line_arbitrary_polygon_data;
0016   template <typename T>
0017   class poly_line_arbitrary_hole_data;
0018
0019   template <typename Unit>
0020   struct scanline_base {
0021
0022     typedef point_data<Unit> Point;
0023     typedef std::pair<Point, Point> half_edge;
0024
0025     class less_point {
0026     public:
0027       typedef Point first_argument_type;
0028       typedef Point second_argument_type;
0029       typedef bool result_type;
0030       inline less_point() {}
0031       inline bool operator () (const Point& pt1, const Point& pt2) const {
0032         if(pt1.get(HORIZONTAL) < pt2.get(HORIZONTAL)) return true;
0033         if(pt1.get(HORIZONTAL) == pt2.get(HORIZONTAL)) {
0034           if(pt1.get(VERTICAL) < pt2.get(VERTICAL)) return true;
0035         }
0036         return false;
0037       }
0038     };
0039
0040     static inline bool between(Point pt, Point pt1, Point pt2) {
0041       less_point lp;
0042       if(lp(pt1, pt2))
0043         return lp(pt, pt2) && lp(pt1, pt);
0044       return lp(pt, pt1) && lp(pt2, pt);
0045     }
0046
0047     template <typename area_type>
0048     static inline Unit compute_intercept(const area_type& dy2,
0049                                          const area_type& dx1,
0050                                          const area_type& dx2) {
0051       //intercept = dy2 * dx1 / dx2
0052       //return (Unit)(((area_type)dy2 * (area_type)dx1) / (area_type)dx2);
0053       area_type dx1_q = dx1 / dx2;
0054       area_type dx1_r = dx1 % dx2;
0055       return dx1_q * dy2 + (dy2 * dx1_r)/dx2;
0056     }
0057
0058     template <typename area_type>
0059     static inline bool equal_slope(area_type dx1, area_type dy1, area_type dx2, area_type dy2) {
0060       typedef typename coordinate_traits<Unit>::unsigned_area_type unsigned_product_type;
0061       unsigned_product_type cross_1 = (unsigned_product_type)(dx2 < 0 ? -dx2 :dx2) * (unsigned_product_type)(dy1 < 0 ? -dy1 : dy1);
0062       unsigned_product_type cross_2 = (unsigned_product_type)(dx1 < 0 ? -dx1 :dx1) * (unsigned_product_type)(dy2 < 0 ? -dy2 : dy2);
0063       int dx1_sign = dx1 < 0 ? -1 : 1;
0064       int dx2_sign = dx2 < 0 ? -1 : 1;
0065       int dy1_sign = dy1 < 0 ? -1 : 1;
0066       int dy2_sign = dy2 < 0 ? -1 : 1;
0067       int cross_1_sign = dx2_sign * dy1_sign;
0068       int cross_2_sign = dx1_sign * dy2_sign;
0069       return cross_1 == cross_2 && (cross_1_sign == cross_2_sign || cross_1 == 0);
0070     }
0071
0072     template <typename T>
0073     static inline bool equal_slope_hp(const T& dx1, const T& dy1, const T& dx2, const T& dy2) {
0074       return dx1 * dy2 == dx2 * dy1;
0075     }
0076
0077     static inline bool equal_slope(const Unit& x, const Unit& y,
0078                                    const Point& pt1, const Point& pt2) {
0079       const Point* pts[2] = {&pt1, &pt2};
0080       typedef typename coordinate_traits<Unit>::manhattan_area_type at;
0081       at dy2 = (at)pts[1]->get(VERTICAL) - (at)y;
0082       at dy1 = (at)pts[0]->get(VERTICAL) - (at)y;
0083       at dx2 = (at)pts[1]->get(HORIZONTAL) - (at)x;
0084       at dx1 = (at)pts[0]->get(HORIZONTAL) - (at)x;
0085       return equal_slope(dx1, dy1, dx2, dy2);
0086     }
0087
0088     template <typename area_type>
0089     static inline bool less_slope(area_type dx1, area_type dy1, area_type dx2, area_type dy2) {
0090       //reflext x and y slopes to right hand side half plane
0091       if(dx1 < 0) {
0092         dy1 *= -1;
0093         dx1 *= -1;
0094       } else if(dx1 == 0) {
0095         //if the first slope is vertical the first cannot be less
0096         return false;
0097       }
0098       if(dx2 < 0) {
0099         dy2 *= -1;
0100         dx2 *= -1;
0101       } else if(dx2 == 0) {
0102         //if the second slope is vertical the first is always less unless it is also vertical, in which case they are equal
0103         return dx1 != 0;
0104       }
0105       typedef typename coordinate_traits<Unit>::unsigned_area_type unsigned_product_type;
0106       unsigned_product_type cross_1 = (unsigned_product_type)(dx2 < 0 ? -dx2 :dx2) * (unsigned_product_type)(dy1 < 0 ? -dy1 : dy1);
0107       unsigned_product_type cross_2 = (unsigned_product_type)(dx1 < 0 ? -dx1 :dx1) * (unsigned_product_type)(dy2 < 0 ? -dy2 : dy2);
0108       int dx1_sign = dx1 < 0 ? -1 : 1;
0109       int dx2_sign = dx2 < 0 ? -1 : 1;
0110       int dy1_sign = dy1 < 0 ? -1 : 1;
0111       int dy2_sign = dy2 < 0 ? -1 : 1;
0112       int cross_1_sign = dx2_sign * dy1_sign;
0113       int cross_2_sign = dx1_sign * dy2_sign;
0114       if(cross_1_sign < cross_2_sign) return true;
0115       if(cross_2_sign < cross_1_sign) return false;
0116       if(cross_1_sign == -1) return cross_2 < cross_1;
0117       return cross_1 < cross_2;
0118     }
0119
0120     static inline bool less_slope(const Unit& x, const Unit& y,
0121                                   const Point& pt1, const Point& pt2) {
0122       const Point* pts[2] = {&pt1, &pt2};
0123       //compute y value on edge from pt_ to pts[1] at the x value of pts[0]
0124       typedef typename coordinate_traits<Unit>::manhattan_area_type at;
0125       at dy2 = (at)pts[1]->get(VERTICAL) - (at)y;
0126       at dy1 = (at)pts[0]->get(VERTICAL) - (at)y;
0127       at dx2 = (at)pts[1]->get(HORIZONTAL) - (at)x;
0128       at dx1 = (at)pts[0]->get(HORIZONTAL) - (at)x;
0129       return less_slope(dx1, dy1, dx2, dy2);
0130     }
0131
0132     //return -1 below, 0 on and 1 above line
0133     static inline int on_above_or_below(Point pt, const half_edge& he) {
0134       if(pt == he.first || pt == he.second) return 0;
0135       if(equal_slope(pt.get(HORIZONTAL), pt.get(VERTICAL), he.first, he.second)) return 0;
0136       bool less_result = less_slope(pt.get(HORIZONTAL), pt.get(VERTICAL), he.first, he.second);
0137       int retval = less_result ? -1 : 1;
0138       less_point lp;
0139       if(lp(he.second, he.first)) retval *= -1;
0140       if(!between(pt, he.first, he.second)) retval *= -1;
0141       return retval;
0142     }
0143
0144     //returns true is the segment intersects the integer grid square with lower
0145     //left corner at point
0146     static inline bool intersects_grid(Point pt, const half_edge& he) {
0147       if(pt == he.second) return true;
0148       if(pt == he.first) return true;
0149       rectangle_data<Unit> rect1;
0150       set_points(rect1, he.first, he.second);
0151       if(contains(rect1, pt, true)) {
0152         if(is_vertical(he) || is_horizontal(he)) return true;
0153       } else {
0154         return false; //can't intersect a grid not within bounding box
0155       }
0156       Unit x = pt.get(HORIZONTAL);
0157       Unit y = pt.get(VERTICAL);
0158       if(equal_slope(x, y, he.first, he.second) &&
0159          between(pt, he.first, he.second)) return true;
0160       Point pt01(pt.get(HORIZONTAL), pt.get(VERTICAL) + 1);
0161       Point pt10(pt.get(HORIZONTAL) + 1, pt.get(VERTICAL));
0162       Point pt11(pt.get(HORIZONTAL) + 1, pt.get(VERTICAL) + 1);
0163 //       if(pt01 == he.first) return true;
0164 //       if(pt10 == he.first) return true;
0165 //       if(pt11 == he.first) return true;
0166 //       if(pt01 == he.second) return true;
0167 //       if(pt10 == he.second) return true;
0168 //       if(pt11 == he.second) return true;
0169       //check non-integer intersections
0170       half_edge widget1(pt, pt11);
0171       //intersects but not just at pt11
0172       if(intersects(widget1, he) && on_above_or_below(pt11, he)) return true;
0173       half_edge widget2(pt01, pt10);
0174       //intersects but not just at pt01 or 10
0175       if(intersects(widget2, he) && on_above_or_below(pt01, he) && on_above_or_below(pt10, he)) return true;
0176       return false;
0177     }
0178
0179     static inline Unit evalAtXforYlazy(Unit xIn, Point pt, Point other_pt) {
0180       long double
0181         evalAtXforYret, evalAtXforYxIn, evalAtXforYx1, evalAtXforYy1, evalAtXforYdx1, evalAtXforYdx,
0182         evalAtXforYdy, evalAtXforYx2, evalAtXforYy2, evalAtXforY0;
0183       //y = (x - x1)dy/dx + y1
0184       //y = (xIn - pt.x)*(other_pt.y-pt.y)/(other_pt.x-pt.x) + pt.y
0185       //assert pt.x != other_pt.x
0186       if(pt.y() == other_pt.y())
0187         return pt.y();
0188       evalAtXforYxIn = xIn;
0189       evalAtXforYx1 = pt.get(HORIZONTAL);
0190       evalAtXforYy1 = pt.get(VERTICAL);
0191       evalAtXforYdx1 = evalAtXforYxIn - evalAtXforYx1;
0192       evalAtXforY0 = 0;
0193       if(evalAtXforYdx1 == evalAtXforY0) return (Unit)evalAtXforYy1;
0194       evalAtXforYx2 = other_pt.get(HORIZONTAL);
0195       evalAtXforYy2 = other_pt.get(VERTICAL);
0196
0197       evalAtXforYdx = evalAtXforYx2 - evalAtXforYx1;
0198       evalAtXforYdy = evalAtXforYy2 - evalAtXforYy1;
0199       evalAtXforYret = ((evalAtXforYdx1) * evalAtXforYdy / evalAtXforYdx + evalAtXforYy1);
0200       return (Unit)evalAtXforYret;
0201     }
0202
0203     static inline typename high_precision_type<Unit>::type evalAtXforY(Unit xIn, Point pt, Point other_pt) {
0204       typename high_precision_type<Unit>::type
0205         evalAtXforYret, evalAtXforYxIn, evalAtXforYx1, evalAtXforYy1, evalAtXforYdx1, evalAtXforYdx,
0206         evalAtXforYdy, evalAtXforYx2, evalAtXforYy2, evalAtXforY0;
0207       //y = (x - x1)dy/dx + y1
0208       //y = (xIn - pt.x)*(other_pt.y-pt.y)/(other_pt.x-pt.x) + pt.y
0209       //assert pt.x != other_pt.x
0210       typedef typename high_precision_type<Unit>::type high_precision;
0211       if(pt.y() == other_pt.y())
0212         return (high_precision)pt.y();
0213       evalAtXforYxIn = (high_precision)xIn;
0214       evalAtXforYx1 = pt.get(HORIZONTAL);
0215       evalAtXforYy1 = pt.get(VERTICAL);
0216       evalAtXforYdx1 = evalAtXforYxIn - evalAtXforYx1;
0217       evalAtXforY0 = high_precision(0);
0218       if(evalAtXforYdx1 == evalAtXforY0) return evalAtXforYret = evalAtXforYy1;
0219       evalAtXforYx2 = (high_precision)other_pt.get(HORIZONTAL);
0220       evalAtXforYy2 = (high_precision)other_pt.get(VERTICAL);
0221
0222       evalAtXforYdx = evalAtXforYx2 - evalAtXforYx1;
0223       evalAtXforYdy = evalAtXforYy2 - evalAtXforYy1;
0224       evalAtXforYret = ((evalAtXforYdx1) * evalAtXforYdy / evalAtXforYdx + evalAtXforYy1);
0225       return evalAtXforYret;
0226     }
0227
0228     struct evalAtXforYPack {
0229     typename high_precision_type<Unit>::type
0230     evalAtXforYret, evalAtXforYxIn, evalAtXforYx1, evalAtXforYy1, evalAtXforYdx1, evalAtXforYdx,
0231                            evalAtXforYdy, evalAtXforYx2, evalAtXforYy2, evalAtXforY0;
0232       inline const typename high_precision_type<Unit>::type& evalAtXforY(Unit xIn, Point pt, Point other_pt) {
0233         //y = (x - x1)dy/dx + y1
0234         //y = (xIn - pt.x)*(other_pt.y-pt.y)/(other_pt.x-pt.x) + pt.y
0235         //assert pt.x != other_pt.x
0236         typedef typename high_precision_type<Unit>::type high_precision;
0237         if(pt.y() == other_pt.y()) {
0238           evalAtXforYret = (high_precision)pt.y();
0239           return evalAtXforYret;
0240         }
0241         evalAtXforYxIn = (high_precision)xIn;
0242         evalAtXforYx1 = pt.get(HORIZONTAL);
0243         evalAtXforYy1 = pt.get(VERTICAL);
0244         evalAtXforYdx1 = evalAtXforYxIn - evalAtXforYx1;
0245         evalAtXforY0 = high_precision(0);
0246         if(evalAtXforYdx1 == evalAtXforY0) return evalAtXforYret = evalAtXforYy1;
0247         evalAtXforYx2 = (high_precision)other_pt.get(HORIZONTAL);
0248         evalAtXforYy2 = (high_precision)other_pt.get(VERTICAL);
0249
0250         evalAtXforYdx = evalAtXforYx2 - evalAtXforYx1;
0251         evalAtXforYdy = evalAtXforYy2 - evalAtXforYy1;
0252         evalAtXforYret = ((evalAtXforYdx1) * evalAtXforYdy / evalAtXforYdx + evalAtXforYy1);
0253         return evalAtXforYret;
0254       }
0255     };
0256
0257     static inline bool is_vertical(const half_edge& he) {
0258       return he.first.get(HORIZONTAL) == he.second.get(HORIZONTAL);
0259     }
0260
0261     static inline bool is_horizontal(const half_edge& he) {
0262       return he.first.get(VERTICAL) == he.second.get(VERTICAL);
0263     }
0264
0265     static inline bool is_45_degree(const half_edge& he) {
0266       return euclidean_distance(he.first, he.second, HORIZONTAL) == euclidean_distance(he.first, he.second, VERTICAL);
0267     }
0268
0269     //scanline comparator functor
0270     class less_half_edge {
0271     private:
0272       Unit *x_; //x value at which to apply comparison
0273       int *justBefore_;
0274       evalAtXforYPack * pack_;
0275     public:
0276       typedef half_edge first_argument_type;
0277       typedef half_edge second_argument_type;
0278       typedef bool result_type;
0279       inline less_half_edge() : x_(0), justBefore_(0), pack_(0) {}
0280       inline less_half_edge(Unit *x, int *justBefore, evalAtXforYPack * packIn) : x_(x), justBefore_(justBefore), pack_(packIn) {}
0281       inline less_half_edge(const less_half_edge& that) : x_(that.x_), justBefore_(that.justBefore_),
0282                                                           pack_(that.pack_){}
0283       inline less_half_edge& operator=(const less_half_edge& that) {
0284         x_ = that.x_;
0285         justBefore_ = that.justBefore_;
0286         pack_ = that.pack_;
0287         return *this; }
0288       inline bool operator () (const half_edge& elm1, const half_edge& elm2) const {
0289         if((std::max)(elm1.first.y(), elm1.second.y()) < (std::min)(elm2.first.y(), elm2.second.y()))
0290           return true;
0291         if((std::min)(elm1.first.y(), elm1.second.y()) > (std::max)(elm2.first.y(), elm2.second.y()))
0292           return false;
0293
0294         //check if either x of elem1 is equal to x_
0295         Unit localx = *x_;
0296         Unit elm1y = 0;
0297         bool elm1_at_x = false;
0298         if(localx == elm1.first.get(HORIZONTAL)) {
0299           elm1_at_x = true;
0300           elm1y = elm1.first.get(VERTICAL);
0301         } else if(localx == elm1.second.get(HORIZONTAL)) {
0302           elm1_at_x = true;
0303           elm1y = elm1.second.get(VERTICAL);
0304         }
0305         Unit elm2y = 0;
0306         bool elm2_at_x = false;
0307         if(localx == elm2.first.get(HORIZONTAL)) {
0308           elm2_at_x = true;
0309           elm2y = elm2.first.get(VERTICAL);
0310         } else if(localx == elm2.second.get(HORIZONTAL)) {
0311           elm2_at_x = true;
0312           elm2y = elm2.second.get(VERTICAL);
0313         }
0314         bool retval = false;
0315         if(!(elm1_at_x && elm2_at_x)) {
0316           //at least one of the segments doesn't have an end point a the current x
0317           //-1 below, 1 above
0318           int pt1_oab = on_above_or_below(elm1.first, half_edge(elm2.first, elm2.second));
0319           int pt2_oab = on_above_or_below(elm1.second, half_edge(elm2.first, elm2.second));
0320           if(pt1_oab == pt2_oab) {
0321             if(pt1_oab == -1)
0322               retval = true; //pt1 is below elm2 so elm1 is below elm2
0323           } else {
0324             //the segments can't cross so elm2 is on whatever side of elm1 that one of its ends is
0325             int pt3_oab = on_above_or_below(elm2.first, half_edge(elm1.first, elm1.second));
0326             if(pt3_oab == 1)
0327               retval = true; //elm1's point is above elm1
0328           }
0329         } else {
0330           if(elm1y < elm2y) {
0331             retval = true;
0332           } else if(elm1y == elm2y) {
0333             if(elm1 == elm2)
0334               return false;
0335             retval = less_slope(elm1.second.get(HORIZONTAL) - elm1.first.get(HORIZONTAL),
0336                                      elm1.second.get(VERTICAL) - elm1.first.get(VERTICAL),
0337                                      elm2.second.get(HORIZONTAL) - elm2.first.get(HORIZONTAL),
0338                                      elm2.second.get(VERTICAL) - elm2.first.get(VERTICAL));
0339             retval = ((*justBefore_) != 0) ^ retval;
0340           }
0341         }
0342         return retval;
0343       }
0344     };
0345
0346     template <typename unsigned_product_type>
0347     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) {
0348       int switcher = 0;
0349       if(a_sign < 0) switcher += 1;
0350       if(b_sign < 0) switcher += 2;
0351       if(a < b) switcher += 4;
0352       switch (switcher) {
0353       case 0: //both positive
0354         result = a + b;
0355         result_sign = 1;
0356         break;
0357       case 1: //a is negative
0358         result = a - b;
0359         result_sign = -1;
0360         break;
0361       case 2: //b is negative
0362         result = a - b;
0363         result_sign = 1;
0364         break;
0365       case 3: //both negative
0366         result = a + b;
0367         result_sign = -1;
0368         break;
0369       case 4: //both positive
0370         result = a + b;
0371         result_sign = 1;
0372         break;
0373       case 5: //a is negative
0374         result = b - a;
0375         result_sign = 1;
0376         break;
0377       case 6: //b is negative
0378         result = b - a;
0379         result_sign = -1;
0380         break;
0381       case 7: //both negative
0382         result = b + a;
0383         result_sign = -1;
0384         break;
0385       };
0386     }
0387
0388     struct compute_intersection_pack {
0389       typedef typename high_precision_type<Unit>::type high_precision;
0390       high_precision y_high, dx1, dy1, dx2, dy2, x11, x21, y11, y21, x_num, y_num, x_den, y_den, x, y;
0391       static inline bool compute_lazy_intersection(Point& intersection, const half_edge& he1, const half_edge& he2,
0392                                                    bool projected = false, bool round_closest = false) {
0393         long double y_high, dx1, dy1, dx2, dy2, x11, x21, y11, y21, x_num, y_num, x_den, y_den, x, y;
0394         typedef rectangle_data<Unit> Rectangle;
0395         Rectangle rect1, rect2;
0396         set_points(rect1, he1.first, he1.second);
0397         set_points(rect2, he2.first, he2.second);
0398         if(!projected && !::boost::polygon::intersects(rect1, rect2, true)) return false;
0399         if(is_vertical(he1)) {
0400           if(is_vertical(he2)) return false;
0401           y_high = evalAtXforYlazy(he1.first.get(HORIZONTAL), he2.first, he2.second);
0402           Unit y_local = (Unit)y_high;
0403           if(y_high < y_local) --y_local;
0404           if(projected || contains(rect1.get(VERTICAL), y_local, true)) {
0405             intersection = Point(he1.first.get(HORIZONTAL), y_local);
0406             return true;
0407           } else {
0408             return false;
0409           }
0410         } else if(is_vertical(he2)) {
0411           y_high = evalAtXforYlazy(he2.first.get(HORIZONTAL), he1.first, he1.second);
0412           Unit y_local = (Unit)y_high;
0413           if(y_high < y_local) --y_local;
0414           if(projected || contains(rect2.get(VERTICAL), y_local, true)) {
0415             intersection = Point(he2.first.get(HORIZONTAL), y_local);
0416             return true;
0417           } else {
0418             return false;
0419           }
0420         }
0421         //the bounding boxes of the two line segments intersect, so we check closer to find the intersection point
0422         dy2 = (he2.second.get(VERTICAL)) -
0423           (he2.first.get(VERTICAL));
0424         dy1 = (he1.second.get(VERTICAL)) -
0425           (he1.first.get(VERTICAL));
0426         dx2 = (he2.second.get(HORIZONTAL)) -
0427           (he2.first.get(HORIZONTAL));
0428         dx1 = (he1.second.get(HORIZONTAL)) -
0429           (he1.first.get(HORIZONTAL));
0430         if(equal_slope_hp(dx1, dy1, dx2, dy2)) return false;
0431         //the line segments have different slopes
0432         //we can assume that the line segments are not vertical because such an intersection is handled elsewhere
0433         x11 = (he1.first.get(HORIZONTAL));
0434         x21 = (he2.first.get(HORIZONTAL));
0435         y11 = (he1.first.get(VERTICAL));
0436         y21 = (he2.first.get(VERTICAL));
0437         //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);
0438         //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);
0439         x_num = (x11 * dy1 * dx2 - x21 * dy2 * dx1 + y21 * dx1 * dx2 - y11 * dx1 * dx2);
0440         x_den = (dy1 * dx2 - dy2 * dx1);
0441         y_num = (y11 * dx1 * dy2 - y21 * dx2 * dy1 + x21 * dy1 * dy2 - x11 * dy1 * dy2);
0442         y_den = (dx1 * dy2 - dx2 * dy1);
0443         x = x_num / x_den;
0444         y = y_num / y_den;
0445         //std::cout << "cross1 " << dy1 << " " << dx2 << " " << dy1 * dx2 << "\n";
0446         //std::cout << "cross2 " << dy2 << " " << dx1 << " " << dy2 * dx1 << "\n";
0447         //Unit exp_x = compute_x_intercept<at>(x11, x21, y11, y21, dy1, dy2, dx1, dx2);
0448         //Unit exp_y = compute_x_intercept<at>(y11, y21, x11, x21, dx1, dx2, dy1, dy2);
0449         if(round_closest) {
0450           x = x + 0.5;
0451           y = y + 0.5;
0452         }
0453         Unit x_unit = (Unit)(x);
0454         Unit y_unit = (Unit)(y);
0455         //truncate downward if it went up due to negative number
0456         if(x < x_unit) --x_unit;
0457         if(y < y_unit) --y_unit;
0458         if(is_horizontal(he1))
0459           y_unit = he1.first.y();
0460         if(is_horizontal(he2))
0461           y_unit = he2.first.y();
0462         //if(x != exp_x || y != exp_y)
0463         //  std::cout << exp_x << " " << exp_y << " " << x << " " << y << "\n";
0464         //Unit y1 = evalAtXforY(exp_x, he1.first, he1.second);
0465         //Unit y2 = evalAtXforY(exp_x, he2.first, he2.second);
0466         //std::cout << exp_x << " " << exp_y << " " << y1 << " " << y2 << "\n";
0467         Point result(x_unit, y_unit);
0468         if(!projected && !contains(rect1, result, true)) return false;
0469         if(!projected && !contains(rect2, result, true)) return false;
0470         if(projected) {
0471           rectangle_data<long double> inf_rect(-(long double)(std::numeric_limits<Unit>::max)(),
0472                                                -(long double) (std::numeric_limits<Unit>::max)(),
0473                                                (long double)(std::numeric_limits<Unit>::max)(),
0474                                                (long double) (std::numeric_limits<Unit>::max)() );
0475           if(contains(inf_rect, point_data<long double>(x, y), true)) {
0476             intersection = result;
0477             return true;
0478           } else
0479             return false;
0480         }
0481         intersection = result;
0482         return true;
0483       }
0484
0485       inline bool compute_intersection(Point& intersection, const half_edge& he1, const half_edge& he2,
0486                                        bool projected = false, bool round_closest = false) {
0487         if(!projected && !intersects(he1, he2))
0488            return false;
0489         bool lazy_success = compute_lazy_intersection(intersection, he1, he2, projected);
0490         if(!projected) {
0491           if(lazy_success) {
0492             if(intersects_grid(intersection, he1) &&
0493                intersects_grid(intersection, he2))
0494               return true;
0495           }
0496         } else {
0497           return lazy_success;
0498         }
0499         return compute_exact_intersection(intersection, he1, he2, projected, round_closest);
0500       }
0501
0502       inline bool compute_exact_intersection(Point& intersection, const half_edge& he1, const half_edge& he2,
0503                                              bool projected = false, bool round_closest = false) {
0504         if(!projected && !intersects(he1, he2))
0505            return false;
0506         typedef rectangle_data<Unit> Rectangle;
0507         Rectangle rect1, rect2;
0508         set_points(rect1, he1.first, he1.second);
0509         set_points(rect2, he2.first, he2.second);
0510         if(!::boost::polygon::intersects(rect1, rect2, true)) return false;
0511         if(is_vertical(he1)) {
0512           if(is_vertical(he2)) return false;
0513           y_high = evalAtXforY(he1.first.get(HORIZONTAL), he2.first, he2.second);
0514           Unit y = convert_high_precision_type<Unit>(y_high);
0515           if(y_high < (high_precision)y) --y;
0516           if(contains(rect1.get(VERTICAL), y, true)) {
0517             intersection = Point(he1.first.get(HORIZONTAL), y);
0518             return true;
0519           } else {
0520             return false;
0521           }
0522         } else if(is_vertical(he2)) {
0523           y_high = evalAtXforY(he2.first.get(HORIZONTAL), he1.first, he1.second);
0524           Unit y = convert_high_precision_type<Unit>(y_high);
0525           if(y_high < (high_precision)y) --y;
0526           if(contains(rect2.get(VERTICAL), y, true)) {
0527             intersection = Point(he2.first.get(HORIZONTAL), y);
0528             return true;
0529           } else {
0530             return false;
0531           }
0532         }
0533         //the bounding boxes of the two line segments intersect, so we check closer to find the intersection point
0534         dy2 = (high_precision)(he2.second.get(VERTICAL)) -
0535           (high_precision)(he2.first.get(VERTICAL));
0536         dy1 = (high_precision)(he1.second.get(VERTICAL)) -
0537           (high_precision)(he1.first.get(VERTICAL));
0538         dx2 = (high_precision)(he2.second.get(HORIZONTAL)) -
0539           (high_precision)(he2.first.get(HORIZONTAL));
0540         dx1 = (high_precision)(he1.second.get(HORIZONTAL)) -
0541           (high_precision)(he1.first.get(HORIZONTAL));
0542         if(equal_slope_hp(dx1, dy1, dx2, dy2)) return false;
0543         //the line segments have different slopes
0544         //we can assume that the line segments are not vertical because such an intersection is handled elsewhere
0545         x11 = (high_precision)(he1.first.get(HORIZONTAL));
0546         x21 = (high_precision)(he2.first.get(HORIZONTAL));
0547         y11 = (high_precision)(he1.first.get(VERTICAL));
0548         y21 = (high_precision)(he2.first.get(VERTICAL));
0549         //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);
0550         //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);
0551         x_num = (x11 * dy1 * dx2 - x21 * dy2 * dx1 + y21 * dx1 * dx2 - y11 * dx1 * dx2);
0552         x_den = (dy1 * dx2 - dy2 * dx1);
0553         y_num = (y11 * dx1 * dy2 - y21 * dx2 * dy1 + x21 * dy1 * dy2 - x11 * dy1 * dy2);
0554         y_den = (dx1 * dy2 - dx2 * dy1);
0555         x = x_num / x_den;
0556         y = y_num / y_den;
0557         //std::cout << x << " " << y << "\n";
0558         //std::cout << "cross1 " << dy1 << " " << dx2 << " " << dy1 * dx2 << "\n";
0559         //std::cout << "cross2 " << dy2 << " " << dx1 << " " << dy2 * dx1 << "\n";
0560         //Unit exp_x = compute_x_intercept<at>(x11, x21, y11, y21, dy1, dy2, dx1, dx2);
0561         //Unit exp_y = compute_x_intercept<at>(y11, y21, x11, x21, dx1, dx2, dy1, dy2);
0562         if(round_closest) {
0563           x = x + (high_precision)0.5;
0564           y = y + (high_precision)0.5;
0565         }
0566         Unit x_unit = convert_high_precision_type<Unit>(x);
0567         Unit y_unit = convert_high_precision_type<Unit>(y);
0568         //truncate downward if it went up due to negative number
0569         if(x < (high_precision)x_unit) --x_unit;
0570         if(y < (high_precision)y_unit) --y_unit;
0571         if(is_horizontal(he1))
0572           y_unit = he1.first.y();
0573         if(is_horizontal(he2))
0574           y_unit = he2.first.y();
0575         //if(x != exp_x || y != exp_y)
0576         //  std::cout << exp_x << " " << exp_y << " " << x << " " << y << "\n";
0577         //Unit y1 = evalAtXforY(exp_x, he1.first, he1.second);
0578         //Unit y2 = evalAtXforY(exp_x, he2.first, he2.second);
0579         //std::cout << exp_x << " " << exp_y << " " << y1 << " " << y2 << "\n";
0580         Point result(x_unit, y_unit);
0581         if(!contains(rect1, result, true)) return false;
0582         if(!contains(rect2, result, true)) return false;
0583         if(projected) {
0584           high_precision b1 = (high_precision) (std::numeric_limits<Unit>::min)();
0585           high_precision b2 = (high_precision) (std::numeric_limits<Unit>::max)();
0586           if(x > b2 || y > b2 || x < b1 || y < b1)
0587             return false;
0588         }
0589         intersection = result;
0590         return true;
0591       }
0592     };
0593
0594     static inline bool compute_intersection(Point& intersection, const half_edge& he1, const half_edge& he2) {
0595       typedef typename high_precision_type<Unit>::type high_precision;
0596       typedef rectangle_data<Unit> Rectangle;
0597       Rectangle rect1, rect2;
0598       set_points(rect1, he1.first, he1.second);
0599       set_points(rect2, he2.first, he2.second);
0600       if(!::boost::polygon::intersects(rect1, rect2, true)) return false;
0601       if(is_vertical(he1)) {
0602         if(is_vertical(he2)) return false;
0603         high_precision y_high = evalAtXforY(he1.first.get(HORIZONTAL), he2.first, he2.second);
0604         Unit y = convert_high_precision_type<Unit>(y_high);
0605         if(y_high < (high_precision)y) --y;
0606         if(contains(rect1.get(VERTICAL), y, true)) {
0607           intersection = Point(he1.first.get(HORIZONTAL), y);
0608           return true;
0609         } else {
0610           return false;
0611         }
0612       } else if(is_vertical(he2)) {
0613         high_precision y_high = evalAtXforY(he2.first.get(HORIZONTAL), he1.first, he1.second);
0614         Unit y = convert_high_precision_type<Unit>(y_high);
0615         if(y_high < (high_precision)y) --y;
0616         if(contains(rect2.get(VERTICAL), y, true)) {
0617           intersection = Point(he2.first.get(HORIZONTAL), y);
0618           return true;
0619         } else {
0620           return false;
0621         }
0622       }
0623       //the bounding boxes of the two line segments intersect, so we check closer to find the intersection point
0624       high_precision dy2 = (high_precision)(he2.second.get(VERTICAL)) -
0625         (high_precision)(he2.first.get(VERTICAL));
0626       high_precision dy1 = (high_precision)(he1.second.get(VERTICAL)) -
0627         (high_precision)(he1.first.get(VERTICAL));
0628       high_precision dx2 = (high_precision)(he2.second.get(HORIZONTAL)) -
0629         (high_precision)(he2.first.get(HORIZONTAL));
0630       high_precision dx1 = (high_precision)(he1.second.get(HORIZONTAL)) -
0631         (high_precision)(he1.first.get(HORIZONTAL));
0632       if(equal_slope_hp(dx1, dy1, dx2, dy2)) return false;
0633       //the line segments have different slopes
0634       //we can assume that the line segments are not vertical because such an intersection is handled elsewhere
0635       high_precision x11 = (high_precision)(he1.first.get(HORIZONTAL));
0636       high_precision x21 = (high_precision)(he2.first.get(HORIZONTAL));
0637       high_precision y11 = (high_precision)(he1.first.get(VERTICAL));
0638       high_precision y21 = (high_precision)(he2.first.get(VERTICAL));
0639       //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);
0640       //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);
0641       high_precision x_num = (x11 * dy1 * dx2 - x21 * dy2 * dx1 + y21 * dx1 * dx2 - y11 * dx1 * dx2);
0642       high_precision x_den = (dy1 * dx2 - dy2 * dx1);
0643       high_precision y_num = (y11 * dx1 * dy2 - y21 * dx2 * dy1 + x21 * dy1 * dy2 - x11 * dy1 * dy2);
0644       high_precision y_den = (dx1 * dy2 - dx2 * dy1);
0645       high_precision x = x_num / x_den;
0646       high_precision y = y_num / y_den;
0647       //std::cout << "cross1 " << dy1 << " " << dx2 << " " << dy1 * dx2 << "\n";
0648       //std::cout << "cross2 " << dy2 << " " << dx1 << " " << dy2 * dx1 << "\n";
0649       //Unit exp_x = compute_x_intercept<at>(x11, x21, y11, y21, dy1, dy2, dx1, dx2);
0650       //Unit exp_y = compute_x_intercept<at>(y11, y21, x11, x21, dx1, dx2, dy1, dy2);
0651       Unit x_unit = convert_high_precision_type<Unit>(x);
0652       Unit y_unit = convert_high_precision_type<Unit>(y);
0653       //truncate downward if it went up due to negative number
0654       if(x < (high_precision)x_unit) --x_unit;
0655       if(y < (high_precision)y_unit) --y_unit;
0656       if(is_horizontal(he1))
0657         y_unit = he1.first.y();
0658       if(is_horizontal(he2))
0659         y_unit = he2.first.y();
0660       //if(x != exp_x || y != exp_y)
0661       //  std::cout << exp_x << " " << exp_y << " " << x << " " << y << "\n";
0662       //Unit y1 = evalAtXforY(exp_x, he1.first, he1.second);
0663       //Unit y2 = evalAtXforY(exp_x, he2.first, he2.second);
0664       //std::cout << exp_x << " " << exp_y << " " << y1 << " " << y2 << "\n";
0665       Point result(x_unit, y_unit);
0666       if(!contains(rect1, result, true)) return false;
0667       if(!contains(rect2, result, true)) return false;
0668       intersection = result;
0669       return true;
0670     }
0671
0672     static inline bool intersects(const half_edge& he1, const half_edge& he2) {
0673       typedef rectangle_data<Unit> Rectangle;
0674       Rectangle rect1, rect2;
0675       set_points(rect1, he1.first, he1.second);
0676       set_points(rect2, he2.first, he2.second);
0677       if(::boost::polygon::intersects(rect1, rect2, false)) {
0678         if(he1.first == he2.first) {
0679           if(he1.second != he2.second && equal_slope(he1.first.get(HORIZONTAL), he1.first.get(VERTICAL),
0680                                                      he1.second, he2.second)) {
0681             return true;
0682           } else {
0683             return false;
0684           }
0685         }
0686         if(he1.first == he2.second) {
0687           if(he1.second != he2.first && equal_slope(he1.first.get(HORIZONTAL), he1.first.get(VERTICAL),
0688                                                     he1.second, he2.first)) {
0689             return true;
0690           } else {
0691             return false;
0692           }
0693         }
0694         if(he1.second == he2.first) {
0695           if(he1.first != he2.second && equal_slope(he1.second.get(HORIZONTAL), he1.second.get(VERTICAL),
0696                                                     he1.first, he2.second)) {
0697             return true;
0698           } else {
0699             return false;
0700           }
0701         }
0702         if(he1.second == he2.second) {
0703           if(he1.first != he2.first && equal_slope(he1.second.get(HORIZONTAL), he1.second.get(VERTICAL),
0704                                                    he1.first, he2.first)) {
0705             return true;
0706           } else {
0707             return false;
0708           }
0709         }
0710         int oab1 = on_above_or_below(he1.first, he2);
0711         if(oab1 == 0 && between(he1.first, he2.first, he2.second)) return true;
0712         int oab2 = on_above_or_below(he1.second, he2);
0713         if(oab2 == 0 && between(he1.second, he2.first, he2.second)) return true;
0714         if(oab1 == oab2 && oab1 != 0) return false; //both points of he1 are on same side of he2
0715         int oab3 = on_above_or_below(he2.first, he1);
0716         if(oab3 == 0 && between(he2.first, he1.first, he1.second)) return true;
0717         int oab4 = on_above_or_below(he2.second, he1);
0718         if(oab4 == 0 && between(he2.second, he1.first, he1.second)) return true;
0719         if(oab3 == oab4) return false; //both points of he2 are on same side of he1
0720         return true; //they must cross
0721       }
0722       if(is_vertical(he1) && is_vertical(he2) && he1.first.get(HORIZONTAL) == he2.first.get(HORIZONTAL))
0723         return ::boost::polygon::intersects(rect1.get(VERTICAL), rect2.get(VERTICAL), false) &&
0724           rect1.get(VERTICAL) != rect2.get(VERTICAL);
0725       if(is_horizontal(he1) && is_horizontal(he2) && he1.first.get(VERTICAL) == he2.first.get(VERTICAL))
0726         return ::boost::polygon::intersects(rect1.get(HORIZONTAL), rect2.get(HORIZONTAL), false) &&
0727           rect1.get(HORIZONTAL) != rect2.get(HORIZONTAL);
0728       return false;
0729     }
0730
0731     class vertex_half_edge {
0732     public:
0733       typedef typename high_precision_type<Unit>::type high_precision;
0734       Point pt;
0735       Point other_pt; // 1, 0 or -1
0736       int count; //dxdydTheta
0737       inline vertex_half_edge() : pt(), other_pt(), count() {}
0738       inline vertex_half_edge(const Point& point, const Point& other_point, int countIn) : pt(point), other_pt(other_point), count(countIn) {}
0739       inline vertex_half_edge(const vertex_half_edge& vertex) : pt(vertex.pt), other_pt(vertex.other_pt), count(vertex.count) {}
0740       inline vertex_half_edge& operator=(const vertex_half_edge& vertex){
0741         pt = vertex.pt; other_pt = vertex.other_pt; count = vertex.count; return *this; }
0742       inline bool operator==(const vertex_half_edge& vertex) const {
0743         return pt == vertex.pt && other_pt == vertex.other_pt && count == vertex.count; }
0744       inline bool operator!=(const vertex_half_edge& vertex) const { return !((*this) == vertex); }
0745       inline bool operator<(const vertex_half_edge& vertex) const {
0746         if(pt.get(HORIZONTAL) < vertex.pt.get(HORIZONTAL)) return true;
0747         if(pt.get(HORIZONTAL) == vertex.pt.get(HORIZONTAL)) {
0748           if(pt.get(VERTICAL) < vertex.pt.get(VERTICAL)) return true;
0749           if(pt.get(VERTICAL) == vertex.pt.get(VERTICAL)) { return less_slope(pt.get(HORIZONTAL), pt.get(VERTICAL),
0750                                                                               other_pt, vertex.other_pt);
0751           }
0752         }
0753         return false;
0754       }
0755       inline bool operator>(const vertex_half_edge& vertex) const { return vertex < (*this); }
0756       inline bool operator<=(const vertex_half_edge& vertex) const { return !((*this) > vertex); }
0757       inline bool operator>=(const vertex_half_edge& vertex) const { return !((*this) < vertex); }
0758       inline high_precision evalAtX(Unit xIn) const { return evalAtXforYlazy(xIn, pt, other_pt); }
0759       inline bool is_vertical() const {
0760         return pt.get(HORIZONTAL) == other_pt.get(HORIZONTAL);
0761       }
0762       inline bool is_begin() const {
0763         return pt.get(HORIZONTAL) < other_pt.get(HORIZONTAL) ||
0764           (pt.get(HORIZONTAL) == other_pt.get(HORIZONTAL) &&
0765            (pt.get(VERTICAL) < other_pt.get(VERTICAL)));
0766       }
0767     };
0768
0769     //when scanning Vertex45 for polygon formation we need a scanline comparator functor
0770     class less_vertex_half_edge {
0771     private:
0772       Unit *x_; //x value at which to apply comparison
0773       int *justBefore_;
0774     public:
0775       typedef vertex_half_edge first_argument_type;
0776       typedef vertex_half_edge second_argument_type;
0777       typedef bool result_type;
0778       inline less_vertex_half_edge() : x_(0), justBefore_(0) {}
0779       inline less_vertex_half_edge(Unit *x, int *justBefore) : x_(x), justBefore_(justBefore) {}
0780       inline less_vertex_half_edge(const less_vertex_half_edge& that) : x_(that.x_), justBefore_(that.justBefore_) {}
0781       inline less_vertex_half_edge& operator=(const less_vertex_half_edge& that) { x_ = that.x_; justBefore_ = that.justBefore_; return *this; }
0782       inline bool operator () (const vertex_half_edge& elm1, const vertex_half_edge& elm2) const {
0783         if((std::max)(elm1.pt.y(), elm1.other_pt.y()) < (std::min)(elm2.pt.y(), elm2.other_pt.y()))
0784           return true;
0785         if((std::min)(elm1.pt.y(), elm1.other_pt.y()) > (std::max)(elm2.pt.y(), elm2.other_pt.y()))
0786           return false;
0787         //check if either x of elem1 is equal to x_
0788         Unit localx = *x_;
0789         Unit elm1y = 0;
0790         bool elm1_at_x = false;
0791         if(localx == elm1.pt.get(HORIZONTAL)) {
0792           elm1_at_x = true;
0793           elm1y = elm1.pt.get(VERTICAL);
0794         } else if(localx == elm1.other_pt.get(HORIZONTAL)) {
0795           elm1_at_x = true;
0796           elm1y = elm1.other_pt.get(VERTICAL);
0797         }
0798         Unit elm2y = 0;
0799         bool elm2_at_x = false;
0800         if(localx == elm2.pt.get(HORIZONTAL)) {
0801           elm2_at_x = true;
0802           elm2y = elm2.pt.get(VERTICAL);
0803         } else if(localx == elm2.other_pt.get(HORIZONTAL)) {
0804           elm2_at_x = true;
0805           elm2y = elm2.other_pt.get(VERTICAL);
0806         }
0807         bool retval = false;
0808         if(!(elm1_at_x && elm2_at_x)) {
0809           //at least one of the segments doesn't have an end point a the current x
0810           //-1 below, 1 above
0811           int pt1_oab = on_above_or_below(elm1.pt, half_edge(elm2.pt, elm2.other_pt));
0812           int pt2_oab = on_above_or_below(elm1.other_pt, half_edge(elm2.pt, elm2.other_pt));
0813           if(pt1_oab == pt2_oab) {
0814             if(pt1_oab == -1)
0815               retval = true; //pt1 is below elm2 so elm1 is below elm2
0816           } else {
0817             //the segments can't cross so elm2 is on whatever side of elm1 that one of its ends is
0818             int pt3_oab = on_above_or_below(elm2.pt, half_edge(elm1.pt, elm1.other_pt));
0819             if(pt3_oab == 1)
0820               retval = true; //elm1's point is above elm1
0821           }
0822         } else {
0823           if(elm1y < elm2y) {
0824             retval = true;
0825           } else if(elm1y == elm2y) {
0826             if(elm1.pt == elm2.pt && elm1.other_pt == elm2.other_pt)
0827               return false;
0828             retval = less_slope(elm1.other_pt.get(HORIZONTAL) - elm1.pt.get(HORIZONTAL),
0829                                      elm1.other_pt.get(VERTICAL) - elm1.pt.get(VERTICAL),
0830                                      elm2.other_pt.get(HORIZONTAL) - elm2.pt.get(HORIZONTAL),
0831                                      elm2.other_pt.get(VERTICAL) - elm2.pt.get(VERTICAL));
0832             retval = ((*justBefore_) != 0) ^ retval;
0833           }
0834         }
0835         return retval;
0836       }
0837     };
0838
0839   };
0840
0841   template <typename Unit>
0842   class polygon_arbitrary_formation : public scanline_base<Unit> {
0843   public:
0844     typedef typename scanline_base<Unit>::Point Point;
0845     typedef typename scanline_base<Unit>::half_edge half_edge;
0846     typedef typename scanline_base<Unit>::vertex_half_edge vertex_half_edge;
0847     typedef typename scanline_base<Unit>::less_vertex_half_edge less_vertex_half_edge;
0848
0849     class poly_line_arbitrary {
0850     public:
0851       typedef typename std::list<Point>::const_iterator iterator;
0852
0853       // default constructor of point does not initialize x and y
0854       inline poly_line_arbitrary() : points() {} //do nothing default constructor
0855
0856       // initialize a polygon from x,y values, it is assumed that the first is an x
0857       // and that the input is a well behaved polygon
0858       template<class iT>
0859       inline poly_line_arbitrary& set(iT inputBegin, iT inputEnd) {
0860         points.clear();  //just in case there was some old data there
0861         while(inputBegin != inputEnd) {
0862           points.insert(points.end(), *inputBegin);
0863           ++inputBegin;
0864         }
0865         return *this;
0866       }
0867
0868       // copy constructor (since we have dynamic memory)
0869       inline poly_line_arbitrary(const poly_line_arbitrary& that) : points(that.points) {}
0870
0871       // assignment operator (since we have dynamic memory do a deep copy)
0872       inline poly_line_arbitrary& operator=(const poly_line_arbitrary& that) {
0873         points = that.points;
0874         return *this;
0875       }
0876
0877       // get begin iterator, returns a pointer to a const Unit
0878       inline iterator begin() const { return points.begin(); }
0879
0880       // get end iterator, returns a pointer to a const Unit
0881       inline iterator end() const { return points.end(); }
0882
0883       inline std::size_t size() const { return points.size(); }
0884
0885       //public data member
0886       std::list<Point> points;
0887     };
0888
0889     class active_tail_arbitrary {
0890     protected:
0891       //data
0892       poly_line_arbitrary* tailp_;
0893       active_tail_arbitrary *otherTailp_;
0894       std::list<active_tail_arbitrary*> holesList_;
0895       bool head_;
0896     public:
0897
0898       /**
0899        * @brief iterator over coordinates of the figure
0900        */
0901       typedef typename poly_line_arbitrary::iterator iterator;
0902
0903       /**
0904        * @brief iterator over holes contained within the figure
0905        */
0906       typedef typename std::list<active_tail_arbitrary*>::const_iterator iteratorHoles;
0907
0908       //default constructor
0909       inline active_tail_arbitrary() : tailp_(), otherTailp_(), holesList_(), head_() {}
0910
0911       //constructor
0912       inline active_tail_arbitrary(const vertex_half_edge& vertex, active_tail_arbitrary* otherTailp = 0) : tailp_(), otherTailp_(), holesList_(), head_() {
0913         tailp_ = new poly_line_arbitrary;
0914         tailp_->points.push_back(vertex.pt);
0915         //bool headArray[4] = {false, true, true, true};
0916         bool inverted = vertex.count == -1;
0917         head_ = (!vertex.is_vertical) ^ inverted;
0918         otherTailp_ = otherTailp;
0919       }
0920
0921       inline active_tail_arbitrary(Point point, active_tail_arbitrary* otherTailp, bool head = true) :
0922         tailp_(), otherTailp_(), holesList_(), head_() {
0923         tailp_ = new poly_line_arbitrary;
0924         tailp_->points.push_back(point);
0925         head_ = head;
0926         otherTailp_ = otherTailp;
0927
0928       }
0929       inline active_tail_arbitrary(active_tail_arbitrary* otherTailp) :
0930         tailp_(), otherTailp_(), holesList_(), head_() {
0931         tailp_ = otherTailp->tailp_;
0932         otherTailp_ = otherTailp;
0933       }
0934
0935       //copy constructor
0936       inline active_tail_arbitrary(const active_tail_arbitrary& that) :
0937         tailp_(), otherTailp_(), holesList_(), head_() { (*this) = that; }
0938
0939       //destructor
0940       inline ~active_tail_arbitrary() {
0941         destroyContents();
0942       }
0943
0944       //assignment operator
0945       inline active_tail_arbitrary& operator=(const active_tail_arbitrary& that) {
0946         tailp_ = new poly_line_arbitrary(*(that.tailp_));
0947         head_ = that.head_;
0948         otherTailp_ = that.otherTailp_;
0949         holesList_ = that.holesList_;
0950         return *this;
0951       }
0952
0953       //equivalence operator
0954       inline bool operator==(const active_tail_arbitrary& b) const {
0955         return tailp_ == b.tailp_ && head_ == b.head_;
0956       }
0957
0958       /**
0959        * @brief get the pointer to the polyline that this is an active tail of
0960        */
0961       inline poly_line_arbitrary* getTail() const { return tailp_; }
0962
0963       /**
0964        * @brief get the pointer to the polyline at the other end of the chain
0965        */
0966       inline poly_line_arbitrary* getOtherTail() const { return otherTailp_->tailp_; }
0967
0968       /**
0969        * @brief get the pointer to the activetail at the other end of the chain
0970        */
0971       inline active_tail_arbitrary* getOtherActiveTail() const { return otherTailp_; }
0972
0973       /**
0974        * @brief test if another active tail is the other end of the chain
0975        */
0976       inline bool isOtherTail(const active_tail_arbitrary& b) const { return &b == otherTailp_; }
0977
0978       /**
0979        * @brief update this end of chain pointer to new polyline
0980        */
0981       inline active_tail_arbitrary& updateTail(poly_line_arbitrary* newTail) { tailp_ = newTail; return *this; }
0982
0983       inline bool join(active_tail_arbitrary* tail) {
0984         if(tail == otherTailp_) {
0985           //std::cout << "joining to other tail!\n";
0986           return false;
0987         }
0988         if(tail->head_ == head_) {
0989           //std::cout << "joining head to head!\n";
0990           return false;
0991         }
0992         if(!tailp_) {
0993           //std::cout << "joining empty tail!\n";
0994           return false;
0995         }
0996         if(!(otherTailp_->head_)) {
0997           otherTailp_->copyHoles(*tail);
0998           otherTailp_->copyHoles(*this);
0999         } 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        */
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        */
1037       inline const std::list<active_tail_arbitrary*>& getHoles() const { return holesList_; }
1038
1039       /**
1040        * @brief copy holes from that to this
1041        */
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        */
1047       inline bool solidToRight() const { return !head_; }
1048       inline bool solidToLeft() const { return head_; }
1049
1050       /**
1051        * @brief get vertex
1052        */
1053       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        */
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>
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
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
1176       inline void print() {
1177         //std::cout << this << " " << tailp_ << " " << otherTailp_ << " " << holesList_.size() << " " << head_ << "\n";
1178       }
1179
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;
1214       inline less_half_edge_count() : pt_() {}
1215       inline less_half_edge_count(Point point) : pt_(point) {}
1216       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
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;
1235       inline less_incoming_count() : pt_() {}
1236       inline less_incoming_count(Point point) : pt_(point) {}
1237       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
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
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;
1280       inline vertex_arbitrary_compact() : pt(), count() {}
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       }
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       }
1287       inline vertex_arbitrary_compact(const vertex_arbitrary_compact& vertex) : pt(vertex.pt), count(vertex.count) {}
1288       inline vertex_arbitrary_compact& operator=(const vertex_arbitrary_compact& vertex){
1289         pt = vertex.pt; count = vertex.count; return *this; }
1290       inline bool operator==(const vertex_arbitrary_compact& vertex) const {
1291         return pt == vertex.pt && count == vertex.count; }
1292       inline bool operator!=(const vertex_arbitrary_compact& vertex) const { return !((*this) == vertex); }
1293       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       }
1300       inline bool operator>(const vertex_arbitrary_compact& vertex) const { return vertex < (*this); }
1301       inline bool operator<=(const vertex_arbitrary_compact& vertex) const { return !((*this) > vertex); }
1302       inline bool operator>=(const vertex_arbitrary_compact& vertex) const { return !((*this) < vertex); }
1303       inline bool have_vertex_half_edge(int index) const { return count[index]; }
1304       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:
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     }
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     }
1334     inline polygon_arbitrary_formation(const polygon_arbitrary_formation& that) :
1335       scanData_(), x_((std::numeric_limits<Unit>::min)()), justBefore_(false), fractureHoles_(0) { (*this) = 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>
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>
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
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
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>
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
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>
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>
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>
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>
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>
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>
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>
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>
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>
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;
2138     inline poly_line_arbitrary_hole_data() : p_(0) {}
2139     inline poly_line_arbitrary_hole_data(active_tail_arbitrary* p) : p_(p) {}
2140     //use default copy and assign
2141     inline iterator begin() const { return p_->getTail()->begin(); }
2142     inline iterator end() const { return p_->getTail()->end(); }
2143     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
2171       inline iterator_holes_type() : hole_(), itr_() {}
2172       inline iterator_holes_type(typename active_tail_arbitrary::iteratorHoles itr) : hole_(), itr_(itr) {}
2173       inline iterator_holes_type(const iterator_holes_type& that) : hole_(that.hole_), itr_(that.itr_) {}
2174       inline iterator_holes_type& operator=(const iterator_holes_type& that) {
2175         itr_ = that.itr_;
2176         return *this;
2177       }
2178       inline bool operator==(const iterator_holes_type& that) { return itr_ == that.itr_; }
2179       inline bool operator!=(const iterator_holes_type& that) { return itr_ != that.itr_; }
2180       inline iterator_holes_type& operator++() {
2181         ++itr_;
2182         return *this;
2183       }
2184       inline const iterator_holes_type operator++(int) {
2185         iterator_holes_type tmp = *this;
2186         ++(*this);
2187         return tmp;
2188       }
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
2197     inline poly_line_arbitrary_polygon_data() : p_(0) {}
2198     inline poly_line_arbitrary_polygon_data(active_tail_arbitrary* p) : p_(p) {}
2199     //use default copy and assign
2200     inline iterator_type begin() const { return p_->getTail()->begin(); }
2201     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(); }
2204     inline iterator_holes_type begin_holes() const { return iterator_holes_type(p_->getHoles().begin()); }
2205     inline iterator_holes_type end_holes() const { return iterator_holes_type(p_->getHoles().end()); }
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
2208     inline std::size_t size_holes() const { return 0; }
2209     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:
2242     inline trapezoid_arbitrary_formation() : polygon_arbitrary_formation<Unit>() {}
2243     inline trapezoid_arbitrary_formation(const trapezoid_arbitrary_formation& that) : polygon_arbitrary_formation<Unit>(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>
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
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*>
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
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
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>
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>
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>
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>
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>
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>
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