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_;
};