[Boost-commit] svn:boost r76845 - in sandbox/gtl: boost/polygon libs/polygon/test libs/polygon/voronoi_example

Subject: [Boost-commit] svn:boost r76845 - in sandbox/gtl: boost/polygon libs/polygon/test libs/polygon/voronoi_example
From: sydorchuk.andriy_at_[hidden]
Date: 2012-02-02 18:09:50

Author: asydorchuk
Date: 2012-02-02 18:09:48 EST (Thu, 02 Feb 2012)
New Revision: 76845
URL: http://svn.boost.org/trac/boost/changeset/76845

Log:
Added voronoi_diagram_traits.
Moved voronoi_utils to a separate file unit.
Tests update.
Added:
   sandbox/gtl/boost/polygon/voronoi_utils.hpp   (contents, props changed)
Text files modified: 
   sandbox/gtl/boost/polygon/voronoi_diagram.hpp                   |   610 +++++++-------------------------------- 
   sandbox/gtl/libs/polygon/test/voronoi_benchmark_test.cpp        |     1                                         
   sandbox/gtl/libs/polygon/test/voronoi_builder_test.cpp          |     1                                         
   sandbox/gtl/libs/polygon/test/voronoi_clipping_test.cpp         |    50 +-                                      
   sandbox/gtl/libs/polygon/test/voronoi_test_helper.hpp           |    31 +                                       
   sandbox/gtl/libs/polygon/voronoi_example/voronoi_visualizer.cpp |     9                                         
   6 files changed, 168 insertions(+), 534 deletions(-)

Modified: sandbox/gtl/boost/polygon/voronoi_diagram.hpp
==============================================================================
--- sandbox/gtl/boost/polygon/voronoi_diagram.hpp	(original)
+++ sandbox/gtl/boost/polygon/voronoi_diagram.hpp	2012-02-02 18:09:48 EST (Thu, 02 Feb 2012)
@@ -10,13 +10,11 @@
 #ifndef BOOST_POLYGON_VORONOI_DIAGRAM
 #define BOOST_POLYGON_VORONOI_DIAGRAM
 
-#include <cmath>
 #include <list>
-#include <stack>
 #include <vector>
 
-#include "polygon.hpp"
 #include "detail/voronoi_ctypes.hpp"
+#include "detail/voronoi_structures.hpp"
 
 namespace boost {
 namespace polygon {
@@ -31,64 +29,64 @@
     // Bounding rectangle data structure. Contains coordinates
     // of the bottom left and the upper right points of the rectangle.
     template <typename T>
-    class BRect {
+    class bounding_rectangle {
     public:
         typedef T coordinate_type;
-        typedef point_data<coordinate_type> point_type;
 
-        BRect() {}
+        bounding_rectangle() {}
+
+        bounding_rectangle(T x, T y) {
+            x_min_ = x_max_ = x;
+            y_min_ = y_max_ = y;
+        }
 
         template <typename P>
-        BRect(const P &p) {
-            x_min_ = x_max_ = static_cast<coordinate_type>(p.x());
-            y_min_ = y_max_ = static_cast<coordinate_type>(p.y());
+        bounding_rectangle(const P &p) {
+            x_min_ = x_max_ = p.x();
+            y_min_ = y_max_ = p.y();
+        }
+
+        bounding_rectangle(T x1, T y1, T x2, T y2) {
+            x_min_ = (std::min)(x1, x2);
+            y_min_ = (std::min)(y1, y2);
+            x_max_ = (std::max)(x1, x2);
+            y_max_ = (std::max)(y1, y2);
         }
 
         template <typename P>
-        BRect(const P &p1, const P &p2) {
-            x_min_ = (std::min)(static_cast<coordinate_type>(p1.x()),
-                                static_cast<coordinate_type>(p2.x()));
-            y_min_ = (std::min)(static_cast<coordinate_type>(p1.y()),
-                                static_cast<coordinate_type>(p2.y()));
-            x_max_ = (std::max)(static_cast<coordinate_type>(p1.x()),
-                                static_cast<coordinate_type>(p2.x()));
-            y_max_ = (std::max)(static_cast<coordinate_type>(p1.y()),
-                                static_cast<coordinate_type>(p2.y()));
-        }
-
-        template <typename C>
-        BRect(C x_mn, C y_mn, C x_mx, C y_mx) {
-            x_min_ = (std::min)(static_cast<coordinate_type>(x_mn),
-                                static_cast<coordinate_type>(x_mx));
-            y_min_ = (std::min)(static_cast<coordinate_type>(y_mn),
-                                static_cast<coordinate_type>(y_mx));
-            x_max_ = (std::max)(static_cast<coordinate_type>(x_mn),
-                                static_cast<coordinate_type>(x_mx));
-            y_max_ = (std::max)(static_cast<coordinate_type>(y_mn),
-                                static_cast<coordinate_type>(y_mx));
+        bounding_rectangle(const P &p1, const P &p2) {
+            x_min_ = (std::min)(p1.x(), p2.x());
+            y_min_ = (std::min)(p1.y(), p2.y());
+            x_max_ = (std::max)(p1.x(), p2.x());
+            y_max_ = (std::max)(p1.y(), p2.y());
+        }
+
+        void update(T x, T y) {
+            x_min_ = (std::min)(x_min_, x);
+            y_min_ = (std::min)(y_min_, y);
+            x_max_ = (std::max)(x_max_, x);
+            y_max_ = (std::max)(y_max_, y);
         }
 
         // Extend the rectangle with a new point.
         template <typename P>
         void update(const P &p) {
-            x_min_ = (std::min)(x_min_, static_cast<coordinate_type>(p.x()));
-            y_min_ = (std::min)(y_min_, static_cast<coordinate_type>(p.y()));
-            x_max_ = (std::max)(x_max_, static_cast<coordinate_type>(p.x()));
-            y_max_ = (std::max)(y_max_, static_cast<coordinate_type>(p.y()));
+            x_min_ = (std::min)(x_min_, p.x());
+            y_min_ = (std::min)(y_min_, p.y());
+            x_max_ = (std::max)(x_max_, p.x());
+            y_max_ = (std::max)(y_max_, p.y());
+        }
+
+        bool contains(T x, T y) const {
+            return x >= x_min_ && x <= x_max_ &&
+                   y >= y_min_ && y <= y_max_;
         }
 
         // Check whether a point is situated inside the bounding rectangle.
         template <typename P>
         bool contains(const P &p) const {
-            return static_cast<coordinate_type>(p.x()) >= x_min_ &&
-                   static_cast<coordinate_type>(p.x()) <= x_max_ &&
-                   static_cast<coordinate_type>(p.y()) >= y_min_ &&
-                   static_cast<coordinate_type>(p.y()) <= y_max_;
-        }
-
-        // Check whether the bounding rectangle has a non-zero area.
-        bool verify() const {
-            return (x_min_ <= x_max_) && (y_min_ <= y_max_);
+            return p.x() >= x_min_ && p.x() <= x_max_ &&
+                   p.y() >= y_min_ && p.y() <= y_max_;
         }
 
         // Return the x-coordinate of the bottom left point of the rectangle.
@@ -126,412 +124,6 @@
         coordinate_type y_max_;
     };
 
-    // Voronoi output postprocessing tools.
-    template <typename fpt_>
-    class voronoi_helper {
-    public:
-        typedef fpt_ fpt_type;
-        typedef point_data<fpt_type> point_type;
-        typedef std::vector<point_type> point_set_type;
-        typedef directed_line_segment_data<fpt_type> segment_type;
-        typedef directed_line_segment_set_data<fpt_type> segment_set_type;
-        typedef BRect<fpt_type> brect_type;
-
-        // There are three different types of edges:
-        //   1) Segment edge - has both endpoints;
-        //   2) Ray edge - has one endpoint, infinite
-        //                 in the positive direction;
-        //   3) Line edge - is infinite in both directions.
-        enum kEdgeType {
-            SEGMENT = 0,
-            RAY = 1,
-            LINE = 2
-        };
-
-        // Get a view rectangle based on the voronoi bounding rectangle.
-        template <typename T>
-        static brect_type get_view_rectangle(const BRect<T> &brect) {
-            fpt_type x_len = static_cast<fpt_type>(brect.x_max()) -
-                             static_cast<fpt_type>(brect.x_min());
-            fpt_type y_len = static_cast<fpt_type>(brect.y_max()) -
-                             static_cast<fpt_type>(brect.y_min());
-            fpt_type x_mid = static_cast<fpt_type>(brect.x_max()) +
-                             static_cast<fpt_type>(brect.x_min());
-            fpt_type y_mid = static_cast<fpt_type>(brect.y_max()) +
-                             static_cast<fpt_type>(brect.y_min());
-            x_mid /= 2;
-            y_mid /= 2;
-            fpt_type offset = x_len;
-            if (offset < y_len)
-                offset = y_len;
-            if (offset == static_cast<fpt_type>(0))
-                offset = static_cast<fpt_type>(1);
-            brect_type new_brect(x_mid - offset, y_mid - offset,
-                                 x_mid + offset, y_mid + offset);
-            return new_brect;
-        }
-
-        // Compute intermediate points for the voronoi edge within the given
-        // bounding rectangle. The assumption is made that voronoi rectangle
-        // contains all the finite endpoints of the edge.
-        // Max_error is the maximum distance (relative to the minor of two
-        // rectangle sides) allowed between the parabola and line segments
-        // that interpolate it.
-        template <typename T>
-        static point_set_type get_point_interpolation(
-            const voronoi_edge<T> *edge,
-            const BRect<T> &brect,
-            fpt_type max_error) {
-            // Retrieve the cell to the left of the current edge.
-            const voronoi_cell<T> *cell1 = edge->cell();
-
-            // Retrieve the cell to the right of the current edge.
-            const voronoi_cell<T> *cell2 = edge->twin()->cell();
-
-            // edge_points - contains intermediate points.
-            point_set_type edge_points;
-            if (!(cell1->contains_segment() ^ cell2->contains_segment())) {
-                // Edge is a segment, ray, or line.
-                if (edge->vertex0() != NULL && edge->vertex1() != NULL) {
-                    // Edge is a segment. No further processing is required.
-                    edge_points.push_back(
-                        get_point(edge->vertex0()->vertex()));
-                    edge_points.push_back(
-                        get_point(edge->vertex1()->vertex()));
-                } else {
-                    // Edge is a ray or line.
-                    // Clip it with the bounding rectangle.
-                    const point_type &site1 = get_point(cell1->point0());
-                    const point_type &site2 = get_point(cell2->point0());
-                    point_type origin((site1.x() + site2.x()) / 2,
-                                      (site1.y() + site2.y()) / 2);
-                    point_type direction(site1.y() - site2.y(),
-                                         site2.x() - site1.x());
-
-                    // Find intersection points.
-                    find_intersections(origin, direction, LINE,
-                                       brect, edge_points);
-
-                    // Update endpoints in case edge is a ray.
-                    if (edge->vertex1() != NULL)
-                        edge_points[1] = get_point(edge->vertex1()->vertex());
-                    if (edge->vertex0() != NULL)
-                        edge_points[0] = get_point(edge->vertex0()->vertex());
-                }
-            } else {
-                // point1 - site point;
-                const point_type &point1 = (cell1->contains_segment()) ?
-                    get_point(cell2->point0()) : get_point(cell1->point0());
-
-                // point2 - startpoint of the segment site;
-                const point_type &point2 = (cell1->contains_segment()) ?
-                    get_point(cell1->point0()) : get_point(cell2->point0());
-
-                // point3 - endpoint of the segment site;
-                const point_type &point3 = (cell1->contains_segment()) ?
-                    get_point(cell1->point1()) : get_point(cell2->point1());
-
-                if (edge->vertex0() != NULL && edge->vertex1() != NULL) {
-                    // Edge is a segment or parabolic arc.
-                    edge_points.push_back(
-                        get_point(edge->vertex0()->vertex()));
-                    edge_points.push_back(
-                        get_point(edge->vertex1()->vertex()));
-                    fpt_type max_dist = max_error * brect.min_len();
-                    fill_intermediate_points(point1, point2, point3,
-                                             edge_points, max_dist);
-                } else {
-                    // Edge is a ray or line.
-                    // Clip it with the bounding rectangle.
-                    fpt_type dir_x =
-                        (cell1->contains_segment() ^ (point1 == point3)) ?
-                        point3.y() - point2.y() : point2.y() - point3.y();
-                    fpt_type dir_y =
-                        (cell1->contains_segment() ^ (point1 == point3)) ?
-                        point2.x() - point3.x() : point3.x() - point2.x();
-                    point_type direction(dir_x, dir_y);
-
-                    // Find intersection points.
-                    find_intersections(point1, direction, LINE,
-                                       brect, edge_points);
-
-                    // Update endpoints in case edge is a ray.
-                    if (edge->vertex1() != NULL)
-                        edge_points[1] = get_point(edge->vertex1()->vertex());
-                    if (edge->vertex0() != NULL)
-                        edge_points[0] = get_point(edge->vertex0()->vertex());
-                }
-            }
-            return edge_points;
-        }
-
-        // Interpolate voronoi edge with a set of segments to satisfy maximal
-        // error requirement.
-        template <typename T>
-        static segment_set_type get_segment_interpolation(
-            const voronoi_edge<T> *edge,
-            const BRect<T> &brect,
-            fpt_type max_error) {
-            point_set_type point_interpolation =
-                get_point_interpolcation(edge, brect, max_error);
-            segment_set_type ret_val;
-            for (size_t i = 1; i < point_interpolation.size(); ++i)
-                ret_val.insert(segment_type(point_interpolation[i-1],
-                                            point_interpolation[i]));
-            return ret_val;
-        }
-
-        // Find edge-rectangle intersection points.
-        // Edge is represented by its origin, direction and type.
-        static void find_intersections(
-                const point_type &origin, const point_type &direction,
-                kEdgeType edge_type, const brect_type &brect,
-                point_set_type &intersections) {
-            // Find the center of the rectangle.
-            fpt_type center_x = (brect.x_min() + brect.x_max()) / 2;
-            fpt_type center_y = (brect.y_min() + brect.y_max()) / 2;
-
-            // Find the half-diagonal vector of the rectangle.
-            fpt_type len_x = brect.x_max() - center_x;
-            fpt_type len_y = brect.y_max() - center_y;
-
-            // Find the vector between the origin and the center of the
-            // rectangle.
-            fpt_type diff_x = origin.x() - center_x;
-            fpt_type diff_y = origin.y() - center_y;
-
-            // Find the vector perpendicular to the direction vector.
-            fpt_type perp_x = direction.y();
-            fpt_type perp_y = -direction.x();
-
-            // Projection of the vector between the origin and the center of
-            // the rectangle on the line perpendicular to the direction vector.
-            fpt_type lexpr = magnitude(perp_x * diff_x + perp_y * diff_y);
-
-            // Maximum projection among any of the half-diagonals of the
-            // rectangle on the line perpendicular to the direction vector.
-            fpt_type rexpr = magnitude(perp_x * len_x) +
-                             magnitude(perp_y * len_y);
-
-            // Intersection check. Compare projections.
-            if (lexpr > rexpr)
-                return;
-
-            // Intersection parameters. If fT[i]_used is true than:
-            // origin + fT[i] * direction = intersection point, i = 0 .. 1.
-            // For different edge types next fT values are acceptable:
-            // Segment: [0; 1];
-            // Ray: [0; infinity];
-            // Line: [-infinity; infinity].
-            bool fT0_used = false;
-            bool fT1_used = false;
-            fpt_type fT0 = 0;
-            fpt_type fT1 = 0;
-
-            // Check for the intersections with the lines
-            // going through the sides of the bounding rectangle.
-            clip_line(+direction.x(), -diff_x - len_x,
-                      fT0_used, fT1_used, fT0, fT1);
-            clip_line(-direction.x(), +diff_x - len_x,
-                      fT0_used, fT1_used, fT0, fT1);
-            clip_line(+direction.y(), -diff_y - len_y,
-                      fT0_used, fT1_used, fT0, fT1);
-            clip_line(-direction.y(), +diff_y - len_y,
-                      fT0_used, fT1_used, fT0, fT1);
-
-            // Update intersections vector.
-            if (fT0_used && check_extent(fT0, edge_type))
-                intersections.push_back(point_type(
-                    origin.x() + fT0 * direction.x(),
-                    origin.y() + fT0 * direction.y()));
-            if (fT1_used && fT0 != fT1 && check_extent(fT1, edge_type))
-                intersections.push_back(point_type(
-                    origin.x() + fT1 * direction.x(),
-                    origin.y() + fT1 * direction.y()));
-        }
-
-    private:
-        voronoi_helper();
-
-        template <typename P>
-        static point_type get_point(const P &point) {
-            fpt_type x = static_cast<fpt_type>(point.x());
-            fpt_type y = static_cast<fpt_type>(point.y());
-            return point_type(x, y);
-        }
-
-        // Find intermediate points of the parabola. Number of points
-        // is defined by the value of max_dist parameter - maximum distance
-        // between parabola and line segments that interpolate it.
-        // Parabola is a locus of points equidistant from the point and segment
-        // sites. intermediate_points should contain two initial endpoints
-        // of the edge (voronoi vertices). Intermediate points are inserted
-        // between the given two endpoints.
-        // Max_dist is the maximum distance allowed between parabola and line
-        // segments that interpolate it.
-        static void fill_intermediate_points(
-                point_type point_site,
-                point_type segment_site_start,
-                point_type segment_site_end,
-                point_set_type &intermediate_points,
-                fpt_type max_dist) {
-            // Check if bisector is a segment.
-            if (point_site == segment_site_start ||
-                point_site == segment_site_end)
-                return;
-
-            // Apply the linear transformation to move start point of the
-            // segment to the point with coordinates (0, 0) and the direction
-            // of the segment to coincide the positive direction of the x-axis.
-            fpt_type segm_vec_x = segment_site_end.x() -
-                                  segment_site_start.x();
-            fpt_type segm_vec_y = segment_site_end.y() -
-                                  segment_site_start.y();
-            fpt_type sqr_segment_length = segm_vec_x * segm_vec_x +
-                                          segm_vec_y * segm_vec_y;
-
-            // Compute x-coordinates of the endpoints of the edge
-            // in the transformed space.
-            fpt_type projection_start = sqr_segment_length *
-                get_point_projection(intermediate_points[0],
-                                     segment_site_start,
-                                     segment_site_end);
-            fpt_type projection_end = sqr_segment_length *
-                get_point_projection(intermediate_points[1],
-                                     segment_site_start,
-                                     segment_site_end);
-
-            // Compute parabola parameterers in the transformed space.
-            // Parabola has next representation:
-            // f(x) = ((x-rot_x)^2 + rot_y^2) / (2.0*rot_y).
-            fpt_type point_vec_x = point_site.x() - segment_site_start.x();
-            fpt_type point_vec_y = point_site.y() - segment_site_start.y();
-            fpt_type rot_x = segm_vec_x * point_vec_x +
-                             segm_vec_y * point_vec_y;
-            fpt_type rot_y = segm_vec_x * point_vec_y -
-                             segm_vec_y * point_vec_x;
-
-            // Save the last point.
-            point_type last_point = intermediate_points[1];
-            intermediate_points.pop_back();
-
-            // Use stack to avoid recursion.
-            std::stack<fpt_type> point_stack;
-            point_stack.push(projection_end);
-            fpt_type cur_x = projection_start;
-            fpt_type cur_y = parabola_y(cur_x, rot_x, rot_y);
-
-            // Adjust max_dist parameter in the transformed space.
-            max_dist *= max_dist * sqr_segment_length;
-
-            while (!point_stack.empty()) {
-                fpt_type new_x = point_stack.top();
-                fpt_type new_y = parabola_y(new_x, rot_x, rot_y);
-
-                // Compute coordinates of the point of the parabola that is
-                // furthest from the current line segment.
-                fpt_type mid_x = (new_y - cur_y) / (new_x - cur_x) * rot_y +
-                                 rot_x;
-                fpt_type mid_y = parabola_y(mid_x, rot_x, rot_y);
-
-                // Compute maximum distance between the given parabolic arc
-                // and line segment that interpolates it.
-                fpt_type dist = (new_y - cur_y) * (mid_x - cur_x) -
-                                (new_x - cur_x) * (mid_y - cur_y);
-                dist = dist * dist / ((new_y - cur_y) * (new_y - cur_y) +
-                                      (new_x - cur_x) * (new_x - cur_x));
-                if (dist <= max_dist) {
-                    // Distance between parabola and line segment is
-                    // not greater than max_dist.
-                    point_stack.pop();
-                    fpt_type inter_x =
-                        (segm_vec_x * new_x - segm_vec_y * new_y) /
-                        sqr_segment_length + segment_site_start.x();
-                    fpt_type inter_y =
-                        (segm_vec_x * new_y + segm_vec_y * new_x) /
-                        sqr_segment_length + segment_site_start.y();
-                    intermediate_points.push_back(
-                        point_type(inter_x, inter_y));
-                    cur_x = new_x;
-                    cur_y = new_y;
-                } else {
-                    point_stack.push(mid_x);
-                }
-            }
-
-            // Update last point.
-            intermediate_points.back() = last_point;
-        }
-
-        // Compute y(x) = ((x - a) * (x - a) + b * b) / (2 * b).
-        static fpt_type parabola_y(fpt_type x, fpt_type a, fpt_type b) {
-            return ((x - a) * (x - a) + b * b) / (2.0 * b);
-        }
-
-        // Check whether extent is compatible with the edge type.
-        static bool check_extent(fpt_type extent, kEdgeType etype) {
-            switch (etype) {
-                case SEGMENT:
-                    return extent >= static_cast<fpt_type>(0) &&
-                           extent <= static_cast<fpt_type>(1);
-                case RAY: return extent >= static_cast<fpt_type>(0);
-                case LINE: return true;
-            }
-            return true;
-        }
-
-        // Compute the absolute value.
-        static inline fpt_type magnitude(fpt_type value) {
-            if (value >= static_cast<fpt_type>(0))
-                return value;
-            return -value;
-        }
-
-        // Find fT min and max values: fT = numer / denom.
-        static void clip_line(fpt_type denom, fpt_type numer,
-                              bool &fT0_used, bool &fT1_used,
-                              fpt_type &fT0, fpt_type &fT1) {
-            if (denom > static_cast<fpt_type>(0)) {
-                if (fT1_used && numer > denom * fT1)
-                    return;
-                if (!fT0_used || numer > denom * fT0) {
-                    fT0_used = true;
-                    fT0 = numer / denom;
-                }
-            } else if (denom < static_cast<fpt_type>(0)) {
-                if (fT0_used && numer > denom * fT0)
-                    return;
-                if (!fT1_used || numer > denom * fT1) {
-                    fT1_used = true;
-                    fT1 = numer / denom;
-                }
-            }
-        }
-
-        // Get normalized length of the distance between:
-        //     1) point projection onto the segment;
-        //     2) start point of the segment.
-        // Return this length divided by the segment length.
-        // This is made to avoid sqrt computation during transformation from
-        // the initial space to the transformed one and vice versa.
-        // Assumption is made that projection of the point lies
-        // between the startpoint and endpoint of the segment.
-        static fpt_type get_point_projection(
-                const point_type &point,
-                const point_type &segment_start,
-                const point_type &segment_end) {
-            fpt_type segment_vec_x = segment_end.x() - segment_start.x();
-            fpt_type segment_vec_y = segment_end.y() - segment_start.y();
-            fpt_type point_vec_x = point.x() - segment_start.x();
-            fpt_type point_vec_y = point.y() - segment_start.y();
-            fpt_type sqr_segment_length = segment_vec_x * segment_vec_x +
-                                          segment_vec_y * segment_vec_y;
-            fpt_type vec_dot = segment_vec_x * point_vec_x +
-                               segment_vec_y * point_vec_y;
-            return vec_dot / sqr_segment_length;
-        }
-    };
-
     // Represents voronoi cell.
     // Data members: 1) pointer to the incident edge;
     //               2) site inside the cell;
@@ -541,15 +133,20 @@
     class voronoi_cell {
     public:
         typedef T coordinate_type;
-        typedef point_data<coordinate_type> point_type;
+        typedef detail::point_2d<coordinate_type> point_type;
         typedef voronoi_edge<coordinate_type> voronoi_edge_type;
 
-        template <typename P>
-        voronoi_cell(const P &p1, const P &p2, voronoi_edge_type *edge) :
-            point0_(static_cast<coordinate_type>(p1.x()),
-                    static_cast<coordinate_type>(p1.y())),
-            point1_(static_cast<coordinate_type>(p2.x()),
-                    static_cast<coordinate_type>(p2.y())),
+        voronoi_cell(const point_type &p1, voronoi_edge_type *edge) :
+            point0_(p1),
+            point1_(p1),
+            incident_edge_(edge),
+            num_incident_edges_(0) {}
+
+        voronoi_cell(const point_type &p1,
+                     const point_type &p2,
+                     voronoi_edge_type *edge) :
+            point0_(p1),
+            point1_(p2),
             incident_edge_(edge),
             num_incident_edges_(0) {}
 
@@ -589,13 +186,12 @@
     class voronoi_vertex {
     public:
         typedef T coordinate_type;
-        typedef point_data<T> point_type;
+        typedef detail::point_2d<T> point_type;
         typedef voronoi_edge<coordinate_type> voronoi_edge_type;
 
-        template <typename P>
-        voronoi_vertex(const P &vertex, voronoi_edge_type *edge) :
-            vertex_(static_cast<coordinate_type>(vertex.x()),
-                    static_cast<coordinate_type>(vertex.y())),
+        voronoi_vertex(const point_type &vertex,
+                       voronoi_edge_type *edge) :
+            vertex_(vertex),
             incident_edge_(edge),
             num_incident_edges_(3) {}
 
@@ -724,6 +320,22 @@
         voronoi_edge_type *prev_;
     };
 
+    template <typename T>
+    struct voronoi_diagram_traits {
+        typedef T coordinate_type;
+        typedef struct {
+            template <typename CT>
+            coordinate_type operator()(const CT& value) {
+                return static_cast<coordinate_type>(value);
+            }
+        } ctype_converter_type;
+        typedef detail::point_2d<coordinate_type> point_type;
+        typedef bounding_rectangle<coordinate_type> bounding_rectangle_type;
+        typedef voronoi_cell<coordinate_type> voronoi_cell_type;
+        typedef voronoi_vertex<coordinate_type> voronoi_vertex_type;
+        typedef voronoi_edge<coordinate_type> voronoi_edge_type;
+    };
+
     // Voronoi output data structure.
     // Data members:
     //   1) cell_records_ - vector of the voronoi cells;
@@ -750,27 +362,29 @@
     //     4) implement simplification via copying not degenerate elements
     //        to the new vector as removing elements from a vector takes O(n)
     //        time.
-    template <typename T>
+    template <typename T, typename TRAITS = voronoi_diagram_traits<T> >
     class voronoi_diagram {
     public:
-        typedef T coordinate_type;
-        typedef point_data<coordinate_type> point_type;
+        typedef typename TRAITS::coordinate_type coordinate_type;
+        typedef typename TRAITS::ctype_converter_type ctype_converter_type;
+        typedef typename TRAITS::point_type point_type;
+        typedef typename TRAITS::bounding_rectangle_type bounding_rectangle_type;
+        typedef typename TRAITS::voronoi_cell_type voronoi_cell_type;
+        typedef typename TRAITS::voronoi_vertex_type voronoi_vertex_type;
+        typedef typename TRAITS::voronoi_edge_type voronoi_edge_type;
 
-        typedef voronoi_cell<coordinate_type> voronoi_cell_type;
         typedef std::vector<voronoi_cell_type> voronoi_cells_type;
         typedef typename voronoi_cells_type::iterator
             voronoi_cell_iterator_type;
         typedef typename voronoi_cells_type::const_iterator
             voronoi_cell_const_iterator_type;
 
-        typedef voronoi_vertex<coordinate_type> voronoi_vertex_type;
         typedef std::list<voronoi_vertex_type> voronoi_vertices_type;
         typedef typename voronoi_vertices_type::iterator
             voronoi_vertex_iterator_type;
         typedef typename voronoi_vertices_type::const_iterator
             voronoi_vertex_const_iterator_type;
 
-        typedef voronoi_edge<coordinate_type> voronoi_edge_type;
         typedef std::list<voronoi_edge_type> voronoi_edges_type;
         typedef typename voronoi_edges_type::iterator
             voronoi_edge_iterator_type;
@@ -792,7 +406,7 @@
             num_vertex_records_ = 0;
         }
 
-        const BRect<coordinate_type> &bounding_rectangle() const {
+        const bounding_rectangle_type &bounding_rectangle() const {
             return voronoi_rect_;
         }
 
@@ -834,12 +448,11 @@
         template <typename SEvent>
         void process_single_site(const SEvent &site) {
             // Update bounding rectangle.
-            voronoi_rect_ = BRect<coordinate_type>(site.point0());
+            point_type p = prepare_point(site.point0());
+            voronoi_rect_ = bounding_rectangle_type(p);
 
             // Update cell records.
-            cell_records_.push_back(voronoi_cell_type(site.point0(),
-                                                      site.point1(),
-                                                      NULL));
+            cell_records_.push_back(voronoi_cell_type(p, NULL));
         }
 
         // Insert a new half-edge into the output data structure.
@@ -847,7 +460,7 @@
         // Returns a pointer to a new half-edge.
         template <typename SEvent>
         std::pair<void*, void*> insert_new_edge(const SEvent &site1,
-                              const SEvent &site2) {
+                                                const SEvent &site2) {
             // Get sites' indices.
             int site_index1 = site1.index();
             int site_index2 = site2.index();
@@ -862,25 +475,20 @@
 
             // Add the initial cell during the first edge insertion.
             if (cell_records_.empty()) {
-                cell_records_.push_back(voronoi_cell_type(site1.point0(),
-                                                          site1.point1(),
-                                                          &edge1));
-                voronoi_rect_ = BRect<coordinate_type>(site1.point0(),
-                                                       site1.point0());
+                process_single_site(site1);
+                cell_records_.back().incident_edge(&edge1);
             }
             cell_records_[site_index1].inc_num_incident_edges();
 
             // Update the bounding rectangle.
-            voronoi_rect_.update(site2.point0());
-            if (site2.is_segment()) {
-                voronoi_rect_.update(site2.point1());	
-            }
+            voronoi_rect_.update(prepare_point(site2.point0()));
 
             // The second site represents a new site during site event
             // processing. Add a new cell to the cell records.
-            cell_records_.push_back(voronoi_cell_type(site2.point0(),
-                                                      site2.point1(),
-                                                      &edge2));
+            cell_records_.push_back(voronoi_cell_type(
+                prepare_point(site2.point0()),
+                prepare_point(site2.point1()),
+                &edge2));
             cell_records_.back().inc_num_incident_edges();
 
             // Set up pointers to cells.
@@ -911,8 +519,10 @@
             voronoi_edge_type *edge23 = static_cast<voronoi_edge_type*>(data23);
 
             // Add a new voronoi vertex.
-            vertex_records_.push_back(voronoi_vertex_type(
-                point_type(circle.x(), circle.y()), edge12));
+            coordinate_type x = convert_(circle.x());
+            coordinate_type y = convert_(circle.y());
+            vertex_records_.push_back(
+                voronoi_vertex_type(point_type(x, y), edge12));
             voronoi_vertex_type &new_vertex = vertex_records_.back();
 
             // Update vertex pointers of the old edges.
@@ -1010,9 +620,9 @@
                 const voronoi_vertex_type *v2 = edge_it1->vertex1();
 
                 // Make epsilon robust check.
-                if (ulp_cmp(v1->vertex().x(), v2->vertex().x(), 256) ==
+                if (ulp_cmp_(v1->vertex().x(), v2->vertex().x(), 128) ==
                     detail::ulp_comparison<T>::EQUAL &&
-                    ulp_cmp(v1->vertex().y(), v2->vertex().y(), 256) ==
+                    ulp_cmp_(v1->vertex().y(), v2->vertex().y(), 128) ==
                     detail::ulp_comparison<T>::EQUAL) {
                     // Decrease number of cell's incident edges.
                     edge_it1->cell()->dec_num_incident_edges();
@@ -1044,8 +654,9 @@
                         }
                     }
 
-                    // To guarantee N*lon(N) time we merge vertex with
-                    // less incident edges to the one with more.
+                    // To guarantee O(N) time for all removal operations we
+                    // merge vertex with less incident edges to the one with
+                    // more.
                     if (v1->num_incident_edges() >= v2->num_incident_edges()) {
                         remove_edge(&(*edge_it1));
                     } else {
@@ -1105,6 +716,18 @@
         }
 
     private:
+        template <typename T>
+        point_type prepare_point(const T& x, const T& y) {
+            coordinate_type xc = convert_(x);
+            coordinate_type yc = convert_(y);
+            return point_type(xc, yc);
+        }
+
+        template <typename P>
+        point_type prepare_point(const P& p) {
+            return prepare_point(p.x(), p.y());
+        }
+
         // Remove degenerate edge.
         void remove_edge(voronoi_edge_type *edge) {
             voronoi_vertex_type *vertex1 = edge->vertex0();
@@ -1156,8 +779,9 @@
         int num_edge_records_;
         int num_vertex_records_;
 
-        BRect<coordinate_type> voronoi_rect_;
-        detail::ulp_comparison<T> ulp_cmp;
+        bounding_rectangle_type voronoi_rect_;
+        ctype_converter_type convert_;
+        detail::ulp_comparison<T> ulp_cmp_;
 
         // Disallow copy constructor and operator=
         voronoi_diagram(const voronoi_diagram&);

Added: sandbox/gtl/boost/polygon/voronoi_utils.hpp
==============================================================================
--- (empty file)
+++ sandbox/gtl/boost/polygon/voronoi_utils.hpp	2012-02-02 18:09:48 EST (Thu, 02 Feb 2012)
@@ -0,0 +1,450 @@
+// Boost.Polygon library voronoi_utils.hpp header file
+
+//          Copyright Andrii Sydorchuk 2010-2012.
+// Distributed under the Boost Software License, Version 1.0.
+//    (See accompanying file LICENSE_1_0.txt or copy at
+//          http://www.boost.org/LICENSE_1_0.txt)
+
+// See http://www.boost.org for updates, documentation, and revision history.
+
+#ifndef BOOST_POLYGON_VORONOI_UTILS
+#define BOOST_POLYGON_VORONOI_UTILS
+
+#include <cmath>
+#include <stack>
+#include <vector>
+
+#include "polygon.hpp"
+#include "voronoi_diagram.hpp"
+
+namespace boost {
+namespace polygon {
+
+template <typename fpt_>
+struct voronoi_utils_traits;
+
+template<>
+struct voronoi_utils_traits<double> {
+    typedef double fpt_type;
+    typedef point_data<fpt_type> point_type;
+    typedef std::vector<point_type> point_set_type;
+    typedef directed_line_segment_data<fpt_type> segment_type;
+    typedef directed_line_segment_set_data<fpt_type> segment_set_type;
+    typedef bounding_rectangle<fpt_type> brect_type;
+    typedef struct {
+        template <typename CT>
+        fpt_type operator()(const CT& value) {
+            return static_cast<fpt_type>(value);
+        }
+    } ctype_converter_type;
+};
+
+// Voronoi output postprocessing tools.
+template <typename T, typename TRAITS = voronoi_utils_traits<T> >
+class voronoi_utils {
+public:
+    typedef typename TRAITS::fpt_type fpt_type;
+    typedef typename TRAITS::point_type point_type;
+    typedef typename TRAITS::point_set_type point_set_type;
+    typedef typename TRAITS::segment_type segment_type;
+    typedef typename TRAITS::segment_set_type segment_set_type;
+    typedef typename TRAITS::brect_type brect_type;
+    typedef typename TRAITS::ctype_converter_type ctype_converter_type;
+
+    // There are three different types of edges:
+    //   1) Segment edge - has both endpoints;
+    //   2) Ray edge - has one endpoint, infinite
+    //                 in the positive direction;
+    //   3) Line edge - is infinite in both directions.
+    enum kEdgeType {
+        SEGMENT = 0,
+        RAY = 1,
+        LINE = 2
+    };
+
+    // Get a view rectangle based on the voronoi bounding rectangle.
+    template <typename T>
+    static brect_type get_view_rectangle(const bounding_rectangle<T> &brect,
+                                         fpt_type scale = 1.0) {
+        fpt_type x_len = to_fpt(brect.x_max()) - to_fpt(brect.x_min());
+        fpt_type y_len = to_fpt(brect.y_max()) - to_fpt(brect.y_min());
+        fpt_type x_mid = to_fpt(brect.x_max()) + to_fpt(brect.x_min());
+        fpt_type y_mid = to_fpt(brect.y_max()) + to_fpt(brect.y_min());
+        x_mid *= to_fpt(0.5);
+        y_mid *= to_fpt(0.5);
+        fpt_type offset = (std::max)(x_len, y_len) * scale * to_fpt(0.5);
+        brect_type new_brect(x_mid - offset, y_mid - offset,
+                             x_mid + offset, y_mid + offset);
+        return new_brect;
+    }
+
+    // Compute intermediate points for the voronoi edge within the given
+    // bounding rectangle. The assumption is made that voronoi rectangle
+    // contains all the finite endpoints of the edge.
+    // Max_error is the maximum distance (relative to the minor of two
+    // rectangle sides) allowed between the parabola and line segments
+    // that interpolate it.
+    template <typename T>
+    static point_set_type get_point_interpolation(
+        const voronoi_edge<T> *edge,
+        const bounding_rectangle<T> &brect,
+        fpt_type max_error) {
+        // Retrieve the cell to the left of the current edge.
+        const voronoi_cell<T> *cell1 = edge->cell();
+
+        // Retrieve the cell to the right of the current edge.
+        const voronoi_cell<T> *cell2 = edge->twin()->cell();
+
+        // edge_points - contains intermediate points.
+        point_set_type edge_points;
+        if (!(cell1->contains_segment() ^ cell2->contains_segment())) {
+            // Edge is a segment, ray, or line.
+            if (edge->vertex0() != NULL && edge->vertex1() != NULL) {
+                // Edge is a segment. No further processing is required.
+                edge_points.push_back(
+                    get_point(edge->vertex0()->vertex()));
+                edge_points.push_back(
+                    get_point(edge->vertex1()->vertex()));
+            } else {
+                // Edge is a ray or line.
+                // Clip it with the bounding rectangle.
+                const point_type &site1 = get_point(cell1->point0());
+                const point_type &site2 = get_point(cell2->point0());
+                point_type origin((site1.x() + site2.x()) * to_fpt(0.5),
+                                  (site1.y() + site2.y()) * to_fpt(0.5));
+                point_type direction(site1.y() - site2.y(),
+                                     site2.x() - site1.x());
+
+                // Find intersection points.
+                find_intersections(origin, direction, LINE,
+                                   brect, edge_points);
+
+                // Update endpoints in case edge is a ray.
+                if (edge->vertex1() != NULL)
+                    edge_points[1] = get_point(edge->vertex1()->vertex());
+                if (edge->vertex0() != NULL)
+                    edge_points[0] = get_point(edge->vertex0()->vertex());
+            }
+        } else {
+            // point1 - site point;
+            const point_type &point1 = (cell1->contains_segment()) ?
+                get_point(cell2->point0()) : get_point(cell1->point0());
+
+            // point2 - startpoint of the segment site;
+            const point_type &point2 = (cell1->contains_segment()) ?
+                get_point(cell1->point0()) : get_point(cell2->point0());
+
+            // point3 - endpoint of the segment site;
+            const point_type &point3 = (cell1->contains_segment()) ?
+                get_point(cell1->point1()) : get_point(cell2->point1());
+
+            if (edge->vertex0() != NULL && edge->vertex1() != NULL) {
+                // Edge is a segment or parabolic arc.
+                edge_points.push_back(
+                    get_point(edge->vertex0()->vertex()));
+                edge_points.push_back(
+                    get_point(edge->vertex1()->vertex()));
+                fpt_type max_dist = max_error * brect.min_len();
+                fill_intermediate_points(point1, point2, point3,
+                                         edge_points, max_dist);
+            } else {
+                // Edge is a ray or line.
+                // Clip it with the bounding rectangle.
+                fpt_type dir_x =
+                    (cell1->contains_segment() ^ (point1 == point3)) ?
+                    point3.y() - point2.y() : point2.y() - point3.y();
+                fpt_type dir_y =
+                    (cell1->contains_segment() ^ (point1 == point3)) ?
+                    point2.x() - point3.x() : point3.x() - point2.x();
+                point_type direction(dir_x, dir_y);
+
+                // Find intersection points.
+                find_intersections(point1, direction, LINE,
+                                   brect, edge_points);
+
+                // Update endpoints in case edge is a ray.
+                if (edge->vertex1() != NULL)
+                    edge_points[1] = get_point(edge->vertex1()->vertex());
+                if (edge->vertex0() != NULL)
+                    edge_points[0] = get_point(edge->vertex0()->vertex());
+            }
+        }
+        return edge_points;
+    }
+
+    // Interpolate voronoi edge with a set of segments to satisfy maximal
+    // error requirement.
+    template <typename T>
+    static segment_set_type get_segment_interpolation(
+        const voronoi_edge<T> *edge,
+        const bounding_rectangle<T> &brect,
+        fpt_type max_error) {
+        point_set_type point_interpolation =
+            get_point_interpolcation(edge, brect, max_error);
+        segment_set_type ret_val;
+        for (size_t i = 1; i < point_interpolation.size(); ++i)
+            ret_val.insert(segment_type(point_interpolation[i-1],
+                                        point_interpolation[i]));
+        return ret_val;
+    }
+
+    // Find edge-rectangle intersection points.
+    // Edge is represented by its origin, direction and type.
+    static void find_intersections(
+            const point_type &origin, const point_type &direction,
+            kEdgeType edge_type, const brect_type &brect,
+            point_set_type &intersections) {
+        // Find the center of the rectangle.
+        fpt_type center_x = (brect.x_min() + brect.x_max()) * to_fpt(0.5);
+        fpt_type center_y = (brect.y_min() + brect.y_max()) * to_fpt(0.5);
+
+        // Find the half-diagonal vector of the rectangle.
+        fpt_type len_x = brect.x_max() - center_x;
+        fpt_type len_y = brect.y_max() - center_y;
+
+        // Find the vector between the origin and the center of the
+        // rectangle.
+        fpt_type diff_x = origin.x() - center_x;
+        fpt_type diff_y = origin.y() - center_y;
+
+        // Find the vector perpendicular to the direction vector.
+        fpt_type perp_x = direction.y();
+        fpt_type perp_y = -direction.x();
+
+        // Projection of the vector between the origin and the center of
+        // the rectangle on the line perpendicular to the direction vector.
+        fpt_type lexpr = magnitude(perp_x * diff_x + perp_y * diff_y);
+
+        // Maximum projection among any of the half-diagonals of the
+        // rectangle on the line perpendicular to the direction vector.
+        fpt_type rexpr = magnitude(perp_x * len_x) + magnitude(perp_y * len_y);
+
+        // Intersection check. Compare projections.
+        if (lexpr > rexpr)
+            return;
+
+        // Intersection parameters. If fT[i]_used is true than:
+        // origin + fT[i] * direction = intersection point, i = 0 .. 1.
+        // For different edge types next fT values are acceptable:
+        // Segment: [0; 1];
+        // Ray: [0; infinity];
+        // Line: [-infinity; infinity].
+        bool fT0_used = false;
+        bool fT1_used = false;
+        fpt_type fT0 = 0;
+        fpt_type fT1 = 0;
+
+        // Check for the intersections with the lines
+        // going through the sides of the bounding rectangle.
+        clip_line(+direction.x(), -diff_x - len_x,
+                  fT0_used, fT1_used, fT0, fT1);
+        clip_line(-direction.x(), +diff_x - len_x,
+                  fT0_used, fT1_used, fT0, fT1);
+        clip_line(+direction.y(), -diff_y - len_y,
+                  fT0_used, fT1_used, fT0, fT1);
+        clip_line(-direction.y(), +diff_y - len_y,
+                  fT0_used, fT1_used, fT0, fT1);
+
+        // Update intersections vector.
+        if (fT0_used && check_extent(fT0, edge_type))
+            intersections.push_back(point_type(
+                origin.x() + fT0 * direction.x(),
+                origin.y() + fT0 * direction.y()));
+        if (fT1_used && fT0 != fT1 && check_extent(fT1, edge_type))
+            intersections.push_back(point_type(
+                origin.x() + fT1 * direction.x(),
+                origin.y() + fT1 * direction.y()));
+    }
+
+private:
+    voronoi_utils();
+
+    template <typename P>
+    static point_type get_point(const P &point) {
+        fpt_type x = to_fpt(point.x());
+        fpt_type y = to_fpt(point.y());
+        return point_type(x, y);
+    }
+
+    // Find intermediate points of the parabola. Number of points
+    // is defined by the value of max_dist parameter - maximum distance
+    // between parabola and line segments that interpolate it.
+    // Parabola is a locus of points equidistant from the point and segment
+    // sites. intermediate_points should contain two initial endpoints
+    // of the edge (voronoi vertices). Intermediate points are inserted
+    // between the given two endpoints.
+    // Max_dist is the maximum distance allowed between parabola and line
+    // segments that interpolate it.
+    static void fill_intermediate_points(
+            point_type point_site,
+            point_type segment_site_start,
+            point_type segment_site_end,
+            point_set_type &intermediate_points,
+            fpt_type max_dist) {
+        // Check if bisector is a segment.
+        if (point_site == segment_site_start ||
+            point_site == segment_site_end)
+            return;
+
+        // Apply the linear transformation to move start point of the
+        // segment to the point with coordinates (0, 0) and the direction
+        // of the segment to coincide the positive direction of the x-axis.
+        fpt_type segm_vec_x = segment_site_end.x() -
+                              segment_site_start.x();
+        fpt_type segm_vec_y = segment_site_end.y() -
+                              segment_site_start.y();
+        fpt_type sqr_segment_length = segm_vec_x * segm_vec_x +
+                                      segm_vec_y * segm_vec_y;
+
+        // Compute x-coordinates of the endpoints of the edge
+        // in the transformed space.
+        fpt_type projection_start = sqr_segment_length *
+            get_point_projection(intermediate_points[0],
+                                 segment_site_start,
+                                 segment_site_end);
+        fpt_type projection_end = sqr_segment_length *
+            get_point_projection(intermediate_points[1],
+                                 segment_site_start,
+                                 segment_site_end);
+
+        // Compute parabola parameterers in the transformed space.
+        // Parabola has next representation:
+        // f(x) = ((x-rot_x)^2 + rot_y^2) / (2.0*rot_y).
+        fpt_type point_vec_x = point_site.x() - segment_site_start.x();
+        fpt_type point_vec_y = point_site.y() - segment_site_start.y();
+        fpt_type rot_x = segm_vec_x * point_vec_x +
+                         segm_vec_y * point_vec_y;
+        fpt_type rot_y = segm_vec_x * point_vec_y -
+                         segm_vec_y * point_vec_x;
+
+        // Save the last point.
+        point_type last_point = intermediate_points[1];
+        intermediate_points.pop_back();
+
+        // Use stack to avoid recursion.
+        std::stack<fpt_type> point_stack;
+        point_stack.push(projection_end);
+        fpt_type cur_x = projection_start;
+        fpt_type cur_y = parabola_y(cur_x, rot_x, rot_y);
+
+        // Adjust max_dist parameter in the transformed space.
+        max_dist *= max_dist * sqr_segment_length;
+
+        while (!point_stack.empty()) {
+            fpt_type new_x = point_stack.top();
+            fpt_type new_y = parabola_y(new_x, rot_x, rot_y);
+
+            // Compute coordinates of the point of the parabola that is
+            // furthest from the current line segment.
+            fpt_type mid_x = (new_y - cur_y) / (new_x - cur_x) * rot_y +
+                             rot_x;
+            fpt_type mid_y = parabola_y(mid_x, rot_x, rot_y);
+
+            // Compute maximum distance between the given parabolic arc
+            // and line segment that interpolates it.
+            fpt_type dist = (new_y - cur_y) * (mid_x - cur_x) -
+                            (new_x - cur_x) * (mid_y - cur_y);
+            dist = dist * dist / ((new_y - cur_y) * (new_y - cur_y) +
+                                  (new_x - cur_x) * (new_x - cur_x));
+            if (dist <= max_dist) {
+                // Distance between parabola and line segment is
+                // not greater than max_dist.
+                point_stack.pop();
+                fpt_type inter_x =
+                    (segm_vec_x * new_x - segm_vec_y * new_y) /
+                    sqr_segment_length + segment_site_start.x();
+                fpt_type inter_y =
+                    (segm_vec_x * new_y + segm_vec_y * new_x) /
+                    sqr_segment_length + segment_site_start.y();
+                intermediate_points.push_back(
+                    point_type(inter_x, inter_y));
+                cur_x = new_x;
+                cur_y = new_y;
+            } else {
+                point_stack.push(mid_x);
+            }
+        }
+
+        // Update last point.
+        intermediate_points.back() = last_point;
+    }
+
+    // Compute y(x) = ((x - a) * (x - a) + b * b) / (2 * b).
+    static fpt_type parabola_y(fpt_type x, fpt_type a, fpt_type b) {
+        return ((x - a) * (x - a) + b * b) / (to_fpt(2.0) * b);
+    }
+
+    // Check whether extent is compatible with the edge type.
+    static bool check_extent(fpt_type extent, kEdgeType etype) {
+        switch (etype) {
+            case SEGMENT:
+                return extent >= to_fpt(0.0) &&
+                       extent <= to_fpt(1.0);
+            case RAY: return extent >= to_fpt(0.0);
+            case LINE: return true;
+        }
+        return true;
+    }
+
+    // Compute the absolute value.
+    static inline fpt_type magnitude(fpt_type value) {
+        if (value >= to_fpt(0.0))
+            return value;
+        return -value;
+    }
+
+    // Find fT min and max values: fT = numer / denom.
+    static void clip_line(fpt_type denom, fpt_type numer,
+                          bool &fT0_used, bool &fT1_used,
+                          fpt_type &fT0, fpt_type &fT1) {
+        if (denom > to_fpt(0.0)) {
+            if (fT1_used && numer > denom * fT1)
+                return;
+            if (!fT0_used || numer > denom * fT0) {
+                fT0_used = true;
+                fT0 = numer / denom;
+            }
+        } else if (denom < to_fpt(0.0)) {
+            if (fT0_used && numer > denom * fT0)
+                return;
+            if (!fT1_used || numer > denom * fT1) {
+                fT1_used = true;
+                fT1 = numer / denom;
+            }
+        }
+    }
+
+    // Get normalized length of the distance between:
+    //     1) point projection onto the segment;
+    //     2) start point of the segment.
+    // Return this length divided by the segment length.
+    // This is made to avoid sqrt computation during transformation from
+    // the initial space to the transformed one and vice versa.
+    // Assumption is made that projection of the point lies
+    // between the startpoint and endpoint of the segment.
+    static fpt_type get_point_projection(
+            const point_type &point,
+            const point_type &segment_start,
+            const point_type &segment_end) {
+        fpt_type segment_vec_x = segment_end.x() - segment_start.x();
+        fpt_type segment_vec_y = segment_end.y() - segment_start.y();
+        fpt_type point_vec_x = point.x() - segment_start.x();
+        fpt_type point_vec_y = point.y() - segment_start.y();
+        fpt_type sqr_segment_length = segment_vec_x * segment_vec_x +
+                                      segment_vec_y * segment_vec_y;
+        fpt_type vec_dot = segment_vec_x * point_vec_x +
+                           segment_vec_y * point_vec_y;
+        return vec_dot / sqr_segment_length;
+    }
+
+    template <typename CT>
+    static fpt_type to_fpt(const CT& value) {
+        static ctype_converter_type converter;
+        return converter(value);
+    }
+
+};
+}  // polygon
+}  // boost
+
+#endif  // BOOST_POLYGON_VORONOI_UTILS

Modified: sandbox/gtl/libs/polygon/test/voronoi_benchmark_test.cpp
==============================================================================
--- sandbox/gtl/libs/polygon/test/voronoi_benchmark_test.cpp	(original)
+++ sandbox/gtl/libs/polygon/test/voronoi_benchmark_test.cpp	2012-02-02 18:09:48 EST (Thu, 02 Feb 2012)
@@ -19,6 +19,7 @@
 #include <boost/test/test_case_template.hpp>
 #include <boost/timer.hpp>
 
+#include "boost/polygon/polygon.hpp"
 #include "boost/polygon/voronoi.hpp"
 using namespace boost::polygon;
 

Modified: sandbox/gtl/libs/polygon/test/voronoi_builder_test.cpp
==============================================================================
--- sandbox/gtl/libs/polygon/test/voronoi_builder_test.cpp	(original)
+++ sandbox/gtl/libs/polygon/test/voronoi_builder_test.cpp	2012-02-02 18:09:48 EST (Thu, 02 Feb 2012)
@@ -14,6 +14,7 @@
 #include <boost/random/mersenne_twister.hpp>
 #include <boost/test/test_case_template.hpp>
 
+#include "boost/polygon/polygon.hpp"
 #include "boost/polygon/voronoi.hpp"
 using namespace boost::polygon;
 #include "voronoi_test_helper.hpp"

Modified: sandbox/gtl/libs/polygon/test/voronoi_clipping_test.cpp
==============================================================================
--- sandbox/gtl/libs/polygon/test/voronoi_clipping_test.cpp	(original)
+++ sandbox/gtl/libs/polygon/test/voronoi_clipping_test.cpp	2012-02-02 18:09:48 EST (Thu, 02 Feb 2012)
@@ -10,12 +10,12 @@
 #define BOOST_TEST_MODULE voronoi_clipping_test
 #include <boost/test/test_case_template.hpp>
 
-#include "boost/polygon/voronoi_diagram.hpp"
+#include "boost/polygon/voronoi_utils.hpp"
 using namespace boost::polygon;
 
-typedef voronoi_helper<double> VH;
-typedef VH::brect_type rect_type;
-typedef VH::point_type point_type;
+typedef voronoi_utils<double> VU;
+typedef VU::brect_type rect_type;
+typedef VU::point_type point_type;
 
 #define SZ(st) static_cast<int>(st.size())
 
@@ -32,38 +32,38 @@
     point_type direction5(5, -1);
     std::vector<point_type> intersections;
 
-    VH::find_intersections(origin, direction1_1, VH::SEGMENT, rect, intersections);
+    VU::find_intersections(origin, direction1_1, VU::SEGMENT, rect, intersections);
     BOOST_CHECK_EQUAL(SZ(intersections), 2);
     BOOST_CHECK_EQUAL(intersections[0].x() == 0 && intersections[0].y() == 2, true);
     BOOST_CHECK_EQUAL(intersections[1].x() == 3 && intersections[1].y() == -1, true);
     intersections.clear();
     
-    VH::find_intersections(origin, direction1_2, VH::SEGMENT, rect, intersections);
+    VU::find_intersections(origin, direction1_2, VU::SEGMENT, rect, intersections);
     BOOST_CHECK_EQUAL(SZ(intersections), 1);
     BOOST_CHECK_EQUAL(intersections[0].x() == 0 && intersections[0].y() == 2, true);
     intersections.clear();
 
-    VH::find_intersections(origin, direction1_3, VH::SEGMENT, rect, intersections);
+    VU::find_intersections(origin, direction1_3, VU::SEGMENT, rect, intersections);
     BOOST_CHECK_EQUAL(intersections.empty(), true);
     intersections.clear();
 
-    VH::find_intersections(origin, direction2, VH::SEGMENT, rect, intersections);
+    VU::find_intersections(origin, direction2, VU::SEGMENT, rect, intersections);
     BOOST_CHECK_EQUAL(SZ(intersections), 2);
     BOOST_CHECK_EQUAL(intersections[0].x() == 0 && intersections[0].y() == 1, true);
     BOOST_CHECK_EQUAL(intersections[1].x() == 1 && intersections[1].y() == -1, true);
     intersections.clear();
 
-    VH::find_intersections(origin, direction3, VH::SEGMENT, rect, intersections);
+    VU::find_intersections(origin, direction3, VU::SEGMENT, rect, intersections);
     BOOST_CHECK_EQUAL(SZ(intersections), 1);
     BOOST_CHECK_EQUAL(intersections[0].x() == 1 && intersections[0].y() == 2, true);
     intersections.clear();
 
-    VH::find_intersections(origin, direction4, VH::SEGMENT, rect, intersections);
+    VU::find_intersections(origin, direction4, VU::SEGMENT, rect, intersections);
     BOOST_CHECK_EQUAL(SZ(intersections), 1);
     BOOST_CHECK_EQUAL(intersections[0].x() == 0 && intersections[0].y() == -1, true);
     intersections.clear();
 
-    VH::find_intersections(origin, direction5, VH::SEGMENT, rect, intersections);
+    VU::find_intersections(origin, direction5, VU::SEGMENT, rect, intersections);
     BOOST_CHECK_EQUAL(SZ(intersections), 1);
     BOOST_CHECK_EQUAL(intersections[0].x() == 4 && intersections[0].y() == 2, true);
     intersections.clear();
@@ -78,7 +78,7 @@
         for (int j = -50; j <= 50; j++) {
             intersections.clear();
             point_type direction(i, j);
-            VH::find_intersections(origin, direction, VH::SEGMENT, rect, intersections);
+            VU::find_intersections(origin, direction, VU::SEGMENT, rect, intersections);
             if (abs(i) >= 2 || abs(j) >= 2)
                 BOOST_CHECK_EQUAL(SZ(intersections), 1);
             else
@@ -97,7 +97,7 @@
             double x = 1.0 * i / 26;
             double y = 1.0 * j / 26;
             point_type direction(x, y);
-            VH::find_intersections(origin, direction, VH::SEGMENT, rect, intersections);
+            VU::find_intersections(origin, direction, VU::SEGMENT, rect, intersections);
             BOOST_CHECK_EQUAL(SZ(intersections), 0);
         }
 }
@@ -113,30 +113,30 @@
     point_type direction5(5, -1);
     std::vector<point_type> intersections;
 
-    VH::find_intersections(origin, direction1, VH::RAY, rect, intersections);
+    VU::find_intersections(origin, direction1, VU::RAY, rect, intersections);
     BOOST_CHECK_EQUAL(SZ(intersections), 2);
     BOOST_CHECK_EQUAL(intersections[0].x() == 0 && intersections[0].y() == 2, true);
     BOOST_CHECK_EQUAL(intersections[1].x() == 3 && intersections[1].y() == -1, true);
     intersections.clear();
 
-    VH::find_intersections(origin, direction2, VH::RAY, rect, intersections);
+    VU::find_intersections(origin, direction2, VU::RAY, rect, intersections);
     BOOST_CHECK_EQUAL(SZ(intersections), 2);
     BOOST_CHECK_EQUAL(intersections[0].x() == 0 && intersections[0].y() == 1, true);
     BOOST_CHECK_EQUAL(intersections[1].x() == 1 && intersections[1].y() == -1, true);
     intersections.clear();
 
-    VH::find_intersections(origin, direction3, VH::RAY, rect, intersections);
+    VU::find_intersections(origin, direction3, VU::RAY, rect, intersections);
     BOOST_CHECK_EQUAL(SZ(intersections), 2);
     BOOST_CHECK_EQUAL(intersections[0].x() == 1 && intersections[0].y() == 2, true);
     BOOST_CHECK_EQUAL(intersections[1].x() == 4 && intersections[1].y() == 0.5, true);
     intersections.clear();
 
-    VH::find_intersections(origin, direction4, VH::RAY, rect, intersections);
+    VU::find_intersections(origin, direction4, VU::RAY, rect, intersections);
     BOOST_CHECK_EQUAL(SZ(intersections), 1);
     BOOST_CHECK_EQUAL(intersections[0].x() == 0 && intersections[0].y() == -1, true);
     intersections.clear();
 
-    VH::find_intersections(origin, direction5, VH::RAY, rect, intersections);
+    VU::find_intersections(origin, direction5, VU::RAY, rect, intersections);
     BOOST_CHECK_EQUAL(SZ(intersections), 1);
     BOOST_CHECK_EQUAL(intersections[0].x() == 4 && intersections[0].y() == 2, true);
     intersections.clear();
@@ -154,7 +154,7 @@
             double x = 1.0 * i / 26;
             double y = 1.0 * j / 26;
             point_type direction(x, y);
-            VH::find_intersections(origin, direction, VH::RAY, rect, intersections);
+            VU::find_intersections(origin, direction, VU::RAY, rect, intersections);
             BOOST_CHECK_EQUAL(SZ(intersections), 1);
         }
 }
@@ -170,30 +170,30 @@
     point_type direction5(-5, 1);
     std::vector<point_type> intersections;
     
-    VH::find_intersections(origin, direction1, VH::LINE, rect, intersections);
+    VU::find_intersections(origin, direction1, VU::LINE, rect, intersections);
     BOOST_CHECK_EQUAL(SZ(intersections), 2);
     BOOST_CHECK_EQUAL(intersections[0].x() == 3 && intersections[0].y() == -1, true);
     BOOST_CHECK_EQUAL(intersections[1].x() == 0 && intersections[1].y() == 2, true);
     intersections.clear();
 
-    VH::find_intersections(origin, direction2, VH::LINE, rect, intersections);
+    VU::find_intersections(origin, direction2, VU::LINE, rect, intersections);
     BOOST_CHECK_EQUAL(SZ(intersections), 2);
     BOOST_CHECK_EQUAL(intersections[0].x() == 1 && intersections[0].y() == -1, true);
     BOOST_CHECK_EQUAL(intersections[1].x() == 0 && intersections[1].y() == 1, true);
     intersections.clear();
 
-    VH::find_intersections(origin, direction3, VH::LINE, rect, intersections);
+    VU::find_intersections(origin, direction3, VU::LINE, rect, intersections);
     BOOST_CHECK_EQUAL(SZ(intersections), 2);
     BOOST_CHECK_EQUAL(intersections[0].x() == 4 && intersections[0].y() == 0.5, true);
     BOOST_CHECK_EQUAL(intersections[1].x() == 1 && intersections[1].y() == 2, true);
     intersections.clear();
 
-    VH::find_intersections(origin, direction4, VH::LINE, rect, intersections);
+    VU::find_intersections(origin, direction4, VU::LINE, rect, intersections);
     BOOST_CHECK_EQUAL(SZ(intersections), 1);
     BOOST_CHECK_EQUAL(intersections[0].x() == 0 && intersections[0].y() == -1, true);
     intersections.clear();
 
-    VH::find_intersections(origin, direction5, VH::LINE, rect, intersections);
+    VU::find_intersections(origin, direction5, VU::LINE, rect, intersections);
     BOOST_CHECK_EQUAL(SZ(intersections), 1);
     BOOST_CHECK_EQUAL(intersections[0].x() == 4 && intersections[0].y() == 2, true);
     intersections.clear();
@@ -211,7 +211,7 @@
             double x = 1.0 * i / 26;
             double y = 1.0 * j / 26;
             point_type direction(x, y);
-            VH::find_intersections(origin, direction, VH::LINE, rect, intersections);
+            VU::find_intersections(origin, direction, VU::LINE, rect, intersections);
             BOOST_CHECK_EQUAL(SZ(intersections), 2);
         }
 }

Modified: sandbox/gtl/libs/polygon/test/voronoi_test_helper.hpp
==============================================================================
--- sandbox/gtl/libs/polygon/test/voronoi_test_helper.hpp	(original)
+++ sandbox/gtl/libs/polygon/test/voronoi_test_helper.hpp	2012-02-02 18:09:48 EST (Thu, 02 Feb 2012)
@@ -106,20 +106,29 @@
     return true;
 }
 
+template <typename P>
+struct cmp {
+    bool operator()(const P& p1, const P& p2) const {
+        if (p1.x() != p2.x()) return p1.x() < p2.x();
+        return p1.y() < p2.y();
+    }
+};
+
 template <typename Output>
 bool verfiy_no_line_edge_intersections(const Output &output) {
     // Create map from edges with first point less than the second one.
     // Key is the first point of the edge, value is a vector of second points
     // with the same first point.
     typedef typename Output::point_type point_type;
-    std::map< point_type, std::vector<point_type> > edge_map;
+    cmp<point_type> comparator;
+    std::map< point_type, std::vector<point_type>, cmp<point_type> > edge_map;
     typename Output::voronoi_edge_const_iterator_type edge_it;
     for (edge_it = output.edge_records().begin();
          edge_it != output.edge_records().end(); edge_it++) {
         if (edge_it->is_bounded()) {
             const point_type &vertex0 = edge_it->vertex0()->vertex();
             const point_type &vertex1 = edge_it->vertex1()->vertex();
-            if (vertex0 < vertex1)
+            if (comparator(vertex0, vertex1))
                 edge_map[vertex0].push_back(vertex1);
         }
     }
@@ -127,7 +136,7 @@
 }
 
 template <typename Point2D>
-bool intersection_check(const std::map< Point2D, std::vector<Point2D> > &edge_map) {
+bool intersection_check(const std::map< Point2D, std::vector<Point2D>, cmp<Point2D> > &edge_map) {
     // Iterate over map of edges and check if there are any intersections.
     // All the edges are stored by the low x value. That's why we iterate
     // left to right checking for intersections between all pairs of edges
@@ -135,7 +144,7 @@
     // Complexity. Approximately N*sqrt(N). Worst case N^2.
     typedef Point2D point_type;
     typedef typename point_type::coordinate_type coordinate_type;
-    typedef typename std::map< point_type, std::vector<point_type> >::const_iterator
+    typedef typename std::map< point_type, std::vector<point_type>, cmp<Point2D> >::const_iterator
         edge_map_iterator;
     typedef typename std::vector<point_type>::size_type size_type;
     edge_map_iterator edge_map_it1, edge_map_it2, edge_map_it_bound;
@@ -200,23 +209,21 @@
     return result;
 }
 
-template <typename T>
-void save_voronoi_input(const std::vector< point_data<T> > &points, const char *file_name) {
+template <typename PointIterator>
+void save_points(PointIterator first, PointIterator last, const char *file_name) {
     std::ofstream ofs(file_name);
     ofs << points.size() << std::endl;
-    for (typename std::vector< point_data<T> >::iterator it = points.begin();
-         it != points.end(); ++it) {
+    for (PointIterator it = first; it != last; ++it) {
         ofs << it->x() << " " << it->y() << std::endl;
     }
     ofs.close();
 }
 
-template <typename T>
-void save_voronoi_input(const directed_line_segment_set_data<T> &segments, const char *file_name) {
+template <typename SegmentIterator>
+void save_segments(SegmentIterator first, SegmentIterator last, const char *file_name) {
     std::ofstream ofs(file_name);
     ofs << segments.size() << std::endl;
-    for (typename directed_line_segment_set_data<T>::iterator_type it = segments.begin();
-         it != segments.end(); ++it) {
+    for (SegmentIterator it = first; it != last; ++it) {
         ofs << it->low().x() << " " << it->low().y() << " ";
         ofs << it->high().x() << " " << it->high().y() << std::endl;
     }

Modified: sandbox/gtl/libs/polygon/voronoi_example/voronoi_visualizer.cpp
==============================================================================
--- sandbox/gtl/libs/polygon/voronoi_example/voronoi_visualizer.cpp	(original)
+++ sandbox/gtl/libs/polygon/voronoi_example/voronoi_visualizer.cpp	2012-02-02 18:09:48 EST (Thu, 02 Feb 2012)
@@ -13,6 +13,7 @@
 #include <QtGui/QtGui>
 
 #include "boost/polygon/voronoi.hpp"
+#include "boost/polygon/voronoi_utils.hpp"
 using namespace boost::polygon;
 
 class GLWidget : public QGLWidget {
@@ -59,8 +60,8 @@
         // Build voronoi diagram.
         vd_.clear();
         construct_voronoi(point_sites, segment_sites, &vd_);
-        brect_ = voronoi_helper<coordinate_type>::get_view_rectangle(
-            vd_.bounding_rectangle());
+        brect_ = voronoi_utils<coordinate_type>::get_view_rectangle(
+            vd_.bounding_rectangle(), 2.0);
 
         // Update view.
         update_view_port();
@@ -128,7 +129,7 @@
                 if (!it->is_primary() && primary_edges_only_)
                     continue;
                 std::vector<opoint_type> temp_v =
-                    voronoi_helper<coordinate_type>::get_point_interpolation(
+                    voronoi_utils<coordinate_type>::get_point_interpolation(
                         &(*it), brect_, 1E-3);
                 for (int i = 0; i < static_cast<int>(temp_v.size()) - 1; i++) {
                     glVertex2f(temp_v[i].x(), temp_v[i].y());
@@ -176,7 +177,7 @@
         voronoi_vertex_const_iterator_type;
     typedef voronoi_diagram<coordinate_type>::voronoi_edge_const_iterator_type
         voronoi_edge_const_iterator_type;
-    BRect<coordinate_type> brect_;
+    bounding_rectangle<coordinate_type> brect_;
     voronoi_diagram<coordinate_type> vd_;
     bool primary_edges_only_;
 };