[Boost-commit] svn:boost r55988 - sandbox/gtl/boost/polygon/detail

Subject: [Boost-commit] svn:boost r55988 - sandbox/gtl/boost/polygon/detail
From: lucanus.j.simonson_at_[hidden]
Date: 2009-09-03 00:28:59

Author: ljsimons
Date: 2009-09-03 00:28:56 EDT (Thu, 03 Sep 2009)
New Revision: 55988
URL: http://svn.boost.org/trac/boost/changeset/55988

Log:
improve genericity of the boolean to support floating and exact coordinate types
Text files modified: 
   sandbox/gtl/boost/polygon/detail/polygon_arbitrary_formation.hpp |   173 +++++++++++++++++++++++++++------------ 
   sandbox/gtl/boost/polygon/detail/scan_arbitrary.hpp              |   136 ++++++++++++++++++++-----------         
   2 files changed, 209 insertions(+), 100 deletions(-)

Modified: sandbox/gtl/boost/polygon/detail/polygon_arbitrary_formation.hpp
==============================================================================
--- sandbox/gtl/boost/polygon/detail/polygon_arbitrary_formation.hpp	(original)
+++ sandbox/gtl/boost/polygon/detail/polygon_arbitrary_formation.hpp	2009-09-03 00:28:56 EDT (Thu, 03 Sep 2009)
@@ -8,6 +8,7 @@
 #ifndef BOOST_POLYGON_POLYGON_ARBITRARY_FORMATION_HPP
 #define BOOST_POLYGON_POLYGON_ARBITRARY_FORMATION_HPP
 namespace boost { namespace polygon{
+
   template <typename T, typename T2>
   struct PolyLineArbitraryByConcept {};
 
@@ -41,30 +42,45 @@
       return lp(pt, pt1) && lp(pt2, pt);
     }
     
-    template <typename area_type>
-    static inline Unit compute_intercept(const area_type& dy2,
-                                         const area_type& dx1,
-                                         const area_type& dx2) {
+    template <typename T>
+    static inline Unit compute_intercept(const T& dy2,
+                                         const T& dx1,
+                                         const T& dx2) {
       //intercept = dy2 * dx1 / dx2
-      //return (Unit)(((area_type)dy2 * (area_type)dx1) / (area_type)dx2);
-      area_type dx1_q = dx1 / dx2;
-      area_type dx1_r = dx1 % dx2;
-      return dx1_q * dy2 + (dy2 * dx1_r)/dx2;
-    }
+      //return (Unit)(((T)dy2 * (T)dx1) / (T)dx2);
+      T dx1_q = dx1 / dx2;
+      T dx1_r = dx1 % dx2;
+      return static_cast<Unit>(dx1_q * dy2 + (dy2 * dx1_r)/dx2);
+    }
+
+    template <typename area_type, typename integer = gtl_yes>
+    struct equal_slope_i {
+      static inline bool call(area_type dx1, area_type dy1, area_type dx2, area_type dy2) {
+        typedef typename coordinate_traits<Unit>::unsigned_area_type unsigned_product_type;
+        unsigned_product_type cross_1 = (unsigned_product_type)(dx2 < 0 ? -dx2 :dx2) * (unsigned_product_type)(dy1 < 0 ? -dy1 : dy1);
+        unsigned_product_type cross_2 = (unsigned_product_type)(dx1 < 0 ? -dx1 :dx1) * (unsigned_product_type)(dy2 < 0 ? -dy2 : dy2);
+        int dx1_sign = dx1 < 0 ? -1 : 1;
+        int dx2_sign = dx2 < 0 ? -1 : 1;
+        int dy1_sign = dy1 < 0 ? -1 : 1;
+        int dy2_sign = dy2 < 0 ? -1 : 1;
+        int cross_1_sign = dx2_sign * dy1_sign;
+        int cross_2_sign = dx1_sign * dy2_sign;
+        return cross_1 == cross_2 && (cross_1_sign == cross_2_sign || cross_1 == 0);
+      }
+    };
+
+    template <typename area_type>
+    struct equal_slope_i<area_type, gtl_no> {
+      static inline bool call(area_type dx1, area_type dy1, area_type dx2, area_type dy2) {
+        return dx1* dy2 == dx2 * dy1;
+      }
+    };
 
     template <typename area_type>
     static inline bool equal_slope(area_type dx1, area_type dy1, area_type dx2, area_type dy2) {
-      typedef typename coordinate_traits<Unit>::unsigned_area_type unsigned_product_type;
-      unsigned_product_type cross_1 = (unsigned_product_type)(dx2 < 0 ? -dx2 :dx2) * (unsigned_product_type)(dy1 < 0 ? -dy1 : dy1);
-      unsigned_product_type cross_2 = (unsigned_product_type)(dx1 < 0 ? -dx1 :dx1) * (unsigned_product_type)(dy2 < 0 ? -dy2 : dy2);
-      int dx1_sign = dx1 < 0 ? -1 : 1;
-      int dx2_sign = dx2 < 0 ? -1 : 1;
-      int dy1_sign = dy1 < 0 ? -1 : 1;
-      int dy2_sign = dy2 < 0 ? -1 : 1;
-      int cross_1_sign = dx2_sign * dy1_sign;
-      int cross_2_sign = dx1_sign * dy2_sign;
-      return cross_1 == cross_2 && (cross_1_sign == cross_2_sign || cross_1 == 0);
+      return equal_slope_i<area_type, typename coordinate_traits<Unit>::is_integer_type>::call(dx1, dy1, dx2, dy2);
     }
+      
 
     static inline bool equal_slope(const Unit& x, const Unit& y,
                                    const Point& pt1, const Point& pt2) {
@@ -77,36 +93,71 @@
       return equal_slope(dx1, dy1, dx2, dy2);
     }
 
-    template <typename area_type>
-    static inline bool less_slope(area_type dx1, area_type dy1, area_type dx2, area_type dy2) {
-      //reflext x and y slopes to right hand side half plane
-      if(dx1 < 0) {
-        dy1 *= -1;
-        dx1 *= -1;
-      } else if(dx1 == 0) {
-        //if the first slope is vertical the first cannot be less
-        return false;
+    template <typename area_type, typename integer = gtl_yes>
+    struct less_slope_i {
+      static inline bool call(area_type dx1, area_type dy1, area_type dx2, area_type dy2) {
+        //reflect x and y slopes to right hand side half plane
+        if(dx1 < 0) {
+          dy1 *= -1;
+          dx1 *= -1;
+        } else if(dx1 == 0) {
+          //if the first slope is vertical the first cannot be less
+          return false;
+        }
+        if(dx2 < 0) {
+          dy2 *= -1;
+          dx2 *= -1;
+        } else if(dx2 == 0) {
+          //if the second slope is vertical the first is always less unless it is also vertical, in which case they are equal 
+          return dx1 != 0;
+        }
+        typedef typename coordinate_traits<Unit>::unsigned_area_type unsigned_product_type;
+        unsigned_product_type cross_1 = (unsigned_product_type)(dx2 < 0 ? -dx2 :dx2) * (unsigned_product_type)(dy1 < 0 ? -dy1 : dy1);
+        unsigned_product_type cross_2 = (unsigned_product_type)(dx1 < 0 ? -dx1 :dx1) * (unsigned_product_type)(dy2 < 0 ? -dy2 : dy2);
+        int dx1_sign = dx1 < 0 ? -1 : 1;
+        int dx2_sign = dx2 < 0 ? -1 : 1;
+        int dy1_sign = dy1 < 0 ? -1 : 1;
+        int dy2_sign = dy2 < 0 ? -1 : 1;
+        int cross_1_sign = dx2_sign * dy1_sign;
+        int cross_2_sign = dx1_sign * dy2_sign;
+        if(cross_1_sign < cross_2_sign) return true;
+        if(cross_2_sign < cross_1_sign) return false;
+        if(cross_1_sign == -1) return cross_2 < cross_1;
+        return cross_1 < cross_2;
       }
-      if(dx2 < 0) {
-        dy2 *= -1;
-        dx2 *= -1;
-      } else if(dx2 == 0) {
-        //if the second slope is vertical the first is always less unless it is also vertical, in which case they are equal 
-        return dx1 != 0;
+    };
+
+    template <typename area_type>
+    struct less_slope_i<area_type, gtl_no> {
+      static inline bool call(area_type dx1, area_type dy1, area_type dx2, area_type dy2) {
+        //reflect x and y slopes to right hand side half plane
+        if(dx1 < 0) {
+          dy1 *= -1;
+          dx1 *= -1;
+        } else if(dx1 == 0) {
+          //if the first slope is vertical the first cannot be less
+          return false;
+        }
+        if(dx2 < 0) {
+          dy2 *= -1;
+          dx2 *= -1;
+        } else if(dx2 == 0) {
+          //if the second slope is vertical the first is always less unless it is also vertical, in which case they are equal 
+          return dx1 != 0;
+        }
+        area_type cross_1 = dx2 * dy1;
+        area_type cross_2 = dx1 * dy2;
+        if(cross_1 < 0 && cross_2 >= 0) return true;
+        if(cross_2 < 0 && cross_1 >= 0) return false;
+        if(cross_1 < 0) return cross_2 < cross_1;
+        return cross_1 < cross_2;
       }
-      typedef typename coordinate_traits<Unit>::unsigned_area_type unsigned_product_type;
-      unsigned_product_type cross_1 = (unsigned_product_type)(dx2 < 0 ? -dx2 :dx2) * (unsigned_product_type)(dy1 < 0 ? -dy1 : dy1);
-      unsigned_product_type cross_2 = (unsigned_product_type)(dx1 < 0 ? -dx1 :dx1) * (unsigned_product_type)(dy2 < 0 ? -dy2 : dy2);
-      int dx1_sign = dx1 < 0 ? -1 : 1;
-      int dx2_sign = dx2 < 0 ? -1 : 1;
-      int dy1_sign = dy1 < 0 ? -1 : 1;
-      int dy2_sign = dy2 < 0 ? -1 : 1;
-      int cross_1_sign = dx2_sign * dy1_sign;
-      int cross_2_sign = dx1_sign * dy2_sign;
-      if(cross_1_sign < cross_2_sign) return true;
-      if(cross_2_sign < cross_1_sign) return false;
-      if(cross_1_sign == -1) return cross_2 < cross_1;
-      return cross_1 < cross_2;
+    };
+
+    template <typename area_type>
+    static inline bool less_slope(area_type dx1, area_type dy1, area_type dx2, area_type dy2) {
+      return less_slope_i<area_type, 
+        typename coordinate_traits<Unit>::is_integer_type>::call(dx1, dy1, dx2, dy2);
     }
 
     static inline bool less_slope(const Unit& x, const Unit& y,
@@ -337,6 +388,27 @@
       return q + r;
     }
 
+    template <typename T, typename integer = gtl_yes>
+    struct truncate_down_i {
+      static inline Unit call(const T& t) {
+        Unit retval = static_cast<Unit>(t);
+        if(t < static_cast<T>(retval)) --retval;
+        return retval;
+      }
+    };
+
+    template <typename T>
+    struct truncate_down_i<T, gtl_no> {
+      static inline Unit call(const T& t) {
+        return static_cast<Unit>(t);
+      }
+    };
+
+    template <typename T>
+    static inline Unit truncate_down(const T& t) {
+      return truncate_down_i<T, typename coordinate_traits<Unit>::is_integer_type>::call(t);
+    }
+
     static inline bool compute_intersection(Point& intersection, const half_edge& he1, const half_edge& he2) {
       typedef typename high_precision_type<Unit>::type high_precision;
       typedef rectangle_data<Unit> Rectangle;
@@ -394,11 +466,8 @@
       //std::cout << "cross2 " << dy2 << " " << dx1 << " " << dy2 * dx1 << std::endl;
       //Unit exp_x = compute_x_intercept<at>(x11, x21, y11, y21, dy1, dy2, dx1, dx2);
       //Unit exp_y = compute_x_intercept<at>(y11, y21, x11, x21, dx1, dx2, dy1, dy2);
-      Unit x_unit = (Unit)x;
-      Unit y_unit = (Unit)y;
-      //truncate downward if it went up due to negative number
-      if(x < (high_precision)x_unit) --x_unit;
-      if(y < (high_precision)y_unit) --y_unit;
+      Unit x_unit = truncate_down(x);
+      Unit y_unit = truncate_down(y);
       //if(x != exp_x || y != exp_y)
       //  std::cout << exp_x << " " << exp_y << " " << x << " " << y << std::endl;
       //Unit y1 = evalAtXforY(exp_x, he1.first, he1.second);

Modified: sandbox/gtl/boost/polygon/detail/scan_arbitrary.hpp
==============================================================================
--- sandbox/gtl/boost/polygon/detail/scan_arbitrary.hpp	(original)
+++ sandbox/gtl/boost/polygon/detail/scan_arbitrary.hpp	2009-09-03 00:28:56 EDT (Thu, 03 Sep 2009)
@@ -64,59 +64,99 @@
       return lp(pt, pt1) && lp(pt2, pt);
     }
     
+    static inline bool is_unit_integer() {
+      typedef typename coordinate_traits<Unit>::is_integer_type result_type;
+      return result_type::value;
+    }
+
     //quadratic algorithm to do same work as optimal scan for cross checking
     //assume sorted input
     template <typename iT>
     static inline void validate_scan(std::map<segment_id, std::set<Point> >& intersection_points,
                                      iT begin, iT end) {
-      std::set<Point> pts;
-      std::vector<std::pair<half_edge, segment_id> > data(begin, end);
-      for(std::size_t i = 0; i < data.size(); ++i) {
-        if(data[i].first.second < data[i].first.first) {
-          std::swap(data[i].first.first, data[i].first.second);
-        }
-      }
-      std::sort(data.begin(), data.end());
-      //find all intersection points
-      for(typename std::vector<std::pair<half_edge, segment_id> >::iterator outer = data.begin();
-          outer != data.end(); ++outer) {
-        const half_edge& he1 = (*outer).first;
-        //its own end points
-        pts.insert(he1.first);
-        pts.insert(he1.second);
-        for(typename std::vector<std::pair<half_edge, segment_id> >::iterator inner = outer;
-            inner != data.end(); ++inner) {
-          const half_edge& he2 = (*inner).first;
-          if(he1 == he2) continue;
-          if((std::min)(he2. first.get(HORIZONTAL),
-                        he2.second.get(HORIZONTAL)) > 
-             (std::max)(he1.second.get(HORIZONTAL),
-                        he1.first.get(HORIZONTAL)))
-            break;
-          Point intersection;
-          if(compute_intersection(intersection, he1, he2)) {
-            //their intersection point
-            pts.insert(intersection);
-          } 
-        }
-      }
-      //find all segments that interact with intersection points
-      for(typename std::vector<std::pair<half_edge, segment_id> >::iterator outer = data.begin();
-          outer != data.end(); ++outer) {
-        const half_edge& he1 = (*outer).first;
-        segment_id id1 = (*outer).second;
-        typedef rectangle_data<Unit> Rectangle;
-        Rectangle rect1;
-        set_points(rect1, he1.first, he1.second);
-        typename std::set<Point>::iterator itr = pts.lower_bound((std::min)(he1.first, he1.second));
-        typename std::set<Point>::iterator itr2 = pts.upper_bound((std::max)(he1.first, he1.second));
-        while(itr != pts.end() && itr != pts.begin() && (*itr).get(HORIZONTAL) >= (std::min)(he1.first.get(HORIZONTAL), he1.second.get(HORIZONTAL))) --itr;
-        while(itr2 != pts.end() && (*itr2).get(HORIZONTAL) <= (std::max)(he1.first.get(HORIZONTAL), he1.second.get(HORIZONTAL))) ++itr2;
-        //itr = pts.begin();
-        //itr2 = pts.end();
-        for( ; itr != itr2; ++itr) {
-          if(intersects_grid(*itr, he1))
-            intersection_points[id1].insert(*itr);
+      if(is_unit_integer()) {
+        std::set<Point> pts;
+        std::vector<std::pair<half_edge, segment_id> > data(begin, end);
+        for(std::size_t i = 0; i < data.size(); ++i) {
+          if(data[i].first.second < data[i].first.first) {
+            std::swap(data[i].first.first, data[i].first.second);
+          }
+        }
+        std::sort(data.begin(), data.end());
+        //find all intersection points
+        for(typename std::vector<std::pair<half_edge, segment_id> >::iterator outer = data.begin();
+            outer != data.end(); ++outer) {
+          const half_edge& he1 = (*outer).first;
+          //its own end points
+          pts.insert(he1.first);
+          pts.insert(he1.second);
+          for(typename std::vector<std::pair<half_edge, segment_id> >::iterator inner = outer;
+              inner != data.end(); ++inner) {
+            const half_edge& he2 = (*inner).first;
+            if(he1 == he2) continue;
+            if((std::min)(he2. first.get(HORIZONTAL),
+                          he2.second.get(HORIZONTAL)) > 
+               (std::max)(he1.second.get(HORIZONTAL),
+                          he1.first.get(HORIZONTAL)))
+              break;
+            Point intersection;
+            if(compute_intersection(intersection, he1, he2)) {
+              //their intersection point
+              pts.insert(intersection);
+            } 
+          }
+        }
+        //find all segments that interact with intersection points
+        for(typename std::vector<std::pair<half_edge, segment_id> >::iterator outer = data.begin();
+            outer != data.end(); ++outer) {
+          const half_edge& he1 = (*outer).first;
+          segment_id id1 = (*outer).second;
+          typedef rectangle_data<Unit> Rectangle;
+          Rectangle rect1;
+          set_points(rect1, he1.first, he1.second);
+          typename std::set<Point>::iterator itr = pts.lower_bound((std::min)(he1.first, he1.second));
+          typename std::set<Point>::iterator itr2 = pts.upper_bound((std::max)(he1.first, he1.second));
+          while(itr != pts.end() && itr != pts.begin() && (*itr).get(HORIZONTAL) >= (std::min)(he1.first.get(HORIZONTAL), he1.second.get(HORIZONTAL))) --itr;
+          while(itr2 != pts.end() && (*itr2).get(HORIZONTAL) <= (std::max)(he1.first.get(HORIZONTAL), he1.second.get(HORIZONTAL))) ++itr2;
+          //itr = pts.begin();
+          //itr2 = pts.end();
+          for( ; itr != itr2; ++itr) {
+            if(intersects_grid(*itr, he1))
+              intersection_points[id1].insert(*itr);
+          }
+        }
+      } else {
+        //no snap rounding in floating point case
+        std::vector<std::pair<half_edge, segment_id> > data(begin, end);
+        for(std::size_t i = 0; i < data.size(); ++i) {
+          if(data[i].first.second < data[i].first.first) {
+            std::swap(data[i].first.first, data[i].first.second);
+          }
+        }
+        std::sort(data.begin(), data.end());
+        //find all intersection points
+        for(typename std::vector<std::pair<half_edge, segment_id> >::iterator outer = data.begin();
+            outer != data.end(); ++outer) {
+          const half_edge& he1 = (*outer).first;
+          //its own end points
+          intersection_points[(*outer).second].insert(he1.first);
+          intersection_points[(*outer).second].insert(he1.second);
+          for(typename std::vector<std::pair<half_edge, segment_id> >::iterator inner = outer;
+              inner != data.end(); ++inner) {
+            const half_edge& he2 = (*inner).first;
+            if(he1 == he2) continue;
+            if((std::min)(he2. first.get(HORIZONTAL),
+                          he2.second.get(HORIZONTAL)) > 
+               (std::max)(he1.second.get(HORIZONTAL),
+                          he1.first.get(HORIZONTAL)))
+              break;
+            Point intersection;
+            if(compute_intersection(intersection, he1, he2)) {
+              //their intersection point
+              intersection_points[(*outer).second].insert(intersection);
+              intersection_points[(*inner).second].insert(intersection);
+            } 
+          }
         }
       }
     }