|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r67360 - in sandbox/SOC/2010/sweepline: boost/sweepline boost/sweepline/detail libs/sweepline/example libs/sweepline/test
From: sydorchuk.andriy_at_[hidden]
Date: 2010-12-20 09:30:59
Author: asydorchuk
Date: 2010-12-20 09:30:57 EST (Mon, 20 Dec 2010)
New Revision: 67360
URL: http://svn.boost.org/trac/boost/changeset/67360
Log:
Code refactoring.
Added comments.
Fixed (segment, segment) beach line comparison predicate.
Text files modified:
sandbox/SOC/2010/sweepline/boost/sweepline/detail/voronoi_formation.hpp | 1721 ++++++++++++++++++++++++---------------
sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_output.hpp | 927 +++++++++++++--------
sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_sweepline.hpp | 33
sandbox/SOC/2010/sweepline/libs/sweepline/example/voronoi_visualizer.cpp | 28
sandbox/SOC/2010/sweepline/libs/sweepline/test/benchmark_test.cpp | 2
sandbox/SOC/2010/sweepline/libs/sweepline/test/circle_event_creation_test.cpp | 26
sandbox/SOC/2010/sweepline/libs/sweepline/test/clipping_test.cpp | 32
sandbox/SOC/2010/sweepline/libs/sweepline/test/event_queue_test.cpp | 4
sandbox/SOC/2010/sweepline/libs/sweepline/test/event_types_test.cpp | 12
sandbox/SOC/2010/sweepline/libs/sweepline/test/node_comparer_test.cpp | 10
sandbox/SOC/2010/sweepline/libs/sweepline/test/output_verification.hpp | 128 +-
sandbox/SOC/2010/sweepline/libs/sweepline/test/predicates_test.cpp | 54
sandbox/SOC/2010/sweepline/libs/sweepline/test/sweepline_test.cpp | 97 +
sandbox/SOC/2010/sweepline/libs/sweepline/test/test_type_list.hpp | 2
14 files changed, 1878 insertions(+), 1198 deletions(-)
Modified: sandbox/SOC/2010/sweepline/boost/sweepline/detail/voronoi_formation.hpp
==============================================================================
--- sandbox/SOC/2010/sweepline/boost/sweepline/detail/voronoi_formation.hpp (original)
+++ sandbox/SOC/2010/sweepline/boost/sweepline/detail/voronoi_formation.hpp 2010-12-20 09:30:57 EST (Mon, 20 Dec 2010)
@@ -1,23 +1,15 @@
-// Boost sweepline/voronoi_formation.hpp header file
+// Boost sweepline/voronoi_formation.hpp header file
// Copyright Andrii Sydorchuk 2010.
// 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.
+// See http://www.boost.org for updates, documentation, and revision history.
#ifndef BOOST_SWEEPLINE_VORONOI_FORMATION
#define BOOST_SWEEPLINE_VORONOI_FORMATION
-#define INT_PREDICATE_COMPUTE_DIFFERENCE(a, b, res, sign) \
- if (a >= b) { res = static_cast<unsigned long long>(a - b); sign = true; } \
- else { res = static_cast<unsigned long long>(b - a); sign = false; }
-
-#define INT_PREDICATE_CONVERT_65_BIT(a, res, sign) \
- if (a >= 0) { res = static_cast<unsigned long long>(a); sign = true; } \
- else { res = static_cast<unsigned long long>(-a); sign = false; }
-
namespace boost {
namespace sweepline {
namespace detail {
@@ -33,52 +25,76 @@
template <typename T>
class beach_line_node_data;
- template <typename BeachLineNode>
+ template <typename T>
struct node_comparer;
template <typename T>
class epsilon_robust_comparator;
+ // Represents the result of the orientation test.
enum kOrientation {
RIGHT_ORIENTATION = -1,
- COLINEAR = 0,
+ COLLINEAR = 0,
LEFT_ORIENTATION = 1,
};
-
+
+ // Represents the result of the epsilon robust predicate.
+ // If the result is undefined some further processing is usually required.
enum kPredicateResult {
LESS = -1,
UNDEFINED = 0,
MORE = 1,
};
- // Site event type.
- // Occurs when sweepline sweeps over one of the initial sites.
- // Contains index of a site among the other sorted sites.
+ // Site event type.
+ // Occurs when the sweepline sweeps over one of the initial sites:
+ // 1) point site;
+ // 2) startpoint of the segment site;
+ // 3) endpoint of the segment site.
+ // Implicit segment direction is defined: the startpoint of
+ // the segment compares less than its endpoint.
+ // Each input segment is divided onto two site events:
+ // 1) One going from the startpoint to the endpoint
+ // (is_inverse_ = false);
+ // 2) Another going from the endpoint to the startpoint
+ // (is_inverse_ = true).
+ // In beach line data structure segment sites of the first
+ // type preceed sites of the second type for the same segment.
+ // Variables: point0_ - point site or segment's startpoint;
+ // point1_ - segment's endpoint if site is a segment;
+ // index_ - site event index among other sites;
+ // is_segment_ - flag whether site is a segment;
+ // is_vertical_ - flag whether site is vertical;
+ // is_inverse_ - defines type of the segment site.
+ // Note: for all the point sites is_vertical_ flag is true,
+ // is_inverse_ flag is false.
template <typename T>
class site_event {
public:
typedef T coordinate_type;
- typedef point_2d<T> Point2D;
+ typedef point_2d<T> point_2d_type;
+ // Point site contructors.
site_event() :
is_segment_(false),
is_vertical_(true),
is_inverse_(false) {}
-
+
site_event(coordinate_type x, coordinate_type y, int index) :
point0_(x, y),
- site_index_(index),
+ site_index_(index),
is_segment_(false),
is_vertical_(true),
is_inverse_(false) {}
- site_event(const Point2D &point, int index) :
+ site_event(const point_2d_type &point, int index) :
point0_(point),
site_index_(index),
is_segment_(false),
is_vertical_(true),
is_inverse_(false) {}
+ // Segment site constructors.
site_event(coordinate_type x1, coordinate_type y1,
coordinate_type x2, coordinate_type y2, int index):
point0_(x1, y1),
@@ -91,7 +107,8 @@
is_vertical_ = (point0_.x() == point1_.x());
}
- site_event(const Point2D &point1, const Point2D &point2, int index) :
+ site_event(const point_2d_type &point1,
+ const point_2d_type &point2, int index) :
point0_(point1),
point1_(point2),
site_index_(index),
@@ -102,120 +119,152 @@
is_vertical_ = (point0_.x() == point1_.x());
}
- bool operator==(const site_event &s_event) const {
- return (point0_ == s_event.point0_) &&
- ((!is_segment_ && !s_event.is_segment_) ||
- (is_segment_ && s_event.is_segment_ && (point1_ == s_event.point1_)));
- }
-
- bool operator!=(const site_event &s_event) const {
- return !((*this) == s_event);
- }
-
- // All the sites are sorted due to coordinates of the first point.
- // Point sites preceed vertical segment sites with the same first point.
- // Point sites and vertical segments preceed not vertical segments with
- // the same x coordinate of the first point.
- // Non vertical segments with the same first point are sorted counterclockwise.
- bool operator<(const site_event &s_event) const {
- // If first points have different x coordinates, compare them.
- if (this->point0_.x() != s_event.point0_.x())
- return this->point0_.x() < s_event.point0_.x();
+ bool operator==(const site_event &that) const {
+ return (point0_ == that.point0_) &&
+ ((!is_segment_ && !that.is_segment_) ||
+ (is_segment_ && that.is_segment_ &&
+ point1_ == that.point1_));
+ }
+
+ bool operator!=(const site_event &that) const {
+ return !((*this) == that);
+ }
+
+ // All the sites are sorted by the coordinates of the first point.
+ // Point sites preceed vertical segment sites with equal first points.
+ // Point sites and vertical segments preceed non-vertical segments
+ // with the same x-coordinate of the first point (for any y).
+ // Non-vertical segments with equal first points are ordered using
+ // counterclockwise direction starting from the positive direction of
+ // the x-axis. Such ordering simplifies the initialization step of the
+ // algorithm.
+ bool operator<(const site_event &that) const {
+ // Compare x-coordinates of the first points.
+ if (this->point0_.x() != that.point0_.x())
+ return this->point0_.x() < that.point0_.x();
+ // X-coordinates of the first points are the same.
if (!(this->is_segment_)) {
- if (!s_event.is_segment_) {
- return this->point0_.y() < s_event.point0_.y();
+ // The first site is a point.
+ if (!that.is_segment_) {
+ // The second site is a point.
+ return this->point0_.y() < that.point0_.y();
}
- if (s_event.is_vertical_) {
- return this->point0_.y() <= s_event.point0_.y();
+ if (that.is_vertical_) {
+ // The second site is a vertical segment.
+ return this->point0_.y() <= that.point0_.y();
}
+ // The second site is a non-vertical segment.
return true;
} else {
- if (!s_event.is_segment_ || s_event.is_vertical_) {
+ // The first site is a segment.
+ if (that.is_vertical_) {
+ // The second site is a point or a vertical segment.
if (this->is_vertical_) {
- return this->point0_.y() < s_event.point0_.y();
+ // The first site is a vertical segment.
+ return this->point0_.y() < that.point0_.y();
}
+ // The first site is a non-vertical segment.
return false;
}
- if (this->is_vertical_)
+
+ // The second site is a non-vertical segment.
+ if (this->is_vertical_) {
+ // The first site is a vertical segment.
return true;
- if (this->point0_.y() != s_event.point0_.y())
- return this->point0_.y() < s_event.point0_.y();
- // Sort by angle.
- return orientation_test(this->point1_, this->point0_, s_event.point1_) ==
- LEFT_ORIENTATION;
- }
+ }
+
+ // Compare y-coordinates of the first points.
+ if (this->point0_.y() != that.point0_.y())
+ return this->point0_.y() < that.point0_.y();
+
+ // Y-coordinates of the first points are the same.
+ // Both sites are segment. Sort by angle using CCW ordering
+ // starting at the positive direction of the x axis.
+ return orientation_test(this->point1_, this->point0_,
+ that.point1_) == LEFT_ORIENTATION;
+ }
}
- bool operator<=(const site_event &s_event) const {
- return !(s_event < (*this));
+ bool operator<=(const site_event &that) const {
+ return !(that < (*this));
}
- bool operator>(const site_event &s_event) const {
- return s_event < (*this);
+ bool operator>(const site_event &that) const {
+ return that < (*this);
}
- bool operator>=(const site_event &s_event) const {
- return !((*this) < s_event);
+ bool operator>=(const site_event &that) const {
+ return !((*this) < that);
}
coordinate_type x(bool oriented = false) const {
- if (!oriented)
- return point0_.x();
- return is_inverse_ ? point1_.x() : point0_.x();
+ return x0(oriented);
}
coordinate_type y(bool oriented = false) const {
- if (!oriented)
- return point0_.y();
- return is_inverse_ ? point1_.y() : point0_.y();
+ return y0(oriented);
}
+ // Return x-coordinate of the point for the point sites.
+ // Return x-coordinate of the startpoint for the segment sites.
coordinate_type x0(bool oriented = false) const {
if (!oriented)
return point0_.x();
return is_inverse_ ? point1_.x() : point0_.x();
}
+ // Return y-coordinate of the point for the point sites.
+ // Return y-coordinate of the startpoint for the segment sites.
coordinate_type y0(bool oriented = false) const {
if (!oriented)
return point0_.y();
return is_inverse_ ? point1_.y() : point0_.y();
}
+ // Return x-coordinate of the endpoint of the segment sites.
+ // Doesn't make sense for the point sites.
coordinate_type x1(bool oriented = false) const {
if (!oriented)
return point1_.x();
return is_inverse_ ? point0_.x() : point1_.x();
}
+ // Return y-coordinate of the endpoint of the segment sites.
+ // Doesn't make sense for the point sites.
coordinate_type y1(bool oriented = false) const {
if (!oriented)
return point1_.y();
return is_inverse_ ? point0_.y() : point1_.y();
}
- const Point2D &get_point0(bool oriented = false) const {
+ // Return point for the point sites.
+ // Return startpoint for the segment sites.
+ const point_2d_type &point0(bool oriented = false) const {
if (!oriented)
return point0_;
return is_inverse_ ? point1_ : point0_;
}
- const Point2D &get_point1(bool oriented = false) const {
+ // Return endpoint of the segment sites.
+ // Doesn't make sense for the point sites.
+ const point_2d_type &point1(bool oriented = false) const {
if (!oriented)
return point1_;
return is_inverse_ ? point0_ : point1_;
}
- void set_site_index(int index) {
+ // Return index of the site among all the other sites.
+ void index(int index) {
site_index_ = index;
}
- void set_inverse() {
+ // Change segment site orientation to the opposite one.
+ void inverse() {
is_inverse_ ^= true;
}
- int get_site_index() const {
+ int index() const {
return site_index_;
}
@@ -232,8 +281,8 @@
}
private:
- Point2D point0_;
- Point2D point1_;
+ point_2d_type point0_;
+ point_2d_type point1_;
int site_index_;
bool is_segment_;
bool is_vertical_;
@@ -241,193 +290,248 @@
};
template <typename T>
- site_event<T> make_site_event(T x, T y, int index) {
+ static inline site_event<T> make_site_event(T x, T y, int index) {
return site_event<T>(x, y, index);
}
template <typename T>
- site_event<T> make_site_event(const point_2d<T> &point, int index) {
+ static inline site_event<T> make_site_event(const point_2d<T> &point,
+ int index) {
return site_event<T>(point, index);
}
template <typename T>
- site_event<T> make_site_event(T x1, T y1, T x2, T y2, int index) {
+ static inline site_event<T> make_site_event(T x1, T y1, T x2, T y2,
+ int index) {
return site_event<T>(x1, y1, x2, y2, index);
}
template <typename T>
- site_event<T> make_site_event(const point_2d<T> &point1,
- const point_2d<T> &point2, int index) {
+ static inline site_event<T> make_site_event(const point_2d<T> &point1,
+ const point_2d<T> &point2,
+ int index) {
return site_event<T>(point1, point2, index);
}
- // Circle event type. Occurs when sweepline sweeps over the bottom point of
- // the voronoi circle (with center at the bisectors intersection point).
- // Circle event contains circle center, lowest x coordinate, event state and
- // iterator to the corresponding bisectors.
+ // Circle event type.
+ // Occrus when the sweepline sweeps over the rightmost point of the voronoi
+ // circle (with the center at the intersection point of the bisectors).
+ // Circle event is made by the two consequtive nodes in the beach line data
+ // structure. In case another node was inserted during algorithm execution
+ // between the given two nodes circle event becomes inactive.
+ // Circle events order is based on the comparison of the rightmost points
+ // of the circles. The order is the same like for the point_2d class.
+ // However as coordinates of the circle events might be not integers extra
+ // comparison checks are required to make the comparison robust.
+ // Next representation is used to store parameters of the circle:
+ // each of the parameters is represented as fraction
+ // numerator / denominator. Denominator is common for each of the
+ // three parameters. Epsilon robust comparator class is used
+ // to represent parameters of the circle events.
+ // Variables: center_x_ - numerator of the center's x-coordinate.
+ // center_y_ - numerator of the center's y-coordinate.
+ // lower_x_ - numerator of the bottom point's x-coordinate.
+ // denom_ - positive denominator for the previous three values.
+ // bisector_node_ - iterator to the second bisector's node.
+ // is_active_ - flag whether circle event is still active.
+ // lower_y coordinate is always equal to center_y.
template <typename T>
class circle_event {
public:
typedef T coordinate_type;
- typedef point_2d<T> Point2D;
+ typedef point_2d<T> point_2d_type;
+ typedef site_event<T> site_event_type;
+ typedef epsilon_robust_comparator<T> erc_type;
typedef beach_line_node<coordinate_type> Key;
typedef beach_line_node_data<coordinate_type> Value;
typedef node_comparer<Key> NodeComparer;
- typedef typename std::map< Key, Value, NodeComparer >::iterator beach_line_iterator;
+ typedef typename std::map< Key, Value, NodeComparer >::iterator
+ beach_line_iterator;
circle_event() : is_active_(true) {}
- circle_event(coordinate_type c_x, coordinate_type c_y, coordinate_type lower_x) :
+ circle_event(coordinate_type c_x,
+ coordinate_type c_y,
+ coordinate_type lower_x) :
center_x_(c_x),
center_y_(c_y),
lower_x_(lower_x),
denom_(1),
is_active_(true) {}
- circle_event(const epsilon_robust_comparator<T> &c_x,
- const epsilon_robust_comparator<T> &c_y,
- const epsilon_robust_comparator<T> &lower_x) :
+ circle_event(const erc_type &c_x,
+ const erc_type &c_y,
+ const erc_type &lower_x) :
center_x_(c_x),
center_y_(c_y),
lower_x_(lower_x),
denom_(1),
is_active_(true) {}
- circle_event(const epsilon_robust_comparator<T> &c_x,
- const epsilon_robust_comparator<T> &c_y,
- const epsilon_robust_comparator<T> &lower_x,
- const epsilon_robust_comparator<T> &denom) :
- center_x_(c_x),
- center_y_(c_y),
- lower_x_(lower_x),
- denom_(denom),
- is_active_(true) {
- if (denom_.abs()) {
- center_x_.swap();
- center_y_.swap();
- lower_x_.swap();
- }
+ circle_event(const erc_type &c_x,
+ const erc_type &c_y,
+ const erc_type &lower_x,
+ const erc_type &denom) :
+ center_x_(c_x),
+ center_y_(c_y),
+ lower_x_(lower_x),
+ denom_(denom),
+ is_active_(true) {
+ if (denom_.abs()) {
+ center_x_.swap();
+ center_y_.swap();
+ lower_x_.swap();
}
-
- bool operator==(const circle_event &c_event) const {
- const epsilon_robust_comparator<T> &lhs1 = center_y_ * c_event.denom_;
- const epsilon_robust_comparator<T> &rhs1 = denom_ * c_event.center_y_;
- const epsilon_robust_comparator<T> &lhs2 = lower_x_ * c_event.denom_;
- const epsilon_robust_comparator<T> &rhs2 = denom_ * c_event.lower_x_;
- return (lhs1.compare(rhs1) == UNDEFINED && lhs2.compare(rhs2) == UNDEFINED);
}
- bool operator!=(const circle_event &c_event) const {
- return !((*this) == c_event);
- }
+ bool operator==(const circle_event &that) const {
+ erc_type lhs1 = lower_x_ * that.denom_;
+ erc_type rhs1 = denom_ * that.lower_x_;
+ erc_type lhs2 = center_y_ * that.denom_;
+ erc_type rhs2 = denom_ * that.center_y_;
+ return (lhs1.compare(rhs1, 64) == UNDEFINED &&
+ lhs2.compare(rhs2, 64) == UNDEFINED);
+ }
+
+ bool operator!=(const circle_event &that) const {
+ return !((*this) == that);
+ }
+
+ // Compare rightmost points of the circle events.
+ // Each rightmost point has next representation:
+ // (lower_x / denom, center_y / denom), denom is always positive.
+ bool operator<(const circle_event &that) const {
+ // Evaluate comparison expressions for x-coordinates.
+ erc_type lhs1 = lower_x_ * that.denom_;
+ erc_type rhs1 = denom_ * that.lower_x_;
- bool operator<(const circle_event &c_event) const {
- const epsilon_robust_comparator<T> &lhs1 = lower_x_ * c_event.denom_;
- const epsilon_robust_comparator<T> &rhs1 = denom_ * c_event.lower_x_;
+ // Compare x-coordinates of the rightmost points. If the result
+ // is undefined we assume that x-coordinates are equal.
kPredicateResult pres = lhs1.compare(rhs1, 64);
if (pres != UNDEFINED)
return (pres == LESS);
- const epsilon_robust_comparator<T> &lhs2 = center_y_ * c_event.denom_;
- const epsilon_robust_comparator<T> &rhs2 = denom_ * c_event.center_y_;
+
+ // Evaluate comparison expressions for y-coordinates.
+ erc_type lhs2 = center_y_ * that.denom_;
+ erc_type rhs2 = denom_ * that.center_y_;
+
+ // Compare y-coordinates of the rightmost points.
return (lhs2.compare(rhs2, 64) == LESS);
}
- bool operator<=(const circle_event &c_event) const {
- return !(c_event < (*this));
+ bool operator<=(const circle_event &that) const {
+ return !(that < (*this));
}
- bool operator>(const circle_event &c_event) const {
- return c_event < (*this);
+ bool operator>(const circle_event &that) const {
+ return that < (*this);
}
- bool operator>=(const circle_event &c_event) const {
- return !((*this) < c_event);
+ bool operator>=(const circle_event &that) const {
+ return !((*this) < that);
}
- // Compares bottom voronoi circle point with site event point.
+ // Compare the rightmost point of the circle event with
+ // the point that represents the site event.
// If circle point is less than site point return -1.
// If circle point is equal to site point return 0.
// If circle point is greater than site point return 1.
- kPredicateResult compare(const site_event<coordinate_type> &s_event) const {
+ kPredicateResult compare(const site_event_type &s_event) const {
+ // Compare x-coordinates.
kPredicateResult pres = lower_x_.compare(denom_ * s_event.x(), 64);
+ // If the comparison result is undefined
+ // x-coordinates are considered to be equal.
if (pres != UNDEFINED)
return pres;
+ // Compare y-coordinates in case of equal x-coordinates.
return center_y_.compare(denom_ * s_event.y(), 64);
}
+ // Evaluate x-coordinate of the center of the circle.
coordinate_type x() const {
return center_x_.dif() / denom_.dif();
}
+ // Evaluate y-coordinate of the center of the circle.
coordinate_type y() const {
return center_y_.dif() / denom_.dif();
}
+ // Evaluate x-coordinate of the rightmost point of the circle.
coordinate_type lower_x() const {
return lower_x_.dif() / denom_.dif();
}
- Point2D get_center() const {
+ point_2d_type center() const {
return make_point_2d(x(), y());
}
- const epsilon_robust_comparator<T> &get_erc_x() const {
+ const erc_type &erc_x() const {
return center_x_;
}
- const epsilon_robust_comparator<T> &get_erc_y() const {
+ const erc_type &erc_y() const {
return center_y_;
}
- const epsilon_robust_comparator<T> &get_erc_denom() const {
+ const erc_type &erc_denom() const {
return denom_;
}
- void set_bisector(beach_line_iterator iterator) {
- bisector_node_ = iterator;
- }
-
- void deactivate() {
- is_active_ = false;
+ // Return iterator to the second bisector from the beach line
+ // data structure that forms current circle event.
+ const beach_line_iterator &bisector() const {
+ return bisector_node_;
}
- const beach_line_iterator &get_bisector() const {
- return bisector_node_;
+ void bisector(beach_line_iterator iterator) {
+ bisector_node_ = iterator;
}
bool is_active() const {
return is_active_;
}
+ void deactivate() {
+ is_active_ = false;
+ }
+
private:
- epsilon_robust_comparator<T> center_x_;
- epsilon_robust_comparator<T> center_y_;
- epsilon_robust_comparator<T> lower_x_;
- epsilon_robust_comparator<T> denom_;
+ erc_type center_x_;
+ erc_type center_y_;
+ erc_type lower_x_;
+ erc_type denom_;
beach_line_iterator bisector_node_;
bool is_active_;
};
template <typename T>
- circle_event<T> make_circle_event(T c_x, T c_y, T lower_x) {
+ static inline circle_event<T> make_circle_event(T c_x, T c_y, T lower_x) {
return circle_event<T>(c_x, c_y, lower_x);
}
///////////////////////////////////////////////////////////////////////////
// VORONOI CIRCLE EVENTS QUEUE ////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////
-
- // Event queue data structure, processes circle events.
+
+ // Event queue data structure, holds circle events.
+ // During algorithm run, some of the circle events disappear(become
+ // inactive). Priority queue data structure by itself doesn't support
+ // iterators (there is no direct ability to modify its elements).
+ // Instead list is used to store all the circle events and priority queue
+ // of the iterators to the list elements is used to keep the correct circle
+ // events ordering.
template <typename T>
class circle_events_queue {
public:
typedef T coordinate_type;
typedef circle_event<T> circle_event_type;
- typedef typename std::list<circle_event_type>::iterator circle_event_type_it;
+ typedef typename std::list<circle_event_type>::iterator
+ circle_event_type_it;
circle_events_queue() {}
- void reset() {
+ void clear() {
while (!circle_events_.empty())
circle_events_.pop();
circle_events_list_.clear();
@@ -463,8 +567,11 @@
}
};
+ // Remove all the inactive events from the top of the priority
+ // queue until the first active event is found.
void remove_not_active_events() {
- while (!circle_events_.empty() && !circle_events_.top()->is_active())
+ while (!circle_events_.empty() &&
+ !circle_events_.top()->is_active())
pop();
}
@@ -481,138 +588,192 @@
///////////////////////////////////////////////////////////////////////////
// GEOMETRY PREDICATES ////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////
- static void avoid_cancellation(double value, double &left_expr, double &right_expr) {
+
+ // Geometry predicates with floating-point variables usually require
+ // high-precision predicates to retrieve the correct result.
+ // Epsilon robust predicates give the result within some epsilon relative
+ // error, but are a lot faster than high-precision predicates.
+ // To make algorithm robust and efficient epsilon robust predicates are
+ // used at the first step. In case of the undefined result high-precision
+ // arithmetic is used to produce required robustness. This approach
+ // requires exact computation of epsilon intervals within which epsilon
+ // robust predicates have undefined value.
+ // There are two ways to measure an error of floating-point calculations:
+ // relative error and ULPs(units in the last place).
+ // Let EPS be machine epsilon, then next inequalities have place:
+ // 1 EPS <= 1 ULP <= 2 EPS (1), 0.5 ULP <= 1 EPS <= 1 ULP (2).
+ // ULPs are good for measuring rounding errors and comparing values.
+ // Relative erros are good for computation of general relative
+ // errors of formulas or expressions. So to calculate epsilon
+ // intervals within which epsilon robust predicates have undefined result
+ // next schema is used:
+ // 1) Compute rounding errors of initial variables using ULPs;
+ // 2) Transform ULPs to epsilons using upper bound of the (1);
+ // 3) Compute relative error of the formula using epsilon arithmetic;
+ // 4) Transform epsilon to ULPs using upper bound of the (2);
+ // In case two values are inside undefined ULP range use high-precision
+ // arithmetic to produce the correct result, else output the result.
+ // Look at almost_equal function to see how two floating-point variables
+ // are checked to fit in the ULP range.
+ // If A has relative error of r(A) and B has relative error of r(B) then:
+ // 1) r(A + B) <= max(r(A), r(B)), for A * B >= 0;
+ // 2) r(A - B) <= B*r(A)+A*r(B)/(A-B), for A * B >= 0;
+ // 2) r(A * B) <= r(A) + r(B);
+ // 3) r(A / B) <= r(A) + r(B);
+ // In addition rounding error should be added, that is always equal to
+ // 0.5 ULP or atmost 1 epsilon. As you might see from the above formulas
+ // substraction relative error may be extremely large, that's why
+ // epsilon robust comparator class is used to store floating point values
+ // and avoid substraction.
+ // For further information about relative errors and ULPs try this link:
+ // http://docs.sun.com/source/806-3568/ncg_goldberg.html
+
+ // Used during epsilon robust computations. Avoids substraction.
+ static inline void avoid_cancellation(double value, double &left_expr,
+ double &right_expr) {
if (value >= 0)
left_expr += value;
else
right_expr -= value;
}
- // If two floating-point numbers in the same format are ordered (x < y), then
- // they are ordered the same way when their bits are reinterpreted as
- // sign-magnitude integers.
- static bool almost_equal(double a, double b, int maxUlps) {
+ // Convert value to 64-bit unsigned integer.
+ // Return true if the value is positive, else false.
+ template <typename T>
+ static inline bool convert_to_65_bit(T value, unsigned long long &res) {
+ if (value >= 0) {
+ res = static_cast<unsigned long long>(value);
+ return true;
+ } else {
+ res = static_cast<unsigned long long>(-value);
+ return false;
+ }
+ }
+
+ template <typename T>
+ static inline bool compute_dif_65_bit(T value1, T value2,
+ unsigned long long &res) {
+ if (value1 >= value2) {
+ res = static_cast<unsigned long long>(value1) -
+ static_cast<unsigned long long>(value2);
+ return true;
+ } else {
+ res = static_cast<unsigned long long>(value2) -
+ static_cast<unsigned long long>(value1);
+ return false;
+ }
+ }
+
+ // If two floating-point numbers in the same format are ordered (x < y),
+ // then they are ordered the same way when their bits are reinterpreted as
+ // sign-magnitude integers. Values are considered to be almost equal if
+ // their integer reinterpretatoins differ in not more than maxUlps units.
+ static inline bool almost_equal(double a, double b, int maxUlps) {
long long ll_a, ll_b;
- // Reinterpret double bits as long long.
+
+ // Reinterpret double bits as 64-bit signed integer.
memcpy(&ll_a, &a, sizeof(double));
memcpy(&ll_b, &b, sizeof(double));
// Positive 0.0 is integer zero. Negative 0.0 is 0x8000000000000000.
// Map negative zero to an integer zero representation - making it
- // identical to positive zero - and it makes it so that the smallest
- // negative number is represented by negative one, and downwards from there.
+ // identical to positive zero - the smallest negative number is
+ // represented by negative one, and downwards from there.
if (ll_a < 0)
ll_a = 0x8000000000000000LL - ll_a;
if (ll_b < 0)
ll_b = 0x8000000000000000LL - ll_b;
- // Compare long long representations of input values.
+ // Compare 64-bit signed integer representations of input values.
// Difference in 1 Ulp is equivalent to a relative error of between
// 1/4,000,000,000,000,000 and 1/8,000,000,000,000,000.
long long dif = ll_a - ll_b;
return (dif <= maxUlps) && (dif >= -maxUlps);
}
- // TODO(asydorchuk): Make templates specification for integer coordinate type,
- // as it is actually working with integer input.
+ // Robust orientation test. Works correctly for any input type that
+ // can be casted without lose of data to the integer type within a range
+ // [-2^32, 2^32-1].
+ // Arguments: dif_x1_, dif_y1 - coordinates of the first vector.
+ // dif_x2_, dif_y2 - coordinates of the second vector.
+ // Returns orientation test result for input vectors.
template <typename T>
- static kOrientation orientation_test(const point_2d<T> &point1,
- const point_2d<T> &point2,
- const point_2d<T> &point3) {
- typedef long long ll;
+ static kOrientation orientation_test(T dif_x1_, T dif_y1_,
+ T dif_x2_, T dif_y2_) {
typedef unsigned long long ull;
- ull dif_x1, dif_x2, dif_y1, dif_y2;
+ ull dif_x1, dif_y1, dif_x2, dif_y2;
bool dif_x1_plus, dif_x2_plus, dif_y1_plus, dif_y2_plus;
- INT_PREDICATE_COMPUTE_DIFFERENCE(static_cast<ll>(point1.x()),
- static_cast<ll>(point2.x()),
- dif_x1, dif_x1_plus);
- INT_PREDICATE_COMPUTE_DIFFERENCE(static_cast<ll>(point2.x()),
- static_cast<ll>(point3.x()),
- dif_x2, dif_x2_plus);
- INT_PREDICATE_COMPUTE_DIFFERENCE(static_cast<ll>(point1.y()),
- static_cast<ll>(point2.y()),
- dif_y1, dif_y1_plus);
- INT_PREDICATE_COMPUTE_DIFFERENCE(static_cast<ll>(point2.y()),
- static_cast<ll>(point3.y()),
- dif_y2, dif_y2_plus);
+ dif_x1_plus = convert_to_65_bit(dif_x1_, dif_x1);
+ dif_y1_plus = convert_to_65_bit(dif_y1_, dif_y1);
+ dif_x2_plus = convert_to_65_bit(dif_x2_, dif_x2);
+ dif_y2_plus = convert_to_65_bit(dif_y2_, dif_y2);
+
ull expr_l = dif_x1 * dif_y2;
- bool expr_l_plus = (dif_x1_plus == dif_y2_plus) ? true : false;
ull expr_r = dif_x2 * dif_y1;
+
+ bool expr_l_plus = (dif_x1_plus == dif_y2_plus) ? true : false;
bool expr_r_plus = (dif_x2_plus == dif_y1_plus) ? true : false;
if (expr_l == 0)
expr_l_plus = true;
if (expr_r == 0)
expr_r_plus = true;
-
+
if ((expr_l_plus == expr_r_plus) && (expr_l == expr_r))
- return COLINEAR;
+ return COLLINEAR;
if (!expr_l_plus) {
if (expr_r_plus)
return RIGHT_ORIENTATION;
else
- return (expr_l > expr_r) ? RIGHT_ORIENTATION : LEFT_ORIENTATION;
+ return (expr_l > expr_r) ?
+ RIGHT_ORIENTATION : LEFT_ORIENTATION;
} else {
if (!expr_r_plus)
return LEFT_ORIENTATION;
else
- return (expr_l < expr_r) ? RIGHT_ORIENTATION : LEFT_ORIENTATION;
+ return (expr_l < expr_r) ?
+ RIGHT_ORIENTATION : LEFT_ORIENTATION;
}
}
+ // Robust orientation test. Works correctly for any input coordinate type
+ // that can be casted without lose of data to integer type within a range
+ // [-2^31, 2^31 - 1] - signed integer type.
+ // Arguments: point1, point2 - represent the first vector;
+ // point2, point3 - represent the second vector;
+ // Returns orientation test result for input vectors.
template <typename T>
- static kOrientation orientation_test(T dif_x1_, T dif_y1_, T dif_x2_, T dif_y2_) {
- typedef unsigned long long ull;
- ull dif_x1, dif_y1, dif_x2, dif_y2;
- bool dif_x1_plus, dif_x2_plus, dif_y1_plus, dif_y2_plus;
- INT_PREDICATE_CONVERT_65_BIT(dif_x1_, dif_x1, dif_x1_plus);
- INT_PREDICATE_CONVERT_65_BIT(dif_y1_, dif_y1, dif_y1_plus);
- INT_PREDICATE_CONVERT_65_BIT(dif_x2_, dif_x2, dif_x2_plus);
- INT_PREDICATE_CONVERT_65_BIT(dif_y2_, dif_y2, dif_y2_plus);
-
- ull expr_l = dif_x1 * dif_y2;
- bool expr_l_plus = (dif_x1_plus == dif_y2_plus) ? true : false;
- ull expr_r = dif_x2 * dif_y1;
- bool expr_r_plus = (dif_x2_plus == dif_y1_plus) ? true : false;
-
- if (expr_l == 0)
- expr_l_plus = true;
- if (expr_r == 0)
- expr_r_plus = true;
-
- if ((expr_l_plus == expr_r_plus) && (expr_l == expr_r))
- return COLINEAR;
-
- if (!expr_l_plus) {
- if (expr_r_plus)
- return RIGHT_ORIENTATION;
- else
- return (expr_l > expr_r) ? RIGHT_ORIENTATION : LEFT_ORIENTATION;
- } else {
- if (!expr_r_plus)
- return LEFT_ORIENTATION;
- else
- return (expr_l < expr_r) ? RIGHT_ORIENTATION : LEFT_ORIENTATION;
- }
+ static inline kOrientation orientation_test(const point_2d<T> &point1,
+ const point_2d<T> &point2,
+ const point_2d<T> &point3) {
+ return orientation_test(point1.x() - point2.x(),
+ point1.y() - point2.y(),
+ point2.x() - point3.x(),
+ point2.y() - point3.y());
}
+ // Value is a determinant of two vectors.
+ // Return orientation based on the sign of the determinant.
template <typename T>
- static kOrientation orientation_test(T value) {
+ static inline kOrientation orientation_test(T value) {
if (value == static_cast<T>(0.0))
- return COLINEAR;
- return (value < static_cast<T>(0.0)) ? RIGHT_ORIENTATION : LEFT_ORIENTATION;
+ return COLLINEAR;
+ return (value < static_cast<T>(0.0)) ?
+ RIGHT_ORIENTATION : LEFT_ORIENTATION;
}
+ // Compute robust cross_product: a1 * b2 - b1 * a2.
+ // The result is correct with epsilon relative error equal to 1EPS.
template <typename T>
static double robust_cross_product(T a1_, T b1_, T a2_, T b2_) {
typedef unsigned long long ull;
ull a1, b1, a2, b2;
bool a1_plus, a2_plus, b1_plus, b2_plus;
- INT_PREDICATE_CONVERT_65_BIT(a1_, a1, a1_plus);
- INT_PREDICATE_CONVERT_65_BIT(b1_, b1, b1_plus);
- INT_PREDICATE_CONVERT_65_BIT(a2_, a2, a2_plus);
- INT_PREDICATE_CONVERT_65_BIT(b2_, b2, b2_plus);
+ a1_plus = convert_to_65_bit(a1_, a1);
+ b1_plus = convert_to_65_bit(b1_, b1);
+ a2_plus = convert_to_65_bit(a2_, a2);
+ b2_plus = convert_to_65_bit(b2_, b2);
ull expr_l = a1 * b2;
bool expr_l_plus = (a1_plus == b2_plus) ? true : false;
@@ -623,33 +784,56 @@
expr_l_plus = true;
if (expr_r == 0)
expr_r_plus = true;
-
+
if ((expr_l_plus == expr_r_plus) && (expr_l == expr_r))
return static_cast<T>(0.0);
+ // Produce the result with epsilon relative error.
if (!expr_l_plus) {
if (expr_r_plus)
- return -static_cast<double>(expr_l) - static_cast<double>(expr_r);
- else
- return (expr_l > expr_r) ? -static_cast<double>(expr_l - expr_r) :
- static_cast<double>(expr_r - expr_l);
+ return -static_cast<double>(expr_l) -
+ static_cast<double>(expr_r);
+ else return (expr_l > expr_r) ?
+ -static_cast<double>(expr_l - expr_r) :
+ static_cast<double>(expr_r - expr_l);
} else {
if (!expr_r_plus)
- return static_cast<double>(expr_l) + static_cast<double>(expr_r);
+ return static_cast<double>(expr_l) +
+ static_cast<double>(expr_r);
else
- return (expr_l < expr_r) ? -static_cast<double>(expr_r - expr_l) :
- static_cast<double>(expr_l - expr_r);
+ return (expr_l < expr_r) ?
+ -static_cast<double>(expr_r - expr_l) :
+ static_cast<double>(expr_l - expr_r);
}
}
+ // Class used to make computations with epsilon relative error.
+ // ERC consists of two values: value1 and value2, both positive.
+ // The resulting expression is equal to the value1 - value2.
+ // The main idea is to represent any expression that consists of
+ // addition, substraction, multiplication and division operations
+ // to avoid using substraction. Substraction of a positive value
+ // is equivalent to the addition to value2 and substraction of
+ // a negative value is equivalent to the addition to value1.
+ // Cons: ERC gives error relative not to the resulting value,
+ // but relative to some expression instead. Example:
+ // center_x = 100, ERC's value1 = 10^20, value2 = 10^20,
+ // center_x = 1000, ERC's value3 = 10^21, value4 = 10^21,
+ // such two centers are considered equal(
+ // value1 + value4 = value2 + value3), while they are not.
+ // Pros: ERC is much faster then approaches based on the use
+ // of high-precision libraries. However this will give correct
+ // answer for the previous example.
+ // Solution: Use ERCs in case of defined comparison results and use
+ // high-precision libraries for undefined results.
template <typename T>
class epsilon_robust_comparator {
public:
- epsilon_robust_comparator() :
+ epsilon_robust_comparator() :
positive_sum_(0),
negative_sum_(0) {}
- epsilon_robust_comparator(T value) :
+ epsilon_robust_comparator(T value) :
positive_sum_((value>0)?value:0),
negative_sum_((value<0)?-value:0) {}
@@ -669,6 +853,7 @@
return negative_sum_;
}
+ // Equivalent to the unary minus.
void swap() {
(std::swap)(positive_sum_, negative_sum_);
}
@@ -681,29 +866,31 @@
return false;
}
- epsilon_robust_comparator<T>& operator+=(const T& val) {
+ epsilon_robust_comparator<T> &operator+=(const T &val) {
avoid_cancellation(val, positive_sum_, negative_sum_);
return *this;
}
- epsilon_robust_comparator<T>& operator+=(const epsilon_robust_comparator<T>& erc) {
+ epsilon_robust_comparator<T> &operator+=(
+ const epsilon_robust_comparator<T> &erc) {
positive_sum_ += erc.positive_sum_;
negative_sum_ += erc.negative_sum_;
return *this;
}
- epsilon_robust_comparator<T>& operator-=(const T& val) {
+ epsilon_robust_comparator<T> &operator-=(const T &val) {
avoid_cancellation(val, negative_sum_, positive_sum_);
return *this;
}
- epsilon_robust_comparator<T>& operator-=(const epsilon_robust_comparator<T>& erc) {
+ epsilon_robust_comparator<T> &operator-=(
+ const epsilon_robust_comparator<T> &erc) {
positive_sum_ += erc.negative_sum_;
negative_sum_ += erc.positive_sum_;
return *this;
}
- epsilon_robust_comparator<T>& operator*=(const T& val) {
+ epsilon_robust_comparator<T> &operator*=(const T &val) {
positive_sum_ *= fabs(val);
negative_sum_ *= fabs(val);
if (val < 0) {
@@ -712,7 +899,7 @@
return *this;
}
- epsilon_robust_comparator<T>& operator/=(const T& val) {
+ epsilon_robust_comparator<T> &operator/=(const T &val) {
positive_sum_ /= fabs(val);
negative_sum_ /= fabs(val);
if (val < 0) {
@@ -721,6 +908,7 @@
return *this;
}
+ // Compare predicate with value using ulp precision.
kPredicateResult compare(T value, int ulp = 0) const {
T lhs = positive_sum_ - ((value < 0) ? value : 0);
T rhs = negative_sum_ + ((value > 0) ? value : 0);
@@ -729,7 +917,9 @@
return (lhs < rhs) ? LESS : MORE;
}
- kPredicateResult compare(const epsilon_robust_comparator &rc, int ulp = 0) const {
+ // Compare two predicats using ulp precision.
+ kPredicateResult compare(const epsilon_robust_comparator &rc,
+ int ulp = 0) const {
T lhs = positive_sum_ + rc.neg();
T rhs = negative_sum_ + rc.pos();
if (almost_equal(lhs, rhs, ulp))
@@ -737,7 +927,9 @@
return (lhs < rhs) ? LESS : MORE;
}
- bool compares_undefined(const epsilon_robust_comparator &rc, int ulp) const {
+ // Check whether comparison has undefined result.
+ bool compares_undefined(const epsilon_robust_comparator &rc,
+ int ulp) const {
return this->compare(rc, ulp) == UNDEFINED;
}
@@ -747,116 +939,179 @@
};
template<typename T>
- epsilon_robust_comparator<T> operator+(const epsilon_robust_comparator<T>& lhs,
- const epsilon_robust_comparator<T>& rhs) {
- return epsilon_robust_comparator<T>(lhs.pos() + rhs.pos(), lhs.neg() + rhs.neg());
+ inline epsilon_robust_comparator<T> operator+(
+ const epsilon_robust_comparator<T>& lhs,
+ const epsilon_robust_comparator<T>& rhs) {
+ return epsilon_robust_comparator<T>(lhs.pos() + rhs.pos(),
+ lhs.neg() + rhs.neg());
}
template<typename T>
- epsilon_robust_comparator<T> operator-(const epsilon_robust_comparator<T>& lhs,
- const epsilon_robust_comparator<T>& rhs) {
- return epsilon_robust_comparator<T>(lhs.pos() - rhs.neg(), lhs.neg() + rhs.pos());
+ inline epsilon_robust_comparator<T> operator-(
+ const epsilon_robust_comparator<T>& lhs,
+ const epsilon_robust_comparator<T>& rhs) {
+ return epsilon_robust_comparator<T>(lhs.pos() - rhs.neg(),
+ lhs.neg() + rhs.pos());
}
template<typename T>
- epsilon_robust_comparator<T> operator*(const epsilon_robust_comparator<T>& lhs,
- const epsilon_robust_comparator<T>& rhs) {
+ inline epsilon_robust_comparator<T> operator*(
+ const epsilon_robust_comparator<T>& lhs,
+ const epsilon_robust_comparator<T>& rhs) {
double res_pos = lhs.pos() * rhs.pos() + lhs.neg() * rhs.neg();
double res_neg = lhs.pos() * rhs.neg() + lhs.neg() * rhs.pos();
return epsilon_robust_comparator<T>(res_pos, res_neg);
}
template<typename T>
- epsilon_robust_comparator<T> operator*(const epsilon_robust_comparator<T>& lhs,
- const T& val) {
+ inline epsilon_robust_comparator<T> operator*(
+ const epsilon_robust_comparator<T>& lhs, const T& val) {
if (val >= 0)
- return epsilon_robust_comparator<T>(lhs.pos() * val, lhs.neg() * val);
- return epsilon_robust_comparator<T>(-lhs.neg() * val, -lhs.pos() * val);
+ return epsilon_robust_comparator<T>(lhs.pos() * val,
+ lhs.neg() * val);
+ return epsilon_robust_comparator<T>(-lhs.neg() * val,
+ -lhs.pos() * val);
}
template<typename T>
- epsilon_robust_comparator<T> operator*(const T& val,
- const epsilon_robust_comparator<T>& rhs) {
+ inline epsilon_robust_comparator<T> operator*(
+ const T& val, const epsilon_robust_comparator<T>& rhs) {
if (val >= 0)
- return epsilon_robust_comparator<T>(val * rhs.pos(), val * rhs.neg());
- return epsilon_robust_comparator<T>(-val * rhs.neg(), -val * rhs.pos());
+ return epsilon_robust_comparator<T>(val * rhs.pos(),
+ val * rhs.neg());
+ return epsilon_robust_comparator<T>(-val * rhs.neg(),
+ -val * rhs.pos());
}
- // Returns true if horizontal line going through new site intersects
- // right arc at first, else returns false. If horizontal line goes
- // through intersection point of the given two arcs returns false also.
- // Used during nodes comparison.
- // Let x0 be sweepline coordinate, left site coordinates be (x1, y1),
+ // Robust voronoi vertex data structure. Used during removing degenerate
+ // edges(zero-length).
+ // Vertex coordinates are: center_x / denom, center_y / denom.
+ // Variables: center_x - x-coordinate of the circle's center;
+ // center_y - y-coordinate of the circle's center;
+ // denom - denominator for the previous two values.
+ template <typename T>
+ class robust_voronoi_vertex {
+ public:
+ typedef T coordinate_type;
+ typedef epsilon_robust_comparator<coordinate_type> erc_type;
+
+ robust_voronoi_vertex(const erc_type &c_x,
+ const erc_type &c_y,
+ const erc_type &den) :
+ center_x(c_x),
+ center_y(c_y),
+ denom(den) {}
+
+ coordinate_type x() const { return center_x.dif() / denom.dif(); }
+
+ coordinate_type y() const { return center_y.dif() / denom.dif(); }
+
+ // Compare two vertices with the given ulp precision.
+ bool equals(const robust_voronoi_vertex *that, int ulp) const {
+ erc_type lhs1(this->center_x * that->denom);
+ erc_type rhs1(this->denom * that->center_x);
+ erc_type lhs2(this->center_y * that->denom);
+ erc_type rhs2(this->denom * that->center_y);
+ return lhs1.compares_undefined(rhs1, ulp) &&
+ lhs2.compares_undefined(rhs2, ulp);
+ }
+
+ private:
+ erc_type center_x;
+ erc_type center_y;
+ erc_type denom;
+ };
+
+ // Epsilon robust predicate.
+ // Let x0 be the sweepline coordinate, left site coordinates be (x1, y1),
// right site coordinates be (x2, y2). Equations of the arcs will be:
// x1(y) = ((y - y1)^2 + x1^2 - x0^2) / (2*(x1 - x0));
// x2(y) = ((y - y2)^2 + x2^2 - x0^2) / (2*(x2 - x0)).
- // Horizontal line going throught site (x*, y*) intersects second arc
+ // Horizontal line going throught site (x0, y0) intersects second arc
// at first if x2(y*) > x1(y*) or:
- // (x0-x2)*(x0-x1)*(x1-x2) + (x0-x2)*(y*-y1)^2 < (x0-x1)*(y*-y2)^2.
- template <typename T>
- static kPredicateResult fast_less_predicate(const point_2d<T> &left_point,
- const point_2d<T> &right_point,
- const point_2d<T> &new_point) {
- double fast_a1 = static_cast<double>(new_point.x()) - static_cast<double>(left_point.x());
- double fast_a2 = static_cast<double>(new_point.x()) - static_cast<double>(right_point.x());
- double fast_b1 = static_cast<double>(new_point.y()) - static_cast<double>(left_point.y());
- double fast_b2 = static_cast<double>(new_point.y()) - static_cast<double>(right_point.y());
- double fast_c = static_cast<double>(left_point.x()) - static_cast<double>(right_point.x());
+ // (x0-x2)*(x0-x1)*(x1-x2) + (x0-x1)*(y0-y2)^2 < (x0-x2)*(y0-y1)^2.
+ // Works correctly for any input coordinate type that can be casted without
+ // lose of data to the integer type within a range [-2^31, 2^31-1].
+ template <typename T>
+ static kPredicateResult eps_less_predicate(const point_2d<T> &left_point,
+ const point_2d<T> &right_point,
+ const point_2d<T> &new_point) {
+ // Compute x0 - x1.
+ double fast_a1 = new_point.x() - left_point.x();
+
+ // Compute x0 - x2.
+ double fast_a2 = new_point.x() - right_point.x();
+
+ // Compute y0 - y1.
+ double fast_b1 = new_point.y() - left_point.y();
+
+ // Compute y0 - y2.
+ double fast_b2 = new_point.y() - right_point.y();
+
+ // Compute x1 - x2.
+ double fast_c = left_point.x() - right_point.x();
+
+ // Compute (x0-x2)*(y0-y1)^2. This value is always positive(x0 > x2).
+ // Relative error is equal to 2*EPS.
double fast_left_expr = fast_a1 * fast_b2 * fast_b2;
+
+ // Compute (x0-x1)*(y0-y2)^2. This value is always positive(x0 > x1).
+ // Relative error is equal to 2*EPS.
double fast_right_expr = fast_a2 * fast_b1 * fast_b1;
-
- // Avoid cancellation.
- avoid_cancellation(fast_a1 * fast_a2 * fast_c, fast_left_expr, fast_right_expr);
+
+ // Avoid substraction while adding (x0-x2)*(x0-x1)*(x1-x2).
+ // In case its value is positive add it to the left expression,
+ // else add its module to the right expression.
+ // Relative error of the (x0-x2)*(x0-x1)*(x1-x2) is 2*EPS.
+ // Relative error of the addition to any of the expression is:
+ // max(2*EPS, 2*EPS) + EPS = 3*EPS.
+ avoid_cancellation(fast_a1 * fast_a2 * fast_c,
+ fast_left_expr, fast_right_expr);
+
+ // If expression are outside undefined ULP range, return exact result.
+ // Relative error of one expression is 2*EPS and 3*EPS for another.
+ // 2*EPS + 3*EPS = 5*EPS <= 5*ULP. That's why predicate may produce
+ // exact result in case expressions are out of the 5*ULP range.
if (!almost_equal(fast_left_expr, fast_right_expr, 5))
return (fast_left_expr < fast_right_expr) ? LESS : MORE;
+
+ // Expressions are inside of the 5*ULP range. Result is undefined.
return UNDEFINED;
}
- template <typename T>
- static bool less_predicate(const point_2d<T> &left_point,
- const point_2d<T> &right_point,
- const point_2d<T> &new_point) {
- if (left_point.x() > right_point.x()) {
- if (new_point.y() <= left_point.y())
- return false;
- } else if (left_point.x() < right_point.x()) {
- if (new_point.y() >= right_point.y())
- return true;
- } else {
- return left_point.y() + right_point.y() < 2.0 * new_point.y();
- }
-
- kPredicateResult fast_res = fast_less_predicate(left_point, right_point, new_point);
- if (fast_res != UNDEFINED)
- return (fast_res == LESS);
-
+ // Robust 65-bit predicate, avoids using high-precision libraries.
+ // Works correctly for any input coordinate type that can be casted without
+ // lose of data to the integer type within a range [-2^31, 2^31 - 1].
+ // (x0-x2)*(x0-x1)*(x1-x2) + (x0-x1)*(y0-y2)^2 < (x0-x2)*(y0-y1)^2 check
+ // is equal to: (x1-x2) + (y0-y2)^2/(x0-x2) < (y0-y1)^2 / (x0-x1).
+ template <typename T>
+ static bool robust_65bit_less_predicate(const point_2d<T> &left_point,
+ const point_2d<T> &right_point,
+ const point_2d<T> &new_point) {
typedef long long ll;
typedef unsigned long long ull;
ull a1, a2, b1, b2, b1_sqr, b2_sqr, l_expr, r_expr;
bool l_expr_plus, r_expr_plus;
- // a1 and a2 are greater than zero.
+ // Compute a1 = (x0-x1), a2 = (x0-x2), a1 and a2 are positive.
a1 = static_cast<ull>(static_cast<ll>(new_point.x()) -
static_cast<ll>(left_point.x()));
a2 = static_cast<ull>(static_cast<ll>(new_point.x()) -
static_cast<ll>(right_point.x()));
- // We don't need to know signs of b1 and b2, because we use their squared values.
- INT_PREDICATE_COMPUTE_DIFFERENCE(static_cast<ll>(new_point.y()),
- static_cast<ll>(left_point.y()),
- b1, l_expr_plus);
- INT_PREDICATE_COMPUTE_DIFFERENCE(static_cast<ll>(new_point.y()),
- static_cast<ll>(right_point.y()),
- b2, l_expr_plus);
+ // Compute b1_sqr = (y0-y1)^2 and b2_sqr = (y0-y2)^2. We don't need
+ // to know signs of b1 and b2, because we use their squared values.
+ l_expr_plus = compute_dif_65_bit(new_point.y(), left_point.y(), b1);
+ l_expr_plus = compute_dif_65_bit(new_point.y(), right_point.y(), b2);
b1_sqr = b1 * b1;
b2_sqr = b2 * b2;
ull b1_sqr_mod = b1_sqr % a1;
ull b2_sqr_mod = b2_sqr % a2;
- // Compute left expression.
- INT_PREDICATE_COMPUTE_DIFFERENCE(static_cast<ll>(left_point.x()),
- static_cast<ll>(right_point.x()),
- l_expr, l_expr_plus);
+ // Compute l_expr = (x1 - x2).
+ l_expr_plus = compute_dif_65_bit(left_point.x(), right_point.x(), l_expr);
+
+ // In case (b2_sqr / a2) < (b1_sqr / a1), decrease left_expr by 1.
if (b2_sqr_mod * a1 < b1_sqr_mod * a2) {
if (!l_expr_plus)
l_expr++;
@@ -869,7 +1124,7 @@
}
// Compute right expression.
- INT_PREDICATE_COMPUTE_DIFFERENCE(b1_sqr / a1, b2_sqr / a2, r_expr, r_expr_plus);
+ r_expr_plus = compute_dif_65_bit(b1_sqr / a1, b2_sqr / a2, r_expr);
// Compare left and right expressions.
if (!l_expr_plus && r_expr_plus)
@@ -881,18 +1136,56 @@
return l_expr > r_expr;
}
+ // Robust predicate, avoids using high-precision libraries.
+ // Returns true if horizontal line going through the new point site
+ // intersects right arc at first, else returns false. If horizontal line
+ // goes through intersection point of the given two arcs returns false.
+ // Works correctly for any input coordinate type that can be casted without
+ // lose of data to the integer type within a range [-2^31, 2^31-1].
+ template <typename T>
+ static bool less_predicate(const point_2d<T> &left_point,
+ const point_2d<T> &right_point,
+ const point_2d<T> &new_point) {
+ // Any two point sites with different x-coordinates create two
+ // bisectors. Each bisector is defined by the order the sites
+ // appear in its representation. Predicates used in this function
+ // won't produce the correct result for the any arrangment of the
+ // input sites. That's why some preprocessing is required to handle
+ // such cases.
+ if (left_point.x() > right_point.x()) {
+ if (new_point.y() <= left_point.y())
+ return false;
+ } else if (left_point.x() < right_point.x()) {
+ if (new_point.y() >= right_point.y())
+ return true;
+ } else {
+ // If x-coordinates of the sites are equal, we may produce the
+ // result without any further computations.
+ return left_point.y() + right_point.y() < 2.0 * new_point.y();
+ }
+
+ // Try epsilon robust predicate at first.
+ kPredicateResult eps_res = eps_less_predicate(left_point, right_point, new_point);
+ if (eps_res != UNDEFINED)
+ return (eps_res == LESS);
+
+ // If the result of the epsilon robust predicate is undefined
+ // use 65-bit robust predicate that produces exact result.
+ return robust_65bit_less_predicate(left_point, right_point, new_point);
+ }
+
template <typename T>
static kPredicateResult fast_less_predicate(point_2d<T> site_point, site_event<T> segment_site,
point_2d<T> new_point, bool reverse_order) {
typedef long long ll;
typedef unsigned long long ull;
- if (orientation_test(segment_site.get_point0(true), segment_site.get_point1(true),
+ if (orientation_test(segment_site.point0(true), segment_site.point1(true),
new_point) != RIGHT_ORIENTATION) {
return (!segment_site.is_inverse()) ? LESS : MORE;
}
- const point_2d<T> &segment_start = segment_site.get_point0();
- const point_2d<T> &segment_end = segment_site.get_point1();
+ const point_2d<T> &segment_start = segment_site.point0();
+ const point_2d<T> &segment_end = segment_site.point1();
ll dif_x = static_cast<ll>(new_point.x()) - static_cast<ll>(site_point.x());
ll dif_y = static_cast<ll>(new_point.y()) - static_cast<ll>(site_point.y());
ll a = static_cast<ll>(segment_end.x()) - static_cast<ll>(segment_start.x());
@@ -906,15 +1199,14 @@
return UNDEFINED;
} else {
kOrientation orientation = orientation_test(a, b, dif_x, dif_y);
- if ((orientation == COLINEAR) ||
+ if ((orientation == COLLINEAR) ||
(!segment_site.is_inverse() ^ (orientation == RIGHT_ORIENTATION))) {
if (!segment_site.is_inverse())
return reverse_order ? LESS : UNDEFINED;
return reverse_order ? UNDEFINED : MORE;
- }
+ }
}
- // dif_x and dif_y are integers, so there will be no cancellation issues.
double fast_left_expr = static_cast<double>(a) *
static_cast<double>(dif_y + dif_x) *
static_cast<double>(dif_y - dif_x);
@@ -928,11 +1220,11 @@
}
ull a_rob, a_sign, b_rob, b_sign, dif_x_rob, dif_x_sign, dif_y_rob, dif_y_sign;
- INT_PREDICATE_CONVERT_65_BIT(a, a_rob, a_sign);
- INT_PREDICATE_CONVERT_65_BIT(b, b_rob, b_sign);
- INT_PREDICATE_CONVERT_65_BIT(dif_x, dif_x_rob, dif_x_sign);
- INT_PREDICATE_CONVERT_65_BIT(dif_y, dif_y_rob, dif_y_sign);
-
+ a_sign = convert_to_65_bit(a, a_rob);
+ b_sign = convert_to_65_bit(b, b_rob);
+ dif_x_sign = convert_to_65_bit(dif_x, dif_x_rob);
+ dif_y_sign = convert_to_65_bit(dif_y, dif_y_rob);
+
ull dif_x_sqr = dif_x_rob * dif_x_rob;
ull dif_y_sqr = dif_y_rob * dif_y_rob;
ull left_expr_div = dif_y_sqr / dif_x_sqr + 1;
@@ -944,7 +1236,7 @@
ull right_expr_mod = right_expr % right_expr_denom;
bool comparison_result;
- if ((b_sign != dif_y_sign) && right_expr_div)
+ if ((b_sign != dif_y_sign) && right_expr_div)
comparison_result = true;
else {
if (b_sign != dif_y_sign && right_expr_mod)
@@ -979,7 +1271,7 @@
}
}
}
-
+
if (segment_site.is_inverse() ^ comparison_result ^ reverse_order)
return reverse_order ? LESS : MORE;
return UNDEFINED;
@@ -1017,14 +1309,11 @@
}
#endif
- // Returns true if horizontal line going through new site intersects
+ // Returns true if a horizontal line going through a new site intersects
// right arc at first, else returns false. If horizontal line goes
- // through intersection point of the given two arcs returns false also.
- // Used during nodes comparison.
- // If reverse order is false we are comparing (point, segment) intersection
- // point and new point, else (segment, point) intersection point.
- // (point, segment) and (segment, point) are two distinct points, except
- // case of vertical segment.
+ // through intersection point of the given two arcs returns false also.
+ // reverse_order flag defines arrangement of the sites. If it's false
+ // the order is (point, segment), else - (segment, point).
template <typename T>
static bool less_predicate(point_2d<T> site_point, site_event<T> segment_site,
point_2d<T> new_point, bool reverse_order) {
@@ -1034,8 +1323,8 @@
return (fast_res == LESS);
}
- const point_2d<T> &segment_start = segment_site.get_point0();
- const point_2d<T> &segment_end = segment_site.get_point1();
+ const point_2d<T> &segment_start = segment_site.point0();
+ const point_2d<T> &segment_end = segment_site.point1();
double dif_x = static_cast<double>(new_point.x()) -
static_cast<double>(site_point.x());
double dif_y = static_cast<double>(new_point.y()) -
@@ -1070,92 +1359,91 @@
#endif
}
- template <typename T>
- static bool less_predicate(site_event<T> left_site,
- site_event<T> right_site,
- point_2d<T> new_point) {
- if (left_site.get_site_index() == right_site.get_site_index()) {
- return orientation_test(left_site.get_point0(), left_site.get_point1(), new_point) ==
- LEFT_ORIENTATION;
- }
-
- if (left_site.get_point1() == new_point && right_site.get_point1() == new_point)
- return false;
-
- const point_2d<T> &segment1_start = left_site.get_point1();
- const point_2d<T> &segment1_end = left_site.get_point0();
- const point_2d<T> &segment2_start = right_site.get_point1();
- const point_2d<T> &segment2_end = right_site.get_point0();
- double intersection_x1 = 0.0;
- double intersection_x2 = 0.0;
-
- double a1 = segment1_end.x() - segment1_start.x();
- if (a1 == 0.0) {
- // Avoid cancellation.
- intersection_x2 += (new_point.x() - segment1_end.x()) * 0.5;
+ // Find the distance between the x-coordinate of the sweepline and the
+ // x-coordinate of the point of the intersection of the horizontal line
+ // going through the new site with the arc corresponding to the segment
+ // site. The relative error is atmost equal to 7EPS.
+ template <typename T>
+ static double find_distance(const site_event<T> &segment,
+ const point_2d<T> &new_point) {
+ if (segment.is_vertical()) {
+ return (segment.point0().x() - new_point.x()) * 0.5;
} else {
- double b1 = segment1_end.y() - segment1_start.y();
- double c1 = b1 * (new_point.x() - segment1_start.x()) + a1 * segment1_start.y();
- double mul1 = sqrt(a1 * a1 + b1 * b1);
- if (left_site.is_inverse()) {
+ const point_2d<T> &segment_start = segment.point1();
+ const point_2d<T> &segment_end = segment.point0();
+ double a1 = segment_end.x() - segment_start.x();
+ double b1 = segment_end.y() - segment_start.y();
+ double a3 = new_point.x() - segment_start.x();
+ double b3 = new_point.y() - segment_start.y();
+ double k = sqrt(a1 * a1 + b1 * b1);
+ // Avoid substraction while computing k.
+ if (segment.is_inverse()) {
if (b1 >= 0.0) {
- mul1 = 1.0 / (b1 + mul1);
+ k = 1.0 / (b1 + k);
} else {
- mul1 = (-b1 + mul1) / (a1 * a1);
+ k = (-b1 + k) / (a1 * a1);
}
} else {
if (b1 >= 0.0) {
- mul1 = (-b1 - mul1) / (a1 * a1);
+ k = (-b1 - k) / (a1 * a1);
} else {
- mul1 = 1.0 / (b1 - mul1);
+ k = 1.0 / (b1 - k);
}
}
- avoid_cancellation(a1 * mul1 * new_point.y(), intersection_x1, intersection_x2);
- avoid_cancellation(-c1 * mul1, intersection_x1, intersection_x2);
+ // Relative error of the robust cross product is 1EPS.
+ // Relative error of the k is atmost 5EPS.
+ // The resulting relative error is atmost 7EPS.
+ return robust_cross_product(a1, b1, a3, b3) * k;
}
+ }
- double a2 = segment2_end.x() - segment2_start.x();
- if (a2 == 0.0) {
- // Avoid cancellation.
- intersection_x1 += (new_point.x() - segment2_end.x()) * 0.5;
- } else {
- double b2 = segment2_end.y() - segment2_start.y();
- double c2 = b2 * (new_point.x() - segment2_start.x()) + a2 * segment2_start.y();
- double mul2 = sqrt(a2 * a2 + b2 * b2);
- if (right_site.is_inverse()) {
- if (b2 >= 0.0) {
- mul2 = 1.0 / (b2 + mul2);
- } else {
- mul2 = (-b2 + mul2) / (a2 * a2);
- }
- } else {
- if (b2 >= 0.0) {
- mul2 = (-b2 - mul2) / (a2 * a2);
- } else {
- mul2 = 1.0 / (b2 - mul2);
- }
- }
- avoid_cancellation(a2 * new_point.y() * mul2, intersection_x2, intersection_x1);
- avoid_cancellation(-c2 * mul2, intersection_x2, intersection_x1);
+ // Returns true if a horizontal line going through a new site intersects
+ // right arc at first, else returns false. If horizontal line goes
+ // through intersection point of the given two arcs returns false also.
+ template <typename T>
+ static bool less_predicate(site_event<T> left_segment,
+ site_event<T> right_segment,
+ point_2d<T> new_point) {
+ // Handle temporary segment sites.
+ if (left_segment.index() == right_segment.index()) {
+ return orientation_test(left_segment.point0(),
+ left_segment.point1(),
+ new_point) == LEFT_ORIENTATION;
+ }
+
+ // Distances between the x-coordinate of the sweepline and
+ // the x-coordinates of the points of the intersections of the
+ // horizontal line going through the new site with arcs corresponding
+ // to the first and to the second segment sites respectively.
+ double dist1 = find_distance(left_segment, new_point);
+ double dist2 = find_distance(right_segment, new_point);
+
+ // The undefined ulp range is equal to 7EPS + 7EPS <= 14ULP.
+ if (!almost_equal(dist1, dist2, 14)) {
+ return dist1 < dist2;
}
- if (!almost_equal(intersection_x1, intersection_x2, 20)) {
- return intersection_x1 < intersection_x2;
- }
-
// TODO(asydorchuk): Add mpl support there.
- return intersection_x1 < intersection_x2;
+ return dist1 < dist2;
}
///////////////////////////////////////////////////////////////////////////
// CIRCLE EVENTS CREATION /////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////
+
+ // At the moment all the circle event creation functions use the
+ // epsilon_robust_comparator class to output the parameters of the
+ // inscribed circle. Such approach may produce incorrect results while
+ // comparing two circle events. In such cases high-precision functions
+ // are required to recalculate circle's parameters.
+
template <typename T>
- static T get_sqr_distance(T dif_x, T dif_y) {
+ static inline T sqr_distance(T dif_x, T dif_y) {
return dif_x * dif_x + dif_y * dif_y;
}
- // Create circle event from three point sites.
+ // Find parameters of the inscribed circle that is tangent to three
+ // point sites.
template <typename T>
static bool create_circle_event_ppp(const site_event<T> &site1,
const site_event<T> &site2,
@@ -1187,13 +1475,14 @@
c_x *= inv_orientation;
c_y *= inv_orientation;
epsilon_robust_comparator<T> lower_x(c_x);
- lower_x += sqrt(get_sqr_distance(dif_x1, dif_y1) * get_sqr_distance(dif_x2, dif_y2) *
- get_sqr_distance(dif_x3, dif_y3)) * fabs(inv_orientation);
+ lower_x += sqrt(sqr_distance(dif_x1, dif_y1) * sqr_distance(dif_x2, dif_y2) *
+ sqr_distance(dif_x3, dif_y3)) * fabs(inv_orientation);
c_event = circle_event<double>(c_x, c_y, lower_x);
return true;
}
- // Create circle event from two point sites and one segment site.
+ // Find parameters of the inscribed circle that is tangent to two
+ // point sites and on segment site.
template <typename T>
static bool create_circle_event_pps(const site_event<T> &site1,
const site_event<T> &site2,
@@ -1201,10 +1490,10 @@
int segment_index,
circle_event<T> &c_event) {
if (segment_index != 2) {
- kOrientation orient1 = orientation_test(site1.get_point0(),
- site2.get_point0(), site3.get_point0(true));
- kOrientation orient2 = orientation_test(site1.get_point0(),
- site2.get_point0(), site3.get_point1(true));
+ kOrientation orient1 = orientation_test(site1.point0(),
+ site2.point0(), site3.point0(true));
+ kOrientation orient2 = orientation_test(site1.point0(),
+ site2.point0(), site3.point1(true));
if (segment_index == 1 && site1.x0() >= site2.x0()) {
if (orient1 != RIGHT_ORIENTATION)
return false;
@@ -1215,26 +1504,26 @@
return false;
}
} else {
- if (site3.get_point0(true) == site1.get_point0() &&
- site3.get_point1(true) == site2.get_point0())
+ if (site3.point0(true) == site1.point0() &&
+ site3.point1(true) == site2.point0())
return false;
}
- double line_a = site3.get_point1().y() - site3.get_point0().y();
- double line_b = site3.get_point0().x() - site3.get_point1().x();
+ double line_a = site3.point1().y() - site3.point0().y();
+ double line_b = site3.point0().x() - site3.point1().x();
double vec_x = site2.y() - site1.y();
double vec_y = site1.x() - site2.x();
double teta = robust_cross_product(line_a, line_b, -vec_y, vec_x);
double A = robust_cross_product(line_a, line_b,
- site3.get_point1().y() - site1.y(),
- site1.x() - site3.get_point1().x());
+ site3.point1().y() - site1.y(),
+ site1.x() - site3.point1().x());
double B = robust_cross_product(line_a, line_b,
- site3.get_point1().y() - site2.y(),
- site2.x() - site3.get_point1().x());
+ site3.point1().y() - site2.y(),
+ site2.x() - site3.point1().x());
double denom = robust_cross_product(vec_x, vec_y, line_a, line_b);
- double inv_segm_len = 1.0 / sqrt(get_sqr_distance(line_a, line_b));
+ double inv_segm_len = 1.0 / sqrt(sqr_distance(line_a, line_b));
epsilon_robust_comparator<double> t;
- if (orientation_test(denom) == COLINEAR) {
+ if (orientation_test(denom) == COLLINEAR) {
t += teta / (4.0 * (A + B));
t -= A * B / (teta * (A + B));
} else {
@@ -1262,27 +1551,28 @@
return true;
}
- // Create circle event from one point site and two segment sites.
+ // Find parameters of the inscribed circle that is tangent to one
+ // point site and two segment sites.
template <typename T>
static bool create_circle_event_pss(const site_event<T> &site1,
const site_event<T> &site2,
const site_event<T> &site3,
int point_index,
circle_event<T> &c_event) {
- // Intersection check.
- if (site2.get_site_index() == site3.get_site_index()) {
+ if (site2.index() == site3.index()) {
return false;
- }
- const point_2d<T> &segm_start1 = site2.get_point1(true);
- const point_2d<T> &segm_end1 = site2.get_point0(true);
- const point_2d<T> &segm_start2 = site3.get_point0(true);
- const point_2d<T> &segm_end2 = site3.get_point1(true);
+ }
+
+ const point_2d<T> &segm_start1 = site2.point1(true);
+ const point_2d<T> &segm_end1 = site2.point0(true);
+ const point_2d<T> &segm_start2 = site3.point0(true);
+ const point_2d<T> &segm_end2 = site3.point1(true);
if (point_index == 2) {
if (!site2.is_inverse() && site3.is_inverse())
return false;
if (site2.is_inverse() == site3.is_inverse() &&
- orientation_test(segm_end1, site1.get_point0(), segm_end2) != RIGHT_ORIENTATION)
+ orientation_test(segm_end1, site1.point0(), segm_end2) != RIGHT_ORIENTATION)
return false;
}
@@ -1291,7 +1581,7 @@
double a2 = segm_end2.x() - segm_start2.x();
double b2 = segm_end2.y() - segm_start2.y();
double orientation = robust_cross_product(b1, a1, b2, a2);
- if (orientation_test(orientation) == COLINEAR) {
+ if (orientation_test(orientation) == COLLINEAR) {
double a = a1 * a1 + b1 * b1;
double c = robust_cross_product(b1, a1, segm_start2.y() - segm_start1.y(),
segm_start2.x() - segm_start1.x());
@@ -1365,19 +1655,19 @@
return true;
}
- // Create circle event from three segment sites.
+ // Find parameters of the inscribed circle that is tangent to three
+ // segment sites.
template <typename T>
static bool create_circle_event_sss(const site_event<T> &site1,
const site_event<T> &site2,
const site_event<T> &site3,
circle_event<T> &c_event) {
- if (site1.get_site_index() == site2.get_site_index() ||
- site2.get_site_index() == site3.get_site_index()) {
+ if (site1.index() == site2.index() || site2.index() == site3.index()) {
return false;
}
double a1 = site1.x1(true) - site1.x0(true);
double b1 = site1.y1(true) - site1.y0(true);
- double c1 = robust_cross_product(b1, a1, site1.y0(true), site1.x0(true));
+ double c1 = robust_cross_product(b1, a1, site1.y0(true), site1.x0(true));
double a2 = site2.x1(true) - site2.x0(true);
double b2 = site2.y1(true) - site2.y0(true);
double c2 = robust_cross_product(b2, a2, site2.y0(true), site2.x0(true));
@@ -1415,78 +1705,78 @@
}
///////////////////////////////////////////////////////////////////////////
- // VORONOI BEACH LINE TYPES ///////////////////////////////////////////////
+ // VORONOI BEACH LINE DATA TYPES //////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////
-
- // Represents bisector made by two arcs that correspond to the left and
- // right sites. Arc is defined as curve with points equidistant from the
- // site and from the sweepline.
- // Let sweepline sweep from left to right and it's current coordinate
- // be x0, site coordinates be (x1, y1). In this case equation of the arc
- // may be written as: (x-x0)^2 = (x-x1)^2 + (y-y1)^2, or
- // x = ((y - y1)^2 + x1^2 - x0^2) / (2*(x1 - x0)).
- // In general case two arcs will create two different bisectors. That's why
- // the order of arcs is important to define unique bisector.
+
+ // Represents a bisector node made by two arcs that correspond to the left
+ // and right sites. Arc is defined as a curve with points equidistant from
+ // the site and from the sweepline. If the site is a point then the arc is
+ // a parabola, otherwise it's a line segment. A segment site event will
+ // produce different bisectors depending on its direction.
+ // In the general case two sites will create two opposite bisectors. That's
+ // why the order of the sites is important to define the unique bisector.
+ // The one site is considered to be newer than the other in case it was
+ // processed by the algorithm later.
template <typename T>
class beach_line_node {
public:
typedef T coordinate_type;
- typedef point_2d<T> Point2D;
+ typedef point_2d<T> point_2d_type;
typedef site_event<T> site_event_type;
beach_line_node() {}
- // Constructs degenerate bisector, used to search arc that is above
- // given site. The input to the constructor is the site event point.
+ // Constructs degenerate bisector, used to search an arc that is above
+ // the given site. The input to the constructor is the new site point.
explicit beach_line_node(const site_event_type &new_point) {
left_site_ = new_point;
right_site_ = new_point;
}
- // Constructs new bisector. The input to the constructor is two sites
- // that create bisector. The order of sites is important.
+ // Constructs a new bisector. The input to the constructor is the two
+ // sites that create the bisector. The order of sites is important.
beach_line_node(const site_event_type &left_point,
const site_event_type &right_point) {
left_site_ = left_point;
right_site_ = right_point;
}
- // Returns the left site of the bisector.
- const site_event_type &get_left_site() const {
+ const site_event_type &left_site() const {
return left_site_;
}
- // Returns the right site of the bisector.
- const site_event_type &get_right_site() const {
+ const site_event_type &right_site() const {
return right_site_;
}
- void set_left_site(const site_event_type &site) {
+ void left_site(const site_event_type &site) {
left_site_ = site;
}
- void set_right_site(const site_event_type &site) {
+ void right_site(const site_event_type &site) {
right_site_ = site;
}
- void set_left_site_inverse() {
- left_site_.set_inverse();
+ void inverse_left_site() {
+ left_site_.inverse();
}
- void set_right_site_inverse() {
- right_site_.set_inverse();
+ void inverse_right_site() {
+ right_site_.inverse();
}
- coordinate_type get_comparison_x() const {
- return (left_site_.get_site_index() >= right_site_.get_site_index()) ?
+ // Get the x coordinate of the newer site.
+ coordinate_type comparison_x() const {
+ return (left_site_.index() >= right_site_.index()) ?
left_site_.x() : right_site_.x();
}
- std::pair<coordinate_type, int> get_comparison_y(bool is_new_node = true) const {
- if (left_site_.get_site_index() == right_site_.get_site_index()) {
+ // Get comparison pair: y coordinate and direction of the newer site.
+ std::pair<coordinate_type, int> comparison_y(bool is_new_node = true) const {
+ if (left_site_.index() == right_site_.index()) {
return std::make_pair(left_site_.y(), 0);
}
- if (left_site_.get_site_index() > right_site_.get_site_index()) {
+ if (left_site_.index() > right_site_.index()) {
if (!is_new_node && left_site_.is_segment() && left_site_.is_vertical()) {
return std::make_pair(left_site_.y1(), 1);
}
@@ -1495,17 +1785,22 @@
return std::make_pair(right_site_.y(), -1);
}
- int get_comparison_index() const {
- return (left_site_.get_site_index() > right_site_.get_site_index()) ?
- left_site_.get_site_index() : right_site_.get_site_index();
+ // Get the index of the newer site.
+ int comparison_index() const {
+ return (left_site_.index() > right_site_.index()) ?
+ left_site_.index() : right_site_.index();
}
- const site_event_type &get_comparison_site() const {
- return (left_site_.get_site_index() >= right_site_.get_site_index()) ?
+ // Get the newer site.
+ const site_event_type &comparison_site() const {
+ return (left_site_.index() >= right_site_.index()) ?
left_site_ : right_site_;
}
- bool less(const Point2D &new_site) const {
+ // Return true if the horizontal line going through the new site
+ // intersects the arc corresponding to the right_site before the
+ // arc corresponding to the left_site.
+ bool less(const point_2d_type &new_site) const {
if (!left_site_.is_segment()) {
return (!right_site_.is_segment()) ? less_pp(new_site) : less_ps(new_site);
} else {
@@ -1513,29 +1808,31 @@
}
}
- bool less_pp(const Point2D &new_site) const {
- return less_predicate(left_site_.get_point0(), right_site_.get_point0(), new_site);
+ private:
+ bool less_pp(const point_2d_type &new_site) const {
+ return less_predicate(left_site_.point0(), right_site_.point0(), new_site);
}
- bool less_ps(const Point2D &new_site) const {
- return less_predicate(left_site_.get_point0(), right_site_, new_site, false);
+ bool less_ps(const point_2d_type &new_site) const {
+ return less_predicate(left_site_.point0(), right_site_, new_site, false);
}
- bool less_sp(const Point2D &new_site) const {
- return less_predicate(right_site_.get_point0(), left_site_, new_site, true);
+ bool less_sp(const point_2d_type &new_site) const {
+ return less_predicate(right_site_.point0(), left_site_, new_site, true);
}
- bool less_ss(const Point2D &new_site) const {
+ bool less_ss(const point_2d_type &new_site) const {
return less_predicate(left_site_, right_site_, new_site);
}
- private:
site_event_type left_site_;
site_event_type right_site_;
};
- // Represents edge data sturcture (bisector segment), which is
- // associated with beach line node key in the binary search tree.
+ // Represents edge data sturcture from the voronoi output, that is
+ // associated as a value with beach line bisector as a key in the beach
+ // line. Contains iterator to the circle event in the circle event
+ // queue in case it's the second bisector that forms a circle event.
template <typename T>
class beach_line_node_data {
public:
@@ -1554,11 +1851,11 @@
contains_circle_event_ = false;
}
- voronoi_edge<T> *get_edge() const {
+ voronoi_edge<T> *edge() const {
return edge_;
}
- void set_edge(voronoi_edge<T> *new_edge) {
+ void edge(voronoi_edge<T> *new_edge) {
edge_ = new_edge;
}
private:
@@ -1567,6 +1864,7 @@
bool contains_circle_event_;
};
+ // Beach line comparison functor.
template <typename BeachLineNode>
struct node_comparer {
public:
@@ -1575,80 +1873,140 @@
// Compares nodes in the balanced binary search tree. Nodes are
// compared based on the y coordinates of the arcs intersection points.
// Nodes with less y coordinate of the intersection point go first.
- // Comparison is only called during site events processing. That's why
- // one of the nodes will always lie on the sweepline. Comparison won't
- // work fine for nodes situated above sweepline.
+ // Comparison is only called during the new site events processing.
+ // That's why one of the nodes will always lie on the sweepline and may
+ // be represented as a straight horizontal line.
bool operator() (const BeachLineNode &node1,
const BeachLineNode &node2) const {
// Get x coordinate of the righmost site from both nodes.
- coordinate_type node1_x = node1.get_comparison_x();
- coordinate_type node2_x = node2.get_comparison_x();
+ coordinate_type node1_x = node1.comparison_x();
+ coordinate_type node2_x = node2.comparison_x();
if (node1_x < node2_x) {
- return node1.less(node2.get_comparison_site().get_point0());
+ // The second node contains a new site.
+ return node1.less(node2.comparison_site().point0());
} else if (node1_x > node2_x) {
- return !node2.less(node1.get_comparison_site().get_point0());
+ // The first node contains a new site.
+ return !node2.less(node1.comparison_site().point0());
} else {
- if (node1.get_comparison_index() == node2.get_comparison_index()) {
- return node1.get_comparison_y() < node2.get_comparison_y();
- } else if (node1.get_comparison_index() < node2.get_comparison_index()) {
- std::pair<coordinate_type, int> y1 = node1.get_comparison_y(false);
- std::pair<coordinate_type, int> y2 = node2.get_comparison_y();
+ // This checks were evaluated experimentally.
+ if (node1.comparison_index() == node2.comparison_index()) {
+ // Both nodes are new (inserted during the same site event
+ // processing).
+ return node1.comparison_y() < node2.comparison_y();
+ } else if (node1.comparison_index() < node2.comparison_index()) {
+ std::pair<coordinate_type, int> y1 = node1.comparison_y(false);
+ std::pair<coordinate_type, int> y2 = node2.comparison_y();
if (y1.first != y2.first) {
return y1.first < y2.first;
}
- return (!node1.get_comparison_site().is_segment()) ? (y1.second < 0) : false;
+ return (!node1.comparison_site().is_segment()) ? (y1.second < 0) : false;
} else {
- std::pair<coordinate_type, int> y1 = node1.get_comparison_y();
- std::pair<coordinate_type, int> y2 = node2.get_comparison_y(false);
+ std::pair<coordinate_type, int> y1 = node1.comparison_y();
+ std::pair<coordinate_type, int> y2 = node2.comparison_y(false);
if (y1.first != y2.first) {
return y1.first < y2.first;
}
- return (!node2.get_comparison_site().is_segment()) ? (y2.second > 0) : true;
+ return (!node2.comparison_site().is_segment()) ? (y2.second > 0) : true;
}
}
}
};
- // Voronoi diagram builder. Sweepline sweeps from left to right.
+ ///////////////////////////////////////////////////////////////////////////
+ // VORONOI BUILDER ////////////////////////////////////////////////////////
+ ///////////////////////////////////////////////////////////////////////////
+
+ // The sweepline algorithm implementation to compute voronoi diagram of
+ // input data sets of points and segments (Fortune's algorithm).
+ // Complexity - O(N*logN), memory usage - O(N), where N is the total
+ // number of input objects.
+ // Sweepline is a vertical line that sweeps from the left to the right
+ // along the x-axis positive direction. Each of the input objects is
+ // wrapped into the site event. Each event is characterized by its
+ // coordinates: the point site is defined by the point itself,
+ // the segment site is defined by its startpoint. At any moment we
+ // consider only the sites that lie to the left of the sweepline. Beach
+ // line is a curve formed by the parabolic arcs and line segments, that
+ // consists of the points equidistant from the sweepline and the nearest
+ // site to the left of the sweepline. The part of the generated diagram to
+ // the left of the beach line remains unchanged until the end of the
+ // algorithm. Each node in the beach line corresponds to some voronoi edge.
+ // Each node is represented by two sites. Two neighboring nodes has a
+ // common site. Circle event appears at the rightmost point of the circle
+ // tangent to the three sites that correspond to the two consecutive
+ // bisectors. At each step algorithm goes over the leftmost event
+ // and processes it appropriately. This is made until there are no events.
+ // At the end output data structure holds the voronoi diagram of the
+ // initial set of objects.
+ // Each point creates one site event. Each segment creates three site
+ // events: two for its endpoints and one for the segment itself (this is
+ // made to simplify output construction). All the site events are
+ // constructed at the algorithm initialization step. After that they
+ // are ordered using quicksort algorithm.
+ // Priority queue is used to dynamically hold circle events. At each step
+ // of the algorithm the leftmost event is retrieved by comparing the
+ // current site event and the topmost element from the circle event queue.
+ // Standard map container was chosen to hold state of the beach line. The
+ // keys of the map correspond to the bisectors and values to the
+ // corresponding edges from the output data structure. Specially defined
+ // comparison functor is used to make the beach line correctly ordered.
+ // Epsilon-based and high-precision approaches are used to guarantee
+ // efficiency and robustness of the algorithm implementation.
+ // Member data: 1) site_events_ - vector of the site events;
+ // 2) site_event_iterator_ - iterator to the next
+ // site event;
+ // 3) circle_events_ - priority queue of circle events,
+ // allows to retrieve topmost event in O(logN) time;
+ // 4) beach_line_ - contains current state of the beach line;
+ // 5) end_points_ - maps endpoints of the segment sites with
+ // temporary nodes from the beach line. While sweepline
+ // sweeps over the second endpoint of the segments
+ // temporary nodes are being removed;
+ // 6) output_ - contains voronoi diagram itself.
template <typename T>
class voronoi_builder {
public:
typedef T coordinate_type;
- typedef point_2d<coordinate_type> Point2D;
- typedef voronoi_output<coordinate_type> Output;
+ typedef point_2d<coordinate_type> point_2d_type;
+ typedef voronoi_output<coordinate_type> output_type;
typedef site_event<coordinate_type> site_event_type;
typedef circle_event<coordinate_type> circle_event_type;
- voronoi_builder(Output &output) : output_(output) {
- // Set sites iterator.
- site_events_iterator_ = site_events_.begin();
+ voronoi_builder(output_type &output) : output_(output) {
+ // Avoid algorithm fails in case we forgot to init the builder.
+ site_event_iterator_ = site_events_.begin();
}
- // Init beach line before sweepline run.
- // In case of a few first sites situated on the same
- // vertical line, we init beach line with all of them.
- // In other case just use the first two sites for the initialization.
+ // Create site events.
+ // There will be one site event for each input point and three site
+ // events for each input segment (both endpoints of a segment and the
+ // segment itself).
template <typename iType>
void init(const std::vector< point_2d<iType> > &points,
const std::vector< std::pair< point_2d<iType>, point_2d<iType> > > &segments) {
typedef std::pair< point_2d<iType>, point_2d<iType> > iSegment2D;
- // Clear all data structures.
+
+ // Clear output.
output_.clear();
- // TODO(asydorchuk): Add segments intersection check.
+ // TODO(asydorchuk): We may add segment intersection checks there.
// Avoid additional memory reallocations.
site_events_.reserve(points.size() + segments.size() * 3);
- // Create site events from point sites.
+ // Create a site event from each input point.
for (typename std::vector< point_2d<iType> >::const_iterator it = points.begin();
it != points.end(); it++) {
site_events_.push_back(make_site_event(static_cast<T>(it->x()),
- static_cast<T>(it->y()), 0));
+ static_cast<T>(it->y()),
+ 0));
}
-
- // Create site events from end points of segment sites and segment itself.
+
+ // Each segment creates three segment sites:
+ // 1) the startpoint of the segment;
+ // 2) the endpoint of the segment;
+ // 3) the segment itself.
for (typename std::vector<iSegment2D>::const_iterator it = segments.begin();
it != segments.end(); it++) {
T x1 = static_cast<T>(it->first.x());
@@ -1660,68 +2018,73 @@
site_events_.push_back(make_site_event(x1, y1, x2, y2, 0));
}
- // Sort site events.
+ // Sort the site events.
sort(site_events_.begin(), site_events_.end());
// Remove duplicates.
- site_events_.erase(unique(site_events_.begin(), site_events_.end()),
- site_events_.end());
+ site_events_.erase(unique(
+ site_events_.begin(), site_events_.end()), site_events_.end());
- // Number sites.
- for (int cur_index = 0; cur_index < static_cast<int>(site_events_.size()); cur_index++)
- site_events_[cur_index].set_site_index(cur_index);
-
- // Set site iterator.
- site_events_iterator_ = site_events_.begin();
-
- // Init output with number of site events.
- // There will be one site event for each input point and three site events
- // for each input segment: both ends of a segment and the segment itself.
+ // Number the sites.
+ for (size_t cur = 0; cur < site_events_.size(); cur++)
+ site_events_[cur].index(cur);
+
+ // Init the site's iterator.
+ site_event_iterator_ = site_events_.begin();
+
+ // Init the output data structure.
output_.init(site_events_.size());
}
void clear() {
- circle_events_.reset();
beach_line_.clear();
+ circle_events_.clear();
site_events_.clear();
}
+ // Run the sweepline algorithm.
void run_sweepline() {
- // Init beach line.
+ // Init the beach line.
if (site_events_.empty()) {
+ // No input sites.
return;
} else if (site_events_.size() == 1) {
// Handle one input site case.
output_.process_single_site(site_events_[0]);
- site_events_iterator_++;
+ site_event_iterator_++;
} else {
int skip = 0;
- // Init beach line.
- while(site_events_iterator_ != site_events_.end() &&
- site_events_iterator_->x() == site_events_.begin()->x() &&
- site_events_iterator_->is_vertical()) {
- site_events_iterator_++;
+
+ // Count the number of the sites to init the beach line.
+ while(site_event_iterator_ != site_events_.end() &&
+ site_event_iterator_->x() == site_events_.begin()->x() &&
+ site_event_iterator_->is_vertical()) {
+ site_event_iterator_++;
skip++;
}
if (skip == 1) {
- // Init beach line with two sites.
+ // Init the beach line with the two first sites.
init_beach_line();
} else {
- // Init beach line with sites situated on the same vertical line.
+ // Init the beach line with the sites situated on the same
+ // vertical line. This could be a set of point and vertical
+ // segment sites.
init_beach_line_collinear_sites();
}
}
- // Algorithm stops when there are no events to process.
- while (!circle_events_.empty() ||
- !(site_events_iterator_ == site_events_.end())) {
+ // The algorithm stops when there are no events to process.
+ // The circle events with the same rightmost point as the next
+ // site event go first.
+ while (!circle_events_.empty() ||
+ !(site_event_iterator_ == site_events_.end())) {
if (circle_events_.empty()) {
process_site_event();
- } else if (site_events_iterator_ == site_events_.end()) {
+ } else if (site_event_iterator_ == site_events_.end()) {
process_circle_event();
} else {
- if (circle_events_.top().compare(*site_events_iterator_) == MORE) {
+ if (circle_events_.top().compare(*site_event_iterator_) == MORE) {
process_site_event();
} else {
process_circle_event();
@@ -1729,220 +2092,270 @@
}
}
- // Simplify output.
+ // Simplify the output (remove zero-length edges).
output_.simplify();
clear();
}
private:
- typedef typename std::vector<site_event_type>::const_iterator site_events_iterator_type;
- typedef typename Output::voronoi_edge_type edge_type;
- typedef circle_events_queue<coordinate_type> CircleEventsQueue;
- typedef beach_line_node<coordinate_type> Key;
- typedef beach_line_node_data<coordinate_type> Value;
- typedef node_comparer<Key> NodeComparer;
- typedef std::map< Key, Value, NodeComparer > BeachLine;
- typedef typename std::map< Key, Value, NodeComparer >::iterator beach_line_iterator;
- typedef typename std::pair<Point2D, beach_line_iterator> end_point_type;
+ typedef typename std::vector<site_event_type>::const_iterator
+ site_event_iterator_type;
+ typedef typename output_type::voronoi_edge_type edge_type;
+ typedef circle_events_queue<coordinate_type> circle_event_queue_type;
+ typedef beach_line_node<coordinate_type> key_type;
+ typedef beach_line_node_data<coordinate_type> value_type;
+ typedef node_comparer<key_type> node_comparer_type;
+ typedef std::map< key_type, value_type, node_comparer_type >
+ beach_line_type;
+ typedef typename beach_line_type::iterator beach_line_iterator;
+ typedef std::pair<point_2d_type, beach_line_iterator> end_point_type;
+ // Init the beach line with the two first sites.
+ // The first site is always a point.
void init_beach_line() {
- // The first site is always a point.
// Get the first and the second site events.
- site_events_iterator_type it_first = site_events_.begin();
- site_events_iterator_type it_second = site_events_.begin();
+ site_event_iterator_type it_first = site_events_.begin();
+ site_event_iterator_type it_second = site_events_.begin();
it_second++;
- // Insert new nodes.
+ // Update the beach line.
insert_new_arc(*it_first, *it_first, *it_second, beach_line_.begin());
// The second site has been already processed.
- site_events_iterator_++;
+ // Move the site's iterator.
+ site_event_iterator_++;
}
// Insert bisectors for all collinear initial sites.
- // There should be at least two colinear sites.
void init_beach_line_collinear_sites() {
- site_events_iterator_type it_first = site_events_.begin();
- site_events_iterator_type it_second = site_events_.begin();
+ site_event_iterator_type it_first = site_events_.begin();
+ site_event_iterator_type it_second = site_events_.begin();
it_second++;
- while (it_second != site_events_iterator_) {
- // Create new beach line node.
- Key new_node(*it_first, *it_second);
-
- // Update output.
+ while (it_second != site_event_iterator_) {
+ // Create a new beach line node.
+ key_type new_node(*it_first, *it_second);
+
+ // Update the output.
edge_type *edge = output_.insert_new_edge(*it_first, *it_second);
- // Insert new node into the binary search tree.
- beach_line_.insert(std::pair<Key, Value>(new_node, Value(edge)));
-
+ // Insert a new bisector into the beach line.
+ beach_line_.insert(
+ std::pair<key_type, value_type>(new_node, value_type(edge)));
+
// Update iterators.
it_first++;
it_second++;
}
}
- // Uses special comparison function for the lower bound and insertion
- // operations.
+ // Process a single site.
void process_site_event() {
- site_event_type site_event = *site_events_iterator_;
- site_events_iterator_++;
+ // Get the site event to process.
+ site_event_type site_event = *site_event_iterator_;
- // If new site is end point of some segment, remove temporary nodes from
- // beach line data structure.
+ // Move the site's iterator.
+ site_event_iterator_++;
+
+ // If a new site is an end point of some segment,
+ // remove temporary nodes from the beach line data structure.
if (!site_event.is_segment()) {
- while (!end_points_.empty() && end_points_.top().first == site_event.get_point0()) {
+ while (!end_points_.empty() &&
+ end_points_.top().first == site_event.point0()) {
beach_line_.erase(end_points_.top().second);
end_points_.pop();
}
}
+ // Create degenerate node.
+ key_type new_key(site_event);
+
// Find the node in the binary search tree with left arc
// lying above the new site point.
- Key new_key(site_event);
beach_line_iterator it = beach_line_.lower_bound(new_key);
int it_dist = site_event.is_segment() ? 2 : 1;
+ // Do further processing depending on the above node position.
+ // For any two neighbouring nodes the second site of the first node
+ // is the same as the first site of the second arc.
if (it == beach_line_.end()) {
+ // The above arc corresponds to the second arc of the last node.
+ // Move the iterator to the last node.
it--;
- const site_event_type &site_arc = it->first.get_right_site();
- // Insert new arc into the sweepline.
- beach_line_iterator new_node_it = insert_new_arc(site_arc, site_arc, site_event, it);
- std::advance(new_node_it, -it_dist);
+ // Get the second site of the last node
+ const site_event_type &site_arc = it->first.right_site();
+
+ // Insert new nodes into the beach line. Update the output.
+ beach_line_iterator new_node_it =
+ insert_new_arc(site_arc, site_arc, site_event, it);
- // Add candidate circle to the event queue.
- activate_circle_event(it->first.get_left_site(), it->first.get_right_site(),
+ // Add a candidate circle to the circle event queue.
+ // There could be only one new circle event formed by
+ // a new bisector and the one on the left.
+ std::advance(new_node_it, -it_dist);
+ activate_circle_event(it->first.left_site(),
+ it->first.right_site(),
site_event, new_node_it);
} else if (it == beach_line_.begin()) {
- const site_event_type &site_arc = it->first.get_left_site();
+ // The above arc corresponds to the first site of the first node.
+ const site_event_type &site_arc = it->first.left_site();
- // Insert new arc into the sweepline.
+ // Insert new nodes into the beach line. Update the output.
insert_new_arc(site_arc, site_arc, site_event, it);
- // Add candidate circle to the event queue.
+ // If the site event is a segment, update its direction.
if (site_event.is_segment()) {
- site_event.set_inverse();
+ site_event.inverse();
}
- activate_circle_event(site_event, it->first.get_left_site(),
- it->first.get_right_site(), it);
+
+ // Add a candidate circle to the circle event queue.
+ // There could be only one new circle event formed by
+ // a new bisector and the one on the right.
+ activate_circle_event(site_event, it->first.left_site(),
+ it->first.right_site(), it);
} else {
- const site_event_type &site_arc2 = it->first.get_left_site();
- const site_event_type &site3 = it->first.get_right_site();
+ // The above arc corresponds neither to the first,
+ // nor to the last site in the beach line.
+ const site_event_type &site_arc2 = it->first.left_site();
+ const site_event_type &site3 = it->first.right_site();
- // Remove candidate circle from the event queue.
+ // Remove the candidate circle from the event queue.
it->second.deactivate_circle_event();
it--;
- const site_event_type &site_arc1 = it->first.get_right_site();
- const site_event_type &site1 = it->first.get_left_site();
+ const site_event_type &site_arc1 = it->first.right_site();
+ const site_event_type &site1 = it->first.left_site();
- // Insert new arc into the sweepline.
+ // Insert new nodes into the beach line. Update the output.
beach_line_iterator new_node_it =
insert_new_arc(site_arc1, site_arc2, site_event, it);
- // Add candidate circles to the event queue.
+ // Add candidate circles to the circle event queue.
+ // There could be up to two circle events formed by
+ // a new bisector and the one on the left or right.
std::advance(new_node_it, -it_dist);
- activate_circle_event(site1, site_arc1, site_event, new_node_it);
+ activate_circle_event(site1, site_arc1, site_event,
+ new_node_it);
+
+ // If the site event is a segment, update its direction.
if (site_event.is_segment()) {
- site_event.set_inverse();
+ site_event.inverse();
}
std::advance(new_node_it, it_dist + 1);
- activate_circle_event(site_event, site_arc2, site3, new_node_it);
+ activate_circle_event(site_event, site_arc2, site3,
+ new_node_it);
}
}
- // Doesn't use special comparison function as it works fine only for
- // the site events processing.
+ // Process a single circle event.
+ // In general case circle event is made of the three consequtive sites
+ // that form two bisector nodes in the beach line data structure.
+ // Let circle event sites be A, B, C, two bisectors that define
+ // circle event be (A, B), (B, C). During circle event processing
+ // we remove (A, B), (B, C) and insert (A, C). As beach line comparer
+ // works correctly only if one of the nodes is a new one we remove
+ // (B, C) bisector and change (A, B) bisector to the (A, C). That's
+ // why we use const_cast there and take all the responsibility that
+ // map data structure keeps correct ordering.
void process_circle_event() {
+ // Get the topmost circle event.
const circle_event_type &circle_event = circle_events_.top();
-
- // Retrieve the second bisector node that corresponds to the given
- // circle event.
- beach_line_iterator it_first = circle_event.get_bisector();
+ beach_line_iterator it_first = circle_event.bisector();
beach_line_iterator it_last = it_first;
- // Get the the third site.
- site_event_type site3 = it_first->first.get_right_site();
+ // Get the C site.
+ site_event_type site3 = it_first->first.right_site();
+
+ // Get the half-edge corresponding to the second bisector - (B, C).
+ edge_type *bisector2 = it_first->second.edge();
- // Get second bisector;
- edge_type *bisector2 = it_first->second.get_edge();
-
- // Get first bisector;
+ // Get the half-edge corresponding to the first bisector - (A, B).
it_first--;
- edge_type *bisector1 = it_first->second.get_edge();
-
- // Get the first site that creates given circle event.
- site_event_type site1 = it_first->first.get_left_site();
-
- // Let circle event sites be A, B, C, two bisectors that define
- // circle event be (A, B), (B, C). During circle event processing
- // we remove (A, B), (B, C) and insert (A, C). As beach line nodes
- // comparer doesn't work fine for the circle events we only remove
- // (B, C) bisector and change (A, B) bisector to the (A, C). That's
- // why we use const_cast there and take all the responsibility that
- // map data structure keeps correct ordering.
+ edge_type *bisector1 = it_first->second.edge();
+
+ // Get the A site.
+ site_event_type site1 = it_first->first.left_site();
+
if (!site1.is_segment() && site3.is_segment() &&
- site3.get_point1(true) == site1.get_point0()) {
- site3.set_inverse();
+ site3.point1(true) == site1.point0()) {
+ site3.inverse();
}
- const_cast<Key &>(it_first->first).set_right_site(site3);
- it_first->second.set_edge(output_.insert_new_edge(site1, site3, circle_event,
- bisector1, bisector2));
+
+ // Change the (A, B) bisector node to the (A, C) bisector node.
+ const_cast<key_type &>(it_first->first).right_site(site3);
+
+ // Insert the new bisector into the beach line.
+ it_first->second.edge(output_.insert_new_edge(site1, site3, circle_event,
+ bisector1, bisector2));
+
+ // Remove the (B, C) bisector node from the beach line.
beach_line_.erase(it_last);
it_last = it_first;
- // Pop circle event from the event queue, before we might
- // insert new events in it.
+ // Pop the topmost circle event from the event queue.
circle_events_.pop();
- // Check the new triplets formed by the neighboring arcs
- // for potential circle events. Check left.
+ // Check new triplets formed by the neighboring arcs
+ // to the left for potential circle events.
if (it_first != beach_line_.begin()) {
it_first->second.deactivate_circle_event();
it_first--;
- const site_event_type &site_l1 = it_first->first.get_left_site();
+ const site_event_type &site_l1 = it_first->first.left_site();
activate_circle_event(site_l1, site1, site3, it_last);
}
- // Check the new triplets formed by the neighboring arcs
- // for potential circle events. Check right.
+ // Check the new triplet formed by the neighboring arcs
+ // to the right for potential circle events.
it_last++;
if (it_last != beach_line_.end()) {
it_last->second.deactivate_circle_event();
- const site_event_type &site_r1 = it_last->first.get_right_site();
+ const site_event_type &site_r1 = it_last->first.right_site();
activate_circle_event(site1, site3, site_r1, it_last);
}
-
}
- // Insert new arc below site arc into the beach line.
+ // Insert new nodes into the beach line. Update the output.
beach_line_iterator insert_new_arc(const site_event_type &site_arc1,
const site_event_type &site_arc2,
const site_event_type &site_event,
const beach_line_iterator &position) {
- // Create two new nodes.
- Key new_left_node(site_arc1, site_event);
- Key new_right_node(site_event, site_arc2);
+ // Create two new bisectors with opposite directions.
+ key_type new_left_node(site_arc1, site_event);
+ key_type new_right_node(site_event, site_arc2);
+
+ // Set correct orientation for the first site of the second node.
if (site_event.is_segment()) {
- new_right_node.set_left_site_inverse();
+ new_right_node.inverse_left_site();
}
-
- // Insert two new nodes into the binary search tree.
- // Update output.
+
+ // Update the output.
edge_type *edge = output_.insert_new_edge(site_arc2, site_event);
+
+ // Update the beach line with the (site_arc1, site_event) bisector.
beach_line_iterator it = beach_line_.insert(position,
- std::pair<Key, Value>(new_right_node, Value(edge->twin())));
+ typename beach_line_type::value_type(new_right_node, value_type(edge->twin())));
+
if (site_event.is_segment()) {
- Key new_node(site_event, site_event);
- new_node.set_right_site_inverse();
+ // Update the beach line with temporary bisector, that will
+ // disappear after processing site event going through the
+ // endpoint of the segment site.
+ key_type new_node(site_event, site_event);
+ new_node.inverse_right_site();
beach_line_iterator temp_it = beach_line_.insert(position,
- std::pair<Key, Value>(new_node, Value(NULL)));
- end_points_.push(std::make_pair(site_event.get_point1(), temp_it));
+ typename beach_line_type::value_type(new_node, value_type(NULL)));
+
+ // Update the data structure that holds temporary bisectors.
+ end_points_.push(std::make_pair(site_event.point1(), temp_it));
}
- beach_line_.insert(position, std::pair<Key, Value>(new_left_node, Value(edge)));
+
+ // Update the beach line with the (site_event, site_arc2) bisector.
+ beach_line_.insert(position,
+ typename beach_line_type::value_type(new_left_node, value_type(edge)));
return it;
}
- // Create circle event from the given three points.
+ // Create a circle event from the given three sites.
+ // Returns true if the circle event exists, else false.
+ // If exists circle event is saved into the c_event variable.
bool create_circle_event(const site_event_type &site1,
const site_event_type &site2,
const site_event_type &site3,
@@ -1950,43 +2363,60 @@
if (!site1.is_segment()) {
if (!site2.is_segment()) {
if (!site3.is_segment()) {
+ // (point, point, point) sites.
return create_circle_event_ppp(site1, site2, site3, c_event);
} else {
+ // (point, point, segment) sites.
return create_circle_event_pps(site1, site2, site3, 3, c_event);
}
} else {
if (!site3.is_segment()) {
+ // (point, segment, point) sites.
return create_circle_event_pps(site1, site3, site2, 2, c_event);
} else {
+ // (point, segment, segment) sites.
return create_circle_event_pss(site1, site2, site3, 1, c_event);
}
}
} else {
if (!site2.is_segment()) {
if (!site3.is_segment()) {
+ // (segment, point, point) sites.
return create_circle_event_pps(site2, site3, site1, 1, c_event);
} else {
+ // (segment, point, segment) sites.
return create_circle_event_pss(site2, site1, site3, 2, c_event);
}
} else {
if (!site3.is_segment()) {
+ // (segment, segment, point) sites.
return create_circle_event_pss(site3, site1, site2, 3, c_event);
} else {
+ // (segment, segment, segment) sites.
return create_circle_event_sss(site1, site2, site3, c_event);
}
}
}
}
- // Add new circle event to the event queue.
+ // Add a new circle event to the event queue.
+ // bisector_node corresponds to the (site2, site3) bisector.
void activate_circle_event(const site_event_type &site1,
const site_event_type &site2,
const site_event_type &site3,
beach_line_iterator bisector_node) {
circle_event_type c_event;
+ // Check if the three input sites create a circle event.
if (create_circle_event(site1, site2, site3, c_event)) {
- c_event.set_bisector(bisector_node);
- bisector_node->second.activate_circle_event(circle_events_.push(c_event));
+ // Update circle event's bisector iterator to point to the
+ // second bisector that forms it in the beach line.
+ c_event.bisector(bisector_node);
+
+ // Add the new circle event to the circle events queue.
+ // Update bisector's circle event iterator to point to the
+ // new circle event in the circle event queue.
+ bisector_node->second.activate_circle_event(
+ circle_events_.push(c_event));
}
}
@@ -1998,12 +2428,12 @@
};
std::vector<site_event_type> site_events_;
- site_events_iterator_type site_events_iterator_;
+ site_event_iterator_type site_event_iterator_;
std::priority_queue< end_point_type, std::vector<end_point_type>,
end_point_comparison > end_points_;
- CircleEventsQueue circle_events_;
- BeachLine beach_line_;
- Output &output_;
+ circle_event_queue_type circle_events_;
+ beach_line_type beach_line_;
+ output_type &output_;
//Disallow copy constructor and operator=
voronoi_builder(const voronoi_builder&);
@@ -2014,7 +2444,4 @@
} // boost
} // detail
-#undef INT_PREDICATE_CONVERT_65_BIT
-#undef INT_PREDICATE_COMPUTE_DIFFERENCE
-
#endif
Modified: sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_output.hpp
==============================================================================
--- sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_output.hpp (original)
+++ sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_output.hpp 2010-12-20 09:30:57 EST (Mon, 20 Dec 2010)
@@ -1,11 +1,11 @@
-// Boost sweepline/voronoi_output.hpp header file
+// Boost sweepline/voronoi_output.hpp header file
// Copyright Andrii Sydorchuk 2010.
// 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.
+// See http://www.boost.org for updates, documentation, and revision history.
#ifndef BOOST_SWEEPLINE_VORONOI_OUTPUT
#define BOOST_SWEEPLINE_VORONOI_OUTPUT
@@ -26,7 +26,7 @@
class circle_event;
template <typename T>
- class epsilon_robust_comparator;
+ class robust_voronoi_vertex;
template <typename T>
class voronoi_builder;
@@ -41,7 +41,7 @@
template <typename T>
class voronoi_output;
- // Represents 2D point.
+ // Cartesian 2D point data structure.
template <typename T>
class point_2d {
public:
@@ -49,37 +49,37 @@
point_2d() {}
- point_2d(T x, T y) {
+ point_2d(coordinate_type x, coordinate_type y) {
x_ = x;
y_ = y;
}
- bool operator==(const point_2d &point) const {
- return (this->x_ == point.x()) && (this->y_ == point.y());
+ bool operator==(const point_2d &that) const {
+ return (this->x_ == that.x()) && (this->y_ == that.y());
}
- bool operator!=(const point_2d &point) const {
- return (this->x_ != point.x()) || (this->y_ != point.y());
+ bool operator!=(const point_2d &that) const {
+ return (this->x_ != that.x()) || (this->y_ != that.y());
}
- // This comparator actually defines sweepline movement direction.
- bool operator<(const point_2d &point) const {
- // Sweepline sweeps from left to right.
- if (this->x_ != point.x_)
- return this->x_ < point.x_;
- return this->y_ < point.y_;
+ // Comparison function.
+ // Defines ordering of the points on the 2D plane.
+ bool operator<(const point_2d &that) const {
+ if (this->x_ != that.x_)
+ return this->x_ < that.x_;
+ return this->y_ < that.y_;
}
- bool operator<=(const point_2d &point) const {
- return !(point < (*this));
+ bool operator<=(const point_2d &that) const {
+ return !(that < (*this));
}
- bool operator>(const point_2d &point) const {
- return point < (*this);
+ bool operator>(const point_2d &that) const {
+ return that < (*this);
}
- bool operator>=(const point_2d &point) const {
- return !((*this) < point);
+ bool operator>=(const point_2d &that) const {
+ return !((*this) < that);
}
coordinate_type x() const {
@@ -104,19 +104,21 @@
};
template <typename T>
- point_2d<T> make_point_2d(T x, T y) {
+ static inline point_2d<T> make_point_2d(T x, T y) {
return point_2d<T>(x,y);
}
- // Bounding rectangle data structure.
+ // Bounding rectangle data structure. Contains coordinates
+ // of the bottom left and the upper right points of the rectangle.
template <typename T>
class BRect {
public:
typedef T coordinate_type;
+ typedef point_2d<coordinate_type> point_2d_type;
BRect() {}
- BRect(const point_2d<T> &p1, const point_2d<T> &p2) {
+ BRect(const point_2d_type &p1, const point_2d_type &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());
@@ -131,33 +133,41 @@
y_max_ = (std::max)(y_mn, y_mx);
}
- void update(const point_2d<T> &p) {
+ // Extend the rectangle with a new point.
+ void update(const point_2d_type &p) {
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(const point_2d<T> &p) const {
- return p.x() > x_min_ && p.x() < x_max_ && p.y() > y_min_ && p.y() < y_max_;
+ // Check whether a point is situated inside the bounding rectangle.
+ bool contains(const point_2d_type &p) const {
+ return p.x() > x_min_ && p.x() < x_max_ &&
+ p.y() > y_min_ && p.y() < y_max_;
}
+ // Check whether the bounding rectangle has a non-zero area.
bool is_valid() const {
return (x_min_ < x_max_) && (y_min_ < y_max_);
}
+ // Return the x-coordinate of the bottom left point of the rectangle.
coordinate_type x_min() const {
return x_min_;
}
+ // Return the y-coordinate of the bottom left point of the rectangle.
coordinate_type y_min() const {
return y_min_;
}
+ // Return the x-coordinate of the upper right point of the rectangle.
coordinate_type x_max() const {
return x_max_;
}
+ // Return the y-coordinate of the upper right point of the rectangle.
coordinate_type y_max() const {
return y_max_;
}
@@ -169,366 +179,486 @@
coordinate_type y_max_;
};
- // Contains output visualization tools.
+ // Voronoi output postprocessing tools.
template <typename T>
class voronoi_helper {
public:
typedef T coordinate_type;
- typedef point_2d<T> Point2D;
+ typedef point_2d<coordinate_type> point_2d_type;
+ typedef BRect<coordinate_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,
};
- // Find edge(segment, ray, line) rectangle intersetion points.
- static bool find_intersections(const Point2D &origin, const Point2D &direction,
- kEdgeType edge_type, const BRect<coordinate_type> &brect,
- std::vector< Point2D > &intersections) {
- coordinate_type half = static_cast<coordinate_type>(0.5);
-
- // Find rectangle center.
- coordinate_type center_x = (brect.x_min() + brect.x_max()) * half;
- coordinate_type center_y = (brect.y_min() + brect.y_max()) * half;
+ // Get a view rectangle based on the voronoi bounding rectangle.
+ static BRect<coordinate_type> get_view_rectangle(
+ const BRect<coordinate_type> &brect) {
+ coordinate_type x_len = (brect.x_max() - brect.x_min());
+ coordinate_type y_len = (brect.y_max() - brect.y_min());
+ coordinate_type x_mid = (brect.x_max() + brect.x_min()) * 0.5;
+ coordinate_type y_mid = (brect.y_max() + brect.y_min()) * 0.5;
+ coordinate_type offset = x_len;
+ if (offset < y_len)
+ offset = y_len;
+ if (offset == 0.0)
+ offset = 0.5;
+ BRect<coordinate_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.
+ static std::vector<point_2d_type> get_intermediate_points(
+ const voronoi_edge<coordinate_type> *edge,
+ const BRect<coordinate_type> &brect) {
+ // Retrieve the cell to the left of the current edge.
+ const voronoi_cell<coordinate_type> *cell1 = edge->cell();
+
+ // Retrieve the cell to the right of the current edge.
+ const voronoi_cell<coordinate_type> *cell2 = edge->twin()->cell();
+
+ // edge_points - contains intermediate points.
+ std::vector<point_2d_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(edge->vertex0()->vertex());
+ edge_points.push_back(edge->vertex1()->vertex());
+ } else {
+ // Edge is a ray or line.
+ // Clip it with the bounding rectangle.
+ const point_2d_type &site1 = cell1->point0();
+ const point_2d_type &site2 = cell2->point0();
+ point_2d_type origin((site1.x() + site2.x()) * 0.5,
+ (site1.y() + site2.y()) * 0.5);
+ point_2d_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] = edge->vertex1()->vertex();
+ if (edge->vertex0() != NULL)
+ edge_points[0] = edge->vertex0()->vertex();
+ }
+ } else {
+ // point1 - site point;
+ const point_2d_type &point1 = (cell1->contains_segment()) ?
+ cell2->point0() : cell1->point0();
+
+ // point2 - startpoint of the segment site;
+ const point_2d_type &point2 = (cell1->contains_segment()) ?
+ cell1->point0() : cell2->point0();
+
+ // point3 - endpoint of the segment site;
+ const point_2d_type &point3 = (cell1->contains_segment()) ?
+ cell1->point1() : cell2->point1();
+
+ if (edge->vertex0() != NULL && edge->vertex1() != NULL) {
+ // Edge is a segment or parabolic arc.
+ edge_points.push_back(edge->vertex0()->vertex());
+ edge_points.push_back(edge->vertex1()->vertex());
+ fill_intermediate_points(point1, point2, point3,
+ edge_points);
+ } else {
+ // Edge is a ray or line.
+ // Clip it with the bounding rectangle.
+ coordinate_type dir_x =
+ (cell1->contains_segment() ^ (point1 == point3)) ?
+ point3.y() - point2.y() : point2.y() - point3.y();
+ coordinate_type dir_y =
+ (cell1->contains_segment() ^ (point1 == point3)) ?
+ point2.x() - point3.x() : point3.x() - point2.x();
+ point_2d_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] = edge->vertex1()->vertex();
+ if (edge->vertex0() != NULL)
+ edge_points[0] = edge->vertex0()->vertex();
+ }
+ }
+ return edge_points;
+ }
+
+ // Find edge-rectangle intersection points.
+ // Edge is represented by its origin, direction and type.
+ static void find_intersections(
+ const point_2d_type &origin, const point_2d_type &direction,
+ kEdgeType edge_type, const brect_type &brect,
+ std::vector<point_2d_type> &intersections) {
+ // Find the center of the rectangle.
+ coordinate_type center_x = (brect.x_min() + brect.x_max()) * 0.5;
+ coordinate_type center_y = (brect.y_min() + brect.y_max()) * 0.5;
- // Find rectangle half-diagonal vector.
+ // Find the half-diagonal vector of the rectangle.
coordinate_type len_x = brect.x_max() - center_x;
coordinate_type len_y = brect.y_max() - center_y;
-
- // Find distance between origin and center of rectangle.
+
+ // Find the vector between the origin and the center of the
+ // rectangle.
coordinate_type diff_x = origin.x() - center_x;
coordinate_type diff_y = origin.y() - center_y;
-
- // Find direction perpendicular to the current.
+
+ // Find the vector perpendicular to the direction vector.
coordinate_type perp_x = direction.y();
coordinate_type perp_y = -direction.x();
- // Compare projections of distances.
- coordinate_type lexpr = magnitude(perp_x * diff_x + perp_y * diff_y);
- coordinate_type rexpr = magnitude(perp_x * len_x) + magnitude(perp_y * len_y);
+ // Projection of the vector between the origin and the center of
+ // the rectangle on the line perpendicular to the direction vector.
+ coordinate_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.
+ coordinate_type rexpr = magnitude(perp_x * len_x) +
+ magnitude(perp_y * len_y);
- // Intersection check.
+ // Intersection check. Compare projections.
if (lexpr > rexpr)
- return false;
+ return;
- // Intersection parameters:
- // origin + fT[i] * direction = intersections point, i = 0 .. 1.
+ // 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;
coordinate_type fT0 = 0;
coordinate_type fT1 = 0;
- // Find intersections with lines going through sides of a 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);
+ // 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(make_point_2d(origin.x() + fT0 * direction.x(),
- origin.y() + fT0 * direction.y()));
+ intersections.push_back(make_point_2d(
+ origin.x() + fT0 * direction.x(),
+ origin.y() + fT0 * direction.y()));
if (fT1_used && fT0 != fT1 && check_extent(fT1, edge_type))
- intersections.push_back(make_point_2d(origin.x() + fT1 * direction.x(),
- origin.y() + fT1 * direction.y()));
- return fT0_used && fT1_used;
+ intersections.push_back(make_point_2d(
+ origin.x() + fT1 * direction.x(),
+ origin.y() + fT1 * direction.y()));
}
- static void fill_intermediate_points(Point2D point_site,
- Point2D segment_site_start,
- Point2D segment_site_end,
- std::vector<Point2D> &intermediate_points) {
+ private:
+ voronoi_helper();
+
+ // Find intermediate points of the parabola.
+ // 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.
+ static void fill_intermediate_points(
+ point_2d_type point_site,
+ point_2d_type segment_site_start,
+ point_2d_type segment_site_end,
+ std::vector<point_2d_type> &intermediate_points) {
+ // Get the number of intermediate points to represent parabola.
int num_inter_points = get_intermediate_points_count(
intermediate_points[0], intermediate_points[1]);
+
+ // Check if there are any intermediate points.
+ // In case bisector is a segment there won't be any.
if (num_inter_points <= 0 ||
point_site == segment_site_start ||
point_site == segment_site_end) {
return;
}
+
+ // Reserve space to avoid vector reallocations.
intermediate_points.reserve(2 + num_inter_points);
- coordinate_type segm_vec_x = segment_site_end.x() - segment_site_start.x();
- coordinate_type segm_vec_y = segment_site_end.y() - segment_site_start.y();
- coordinate_type sqr_segment_length = segm_vec_x * segm_vec_x + segm_vec_y * segm_vec_y;
- coordinate_type projection_start =
- get_point_projection(intermediate_points[0], segment_site_start, segment_site_end);
- coordinate_type projection_end =
- get_point_projection(intermediate_points[1], segment_site_start, segment_site_end);
+
+ // 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.
+ coordinate_type segm_vec_x = segment_site_end.x() -
+ segment_site_start.x();
+ coordinate_type segm_vec_y = segment_site_end.y() -
+ segment_site_start.y();
+ coordinate_type sqr_segment_length = segm_vec_x * segm_vec_x +
+ segm_vec_y * segm_vec_y;
+
+ // Compute coordinates of the endpoints of the edge
+ // in the transformed space.
+ coordinate_type projection_start = get_point_projection(
+ intermediate_points[0], segment_site_start, segment_site_end);
+ coordinate_type projection_end = get_point_projection(
+ intermediate_points[1], segment_site_start, segment_site_end);
coordinate_type step = (projection_end - projection_start) *
- sqr_segment_length / (num_inter_points + 1);
- coordinate_type point_vec_x = point_site.x() - segment_site_start.x();
- coordinate_type point_vec_y = point_site.y() - segment_site_start.y();
- coordinate_type point_rot_x = segm_vec_x * point_vec_x + segm_vec_y * point_vec_y;
- coordinate_type point_rot_y = segm_vec_x * point_vec_y - segm_vec_y * point_vec_x;
- coordinate_type projection_cur_x = projection_start * sqr_segment_length + step;
- point_2d<T> last_point = intermediate_points.back();
+ sqr_segment_length / (num_inter_points + 1);
+
+ // Compute parabola parameterers in the transformed space.
+ // Parabola has next representation:
+ // f(x) = ((x-point_rot_x)^2 + point_rot_y^2) / (2.0*point_rot_y).
+ coordinate_type point_vec_x = point_site.x() -
+ segment_site_start.x();
+ coordinate_type point_vec_y = point_site.y() -
+ segment_site_start.y();
+ coordinate_type point_rot_x = segm_vec_x * point_vec_x +
+ segm_vec_y * point_vec_y;
+ coordinate_type point_rot_y = segm_vec_x * point_vec_y -
+ segm_vec_y * point_vec_x;
+
+ // Compute the x-coordinate of the first intermediate point.
+ coordinate_type projection_cur_x = projection_start *
+ sqr_segment_length + step;
+
+ // Temporary remove the last point.
+ point_2d_type last_point = intermediate_points.back();
intermediate_points.pop_back();
- for (int i = 0; i < num_inter_points; i++, projection_cur_x += step) {
+
+ // Generate intermediate points.
+ for (int i = 0; i < num_inter_points;
+ i++, projection_cur_x += step) {
+ // Compute coordinates of the intermediate points
+ // in the transformed space.
coordinate_type inter_rot_x = projection_cur_x;
coordinate_type inter_rot_y =
- ((inter_rot_x - point_rot_x) * (inter_rot_x - point_rot_x) +
+ ((inter_rot_x - point_rot_x) *
+ (inter_rot_x - point_rot_x) +
point_rot_y * point_rot_y) / (2.0 * point_rot_y);
+
+ // Compute coordinates of the intermediates points
+ // in the initial space.
coordinate_type new_point_x =
(segm_vec_x * inter_rot_x - segm_vec_y * inter_rot_y) /
sqr_segment_length + segment_site_start.x();
coordinate_type new_point_y =
(segm_vec_x * inter_rot_y + segm_vec_y * inter_rot_x) /
sqr_segment_length + segment_site_start.y();
- intermediate_points.push_back(make_point_2d(new_point_x, new_point_y));
- }
- intermediate_points.push_back(last_point);
- }
- static BRect<coordinate_type> get_view_rectangle(const BRect<coordinate_type> &brect) {
- coordinate_type x_len = (brect.x_max() - brect.x_min());
- coordinate_type y_len = (brect.y_max() - brect.y_min());
- coordinate_type x_mid = (brect.x_max() + brect.x_min()) / 2;
- coordinate_type y_mid = (brect.y_max() + brect.y_min()) / 2;
- coordinate_type offset = x_len;
- if (offset < y_len)
- offset = y_len;
- if (offset == 0.0)
- offset = 0.5;
- BRect<coordinate_type> new_brect(x_mid - offset, y_mid - offset,
- x_mid + offset, y_mid + offset);
- return new_brect;
- }
-
- static std::vector< point_2d<coordinate_type> > get_intermediate_points(
- const voronoi_edge<coordinate_type> *edge, const BRect<T> &brect) {
- const voronoi_cell<coordinate_type> *cell1 = edge->cell();
- const voronoi_cell<coordinate_type> *cell2 = edge->twin()->cell();
- std::vector<Point2D> edge_points;
- if (!(cell1->contains_segment() ^ cell2->contains_segment())) {
- if (edge->vertex0() != NULL && edge->vertex1() != NULL) {
- edge_points.push_back(edge->vertex0()->vertex());
- edge_points.push_back(edge->vertex1()->vertex());
- } else {
- const Point2D &site1 = cell1->get_point0();
- const Point2D &site2 = cell2->get_point0();
- Point2D origin((site1.x() + site2.x()) * static_cast<coordinate_type>(0.5),
- (site1.y() + site2.y()) * static_cast<coordinate_type>(0.5));
- Point2D direction(site1.y() - site2.y(), site2.x() - site1.x());
- find_intersections(origin, direction, LINE, brect, edge_points);
- if (edge->vertex1() != NULL)
- edge_points[1] = edge->vertex1()->vertex();
- if (edge->vertex0() != NULL)
- edge_points[0] = edge->vertex0()->vertex();
- }
- } else {
- const Point2D &point1 = (cell1->contains_segment()) ?
- cell2->get_point0() : cell1->get_point0();
- const Point2D &point2 = (cell1->contains_segment()) ?
- cell1->get_point0() : cell2->get_point0();
- const Point2D &point3 = (cell1->contains_segment()) ?
- cell1->get_point1() : cell2->get_point1();
- if (edge->vertex0() != NULL && edge->vertex1() != NULL) {
- edge_points.push_back(edge->vertex0()->vertex());
- edge_points.push_back(edge->vertex1()->vertex());
- fill_intermediate_points(point1, point2, point3, edge_points);
- } else {
- coordinate_type dir_x = (cell1->contains_segment() ^ (point1 == point3)) ?
- point3.y() - point2.y() : point2.y() - point3.y();
- coordinate_type dir_y = (cell1->contains_segment() ^ (point1 == point3)) ?
- point2.x() - point3.x() : point3.x() - point2.x();
- Point2D direction(dir_x, dir_y);
- find_intersections(point1, direction, LINE, brect, edge_points);
- if (edge->vertex1() != NULL)
- edge_points[1] = edge->vertex1()->vertex();
- if (edge->vertex0() != NULL)
- edge_points[0] = edge->vertex0()->vertex();
- }
+ // Append a new intermediate point to the end.
+ intermediate_points.push_back(
+ make_point_2d(new_point_x, new_point_y));
}
- return edge_points;
+
+ // Append last point to the end.
+ intermediate_points.push_back(last_point);
}
-
- private:
- voronoi_helper();
+ // Check whether extent is compatible with the edge type.
static bool check_extent(coordinate_type extent, kEdgeType etype) {
switch (etype) {
- case SEGMENT: return extent >= static_cast<coordinate_type>(0.0) &&
- extent <= static_cast<coordinate_type>(1.0);
+ case SEGMENT:
+ return extent >= static_cast<coordinate_type>(0.0) &&
+ extent <= static_cast<coordinate_type>(1.0);
case RAY: return extent >= static_cast<coordinate_type>(0.0);
case LINE: return true;
}
return true;
}
- static coordinate_type magnitude(coordinate_type value) {
+ // Compute the absolute value.
+ static inline coordinate_type magnitude(coordinate_type value) {
if (value >= static_cast<coordinate_type>(0.0))
return value;
return -value;
}
- static bool clip_line(coordinate_type denom,
- coordinate_type numer,
+ // Find fT min and max values: fT = numer / denom.
+ static void clip_line(coordinate_type denom, coordinate_type numer,
bool &fT0_used, bool &fT1_used,
- coordinate_type &fT0,
- coordinate_type &fT1) {
+ coordinate_type &fT0, coordinate_type &fT1) {
if (denom > static_cast<coordinate_type>(0.0)) {
if (fT1_used && numer > denom * fT1)
- return false;
+ return;
if (!fT0_used || numer > denom * fT0) {
fT0_used = true;
fT0 = numer / denom;
}
- return true;
} else if (denom < static_cast<coordinate_type>(0.0)) {
if (fT0_used && numer > denom * fT0)
- return false;
+ return;
if (!fT1_used || numer > denom * fT1) {
fT1_used = true;
fT1 = numer / denom;
}
- return true;
}
- return false;
}
- static coordinate_type get_point_projection(const Point2D &point,
- const Point2D &segment_start,
- const Point2D &segment_end) {
- coordinate_type segment_vec_x = segment_end.x() - segment_start.x();
- coordinate_type segment_vec_y = segment_end.y() - segment_start.y();
+ // 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 coordinate_type get_point_projection(
+ const point_2d_type &point,
+ const point_2d_type &segment_start,
+ const point_2d_type &segment_end) {
+ coordinate_type segment_vec_x = segment_end.x() -
+ segment_start.x();
+ coordinate_type segment_vec_y = segment_end.y() -
+ segment_start.y();
coordinate_type point_vec_x = point.x() - segment_start.x();
coordinate_type point_vec_y = point.y() - segment_start.y();
- coordinate_type sqr_segment_length = segment_vec_x * segment_vec_x +
- segment_vec_y * segment_vec_y;
+ coordinate_type sqr_segment_length =
+ segment_vec_x * segment_vec_x + segment_vec_y * segment_vec_y;
coordinate_type vec_dot = segment_vec_x * point_vec_x +
segment_vec_y * point_vec_y;
return vec_dot / sqr_segment_length;
}
- static int get_intermediate_points_count(const Point2D &point1,
- const Point2D &point2) {
+ // Return the number of the intermediate points to fill parabolic
+ // arc with. Arc is represented by it's two endpoints. This doesn't
+ // have any strong math background. Parameters were chosen
+ // experimentally.
+ static int get_intermediate_points_count(const point_2d_type &point1,
+ const point_2d_type &point2) {
coordinate_type vec_x = point1.x() - point2.x();
coordinate_type vec_y = point1.y() - point2.y();
coordinate_type sqr_len = vec_x * vec_x + vec_y * vec_y;
return static_cast<int>(log(sqr_len) * 3.4 + 1E-6);
}
};
-
+
+ // Represents voronoi cell.
+ // Data members: 1) pointer to the incident edge;
+ // 2) site inside the cell;
+ // 3) number of incident edges.
+ // The cell may contain point or segment site.
template <typename T>
class voronoi_cell {
public:
typedef T coordinate_type;
- typedef detail::site_event<T> site_event_type;
- typedef voronoi_edge<T> voronoi_edge_type;
-
- voronoi_cell(const site_event_type &new_site, voronoi_edge_type *edge) :
- site_(new_site),
+ typedef point_2d<coordinate_type> point_2d_type;
+ typedef detail::site_event<coordinate_type> site_event_type;
+ typedef voronoi_edge<coordinate_type> voronoi_edge_type;
+
+ voronoi_cell(const point_2d_type &p1, const point_2d_type &p2,
+ bool contains_segment, voronoi_edge_type *edge) :
+ point0_(p1),
+ point1_(p2),
+ contains_segment_(contains_segment),
incident_edge_(edge),
num_incident_edges_(0) {}
- bool contains_segment() const {
- return site_.is_segment();
- }
-
- const point_2d<T> &get_point0() const {
- return site_.get_point0();
- }
+ // Returns true if the cell contains segment site, false else.
+ bool contains_segment() const { return contains_segment_; }
- const point_2d<T> &get_point1() const {
- return site_.get_point1();
- }
-
- voronoi_edge_type *incident_edge() {
- return incident_edge_;
- }
+ // Returns site point in case cell contains point site,
+ // the first point of the segment site else.
+ const point_2d_type &point0() const { return point0_; }
+
+ // Returns the second site of the segment site.
+ // Don't use with the point sites.
+ const point_2d_type &point1() const { return point1_; }
+ voronoi_edge_type *incident_edge() { return incident_edge_; }
const voronoi_edge_type *incident_edge() const {
return incident_edge_;
}
- int num_incident_edges() const {
- return num_incident_edges_;
- }
+ int num_incident_edges() const { return num_incident_edges_; }
private:
- friend class voronoi_output<T>;
+ friend class voronoi_output<coordinate_type>;
void incident_edge(voronoi_edge_type *e) { incident_edge_ = e; }
void inc_num_incident_edges() { num_incident_edges_++; }
void dec_num_incident_edges() { num_incident_edges_--; }
-
- site_event_type site_;
+
+ point_2d_type point0_;
+ point_2d_type point1_;
+ bool contains_segment_;
voronoi_edge_type *incident_edge_;
int num_incident_edges_;
};
-
- template <typename T>
- struct robust_voronoi_vertex {
- detail::epsilon_robust_comparator<T> center_x;
- detail::epsilon_robust_comparator<T> center_y;
- detail::epsilon_robust_comparator<T> denom;
-
- robust_voronoi_vertex(const detail::epsilon_robust_comparator<T> &c_x,
- const detail::epsilon_robust_comparator<T> &c_y,
- const detail::epsilon_robust_comparator<T> &den) :
- center_x(c_x),
- center_y(c_y),
- denom(den) {}
- };
+ // Represents voronoi vertex.
+ // Data members: 1) robust vertex data structure;
+ // 2) vertex point itself;
+ // 3) pointer to the incident edge;
+ // 4) number of incident edges.
template <typename T>
class voronoi_vertex {
public:
typedef T coordinate_type;
- typedef voronoi_edge<T> voronoi_edge_type;
+ typedef point_2d<T> point_2d_type;
+ typedef voronoi_edge<coordinate_type> voronoi_edge_type;
- voronoi_vertex(robust_voronoi_vertex<T> *robust_vertex,
- voronoi_edge_type *edge) :
- robust_vertex_(robust_vertex),
- vertex_(robust_vertex->center_x.dif() / robust_vertex->denom.dif(),
- robust_vertex->center_y.dif() / robust_vertex->denom.dif()),
+ voronoi_vertex(const point_2d_type &vertex, voronoi_edge_type *edge) :
+ vertex_(vertex),
incident_edge_(edge),
num_incident_edges_(3) {}
- const point_2d<T> &vertex() const {
- return vertex_;
- }
-
- const robust_voronoi_vertex<T> *robust_vertex() const {
- return robust_vertex_;
- }
+ const point_2d_type &vertex() const { return vertex_; }
- voronoi_edge_type *incident_edge() {
+ voronoi_edge_type *incident_edge() { return incident_edge_; }
+ const voronoi_edge_type *incident_edge() const {
return incident_edge_;
}
- const voronoi_edge_type *incident_edge() const {
- return incident_edge_;
+ int num_incident_edges() const { return num_incident_edges_; }
+
+ private:
+ typedef detail::robust_voronoi_vertex<coordinate_type>
+ robust_voronoi_vertex_type;
+
+ friend class voronoi_output<coordinate_type>;
+
+ const robust_voronoi_vertex_type *robust_vertex() const {
+ return robust_vertex_;
}
- int num_incident_edges() const {
- return num_incident_edges_;
+ void robust_voronoi_vertex(robust_voronoi_vertex_type *v) {
+ robust_vertex_ = v;
}
- private:
- friend class voronoi_output<T>;
-
void incident_edge(voronoi_edge_type *e) { incident_edge_ = e; }
void num_incident_edges(int n) { num_incident_edges_ = n; }
- robust_voronoi_vertex<T> *robust_vertex_;
- point_2d<T> vertex_;
+ robust_voronoi_vertex_type *robust_vertex_;
+ point_2d_type vertex_;
voronoi_edge_type *incident_edge_;
int num_incident_edges_;
};
// Half-edge data structure. Represents voronoi edge.
- // Contains: 1) pointer to cell records;
- // 2) pointers to start/end vertices of half-edge;
- // 3) pointers to previous/next half-edges(CCW);
- // 4) pointers to previous/next half-edges rotated
- // around start point;
- // 5) pointer to twin half-edge.
+ // Variables: 1) pointer to the corresponding cell;
+ // 2) pointer to the vertex that is the starting
+ // point of the half-edge;
+ // 3) pointer to the twin edge;
+ // 4) pointer to the CCW/CW next edge.
+ // 5) pointer to the CCW/CW prev edge.
template <typename T>
class voronoi_edge {
public:
- typedef voronoi_cell<T> voronoi_cell_type;
- typedef voronoi_vertex<T> voronoi_vertex_type;
- typedef voronoi_edge<T> voronoi_edge_type;
+ typedef T coordinate_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_edge() :
cell_(NULL),
@@ -555,37 +685,51 @@
voronoi_edge_type *prev() { return prev_; }
const voronoi_edge_type *prev() const { return prev_; }
+ // Return a pointer to the rotation next edge
+ // over the starting point of the half-edge.
voronoi_edge_type *rot_next() {
- if (vertex_)
- return prev_->twin();
- return NULL;
+ return (vertex_) ? prev_->twin() : NULL;
}
const voronoi_edge_type *rot_next() const {
- if (vertex_)
- return prev_->twin();
- return NULL;
+ return (vertex_) ? prev_->twin() : NULL;
}
+ // Return a pointer to the rotation prev edge
+ // over the starting point of the half-edge.
voronoi_edge_type *rot_prev() {
- if (vertex_)
- return twin_->next();
- return NULL;
+ return (vertex_) ? twin_->next() : NULL;
}
const voronoi_edge_type *rot_prev() const {
- if (vertex_)
- return twin_->next();
- return NULL;
+ return (vertex_) ? twin_->next() : NULL;
+ }
+
+ // Return true if the edge is finite (segment, parabolic arc).
+ // Return false if the edge is infinite (ray, line).
+ bool is_bounded() const { return vertex0() && vertex1(); }
+
+ // Return true if the edge is linear.
+ // Return false if the edge is curved (parabolic arc).
+ bool is_linear() const {
+ return !(cell()->contains_segment() ^
+ twin()->cell()->contains_segment());
+ }
+
+ // Returns true if the edge is curved (parabolic arc).
+ // Returns false if the edge is linear.
+ bool is_curved() const {
+ return (cell()->contains_segment() ^
+ twin()->cell()->contains_segment());
}
private:
- friend class voronoi_output<T>;
+ friend class voronoi_output<coordinate_type>;
void cell(voronoi_cell_type *c) { cell_ = c; }
void vertex0(voronoi_vertex_type *v) { vertex_ = v; }
void vertex1(voronoi_vertex_type *v) { twin_->vertex0(v); }
void twin(voronoi_edge_type *e) { twin_ = e; }
void next(voronoi_edge_type *e) { next_ = e; }
- void prev(voronoi_edge_type *e) { prev_ = e; }
+ void prev(voronoi_edge_type *e) { prev_ = e; }
voronoi_cell_type *cell_;
voronoi_vertex_type *vertex_;
@@ -594,37 +738,64 @@
voronoi_edge_type *prev_;
};
- // Voronoi output data structure based on the half-edges.
- // Contains vector of voronoi cells, doubly linked list of
- // voronoi vertices and voronoi edges.
+ // Voronoi output data structure.
+ // Data members:
+ // 1) cell_records_ - vector of the voronoi cells;
+ // 2) vertex_records_ - list of the voronoi vertices;
+ // 3) edge_records_ - list of the voronoi edges;
+ // 4) robust_vertices_ - list of the robust vertices;
+ // 5) voronoi_rect_ - bounding rectangle;
+ // 6) num_cell_records_ - number of the voronoi cells;
+ // 7) num_vertex_records_ - number of the voronoi vertices;
+ // 8) num_edge_records_ - number of the voronoi edges.
+ // CCW ordering is used on the faces perimeter and around the vertices.
+ // Robust vertices are used to make the simplification stage epsilon
+ // robust. Vector data structure is used to store voronoi cells as the
+ // number of the cells may be precomputed at the initialization step.
+ // As size() method takes O(n) time on the list data structure three
+ // additional counters are used to count the number of the voronoi cells,
+ // vertices and edges. As we use list data structure to represent voronoi
+ // vertices and edges there is no copy method available, because it will
+ // invalidate all the pointers. Another approach could be used to make
+ // copying available:
+ // 1) use vectors to store voronoi vertices and cells;
+ // 2) store vector indices instead of the pointers;
+ // 3) store additional pointer to the voronoi output data structure
+ // in the voronoi cell, vertex, edge data structure.
+ // 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>
class voronoi_output {
public:
typedef T coordinate_type;
- typedef point_2d<T> Point2D;
- typedef detail::site_event<coordinate_type> site_event_type;
- typedef detail::circle_event<coordinate_type> circle_event_type;
+ typedef point_2d<coordinate_type> point_2d_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 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 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;
- typedef typename voronoi_edges_type::const_iterator voronoi_edge_const_iterator_type;
-
- voronoi_output() {
- num_cell_records_ = 0;
- num_edge_records_ = 0;
- num_vertex_records_ = 0;
- }
+ typedef typename voronoi_edges_type::iterator
+ voronoi_edge_iterator_type;
+ typedef typename voronoi_edges_type::const_iterator
+ voronoi_edge_const_iterator_type;
+
+ voronoi_output() :
+ num_cell_records_(0),
+ num_edge_records_(0),
+ num_vertex_records_(0) {}
void clear() {
robust_vertices_.clear();
@@ -666,108 +837,137 @@
}
private:
+ typedef detail::site_event<coordinate_type> site_event_type;
+ typedef detail::circle_event<coordinate_type> circle_event_type;
+ typedef detail::robust_voronoi_vertex<coordinate_type>
+ robust_voronoi_vertex_type;
+
friend class detail::voronoi_builder<T>;
void init(int num_sites) {
+ // Resize cell vector to avoid reallocations.
cell_records_.reserve(num_sites);
+
+ // Init counters.
num_cell_records_ = 0;
num_edge_records_ = 0;
num_vertex_records_ = 0;
}
- // Update voronoi output in case of single site input.
+ // Update the voronoi output in case of a single point input.
void process_single_site(const site_event_type &site) {
// Update bounding rectangle.
- voronoi_rect_ = BRect<coordinate_type>(site.get_point0(), site.get_point0());
+ voronoi_rect_ = BRect<coordinate_type>(site.point0(),
+ site.point0());
// Update cell records.
- cell_records_.push_back(voronoi_cell_type(site, NULL));
+ cell_records_.push_back(voronoi_cell_type(site.point0(),
+ site.point1(),
+ site.is_segment(),
+ NULL));
}
- // Inserts new half-edge into the output data structure during site
- // event processing. Takes as input left and right sites of the new
- // beach line node and returns pointer to the new half-edge.
+ // Insert a new half-edge into the output data structure.
+ // Takes as input left and right sites that form a new bisector.
+ // Returns a pointer to a new half-edge.
voronoi_edge_type *insert_new_edge(const site_event_type &site1,
const site_event_type &site2) {
- // Get indices of sites.
- int site_index1 = site1.get_site_index();
- int site_index2 = site2.get_site_index();
+ // Get sites' indices.
+ int site_index1 = site1.index();
+ int site_index2 = site2.index();
- // Create new half-edge that belongs to the first site.
+ // Create a new half-edge that belongs to the first site.
edge_records_.push_back(voronoi_edge_type());
voronoi_edge_type &edge1 = edge_records_.back();
- // Create new half-edge that belongs to the second site.
+ // Create a new half-edge that belongs to the second site.
edge_records_.push_back(voronoi_edge_type());
voronoi_edge_type &edge2 = edge_records_.back();
- // Add initial cell during first edge insertion.
+ // Add the initial cell during the first edge insertion.
if (cell_records_.empty()) {
- cell_records_.push_back(voronoi_cell_type(site1, &edge1));
- voronoi_rect_ = BRect<coordinate_type>(site1.get_point0(), site1.get_point0());
+ cell_records_.push_back(voronoi_cell_type(site1.point0(),
+ site1.point1(),
+ site1.is_segment(),
+ &edge1));
+ voronoi_rect_ = BRect<coordinate_type>(site1.point0(),
+ site1.point0());
}
cell_records_[site_index1].inc_num_incident_edges();
- // Update bounding rectangle.
- voronoi_rect_.update(site2.get_point0());
+ // Update the bounding rectangle.
+ voronoi_rect_.update(site2.point0());
if (site2.is_segment()) {
- voronoi_rect_.update(site2.get_point1());
+ voronoi_rect_.update(site2.point1());
}
- // Second site represents new site during site event processing.
- // Add new cell to the cell records vector.
- cell_records_.push_back(voronoi_cell_type(site2, &edge2));
+ // 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(),
+ site2.is_segment(),
+ &edge2));
cell_records_.back().inc_num_incident_edges();
-
- // Update pointers to cells.
+
+ // Set up pointers to cells.
edge1.cell(&cell_records_[site_index1]);
edge2.cell(&cell_records_[site_index2]);
- // Update twin pointers.
+ // Set up twin pointers.
edge1.twin(&edge2);
edge2.twin(&edge1);
+ // Return a pointer to the new half-edge.
return &edge1;
}
+ // Insert a new half-edge into the output data structure with the
+ // start at the point where two previously added half-edges intersect.
+ // Takes as input two sites that create a new bisector, circle event
+ // that correponds to the intersection point of the two old half-edges,
+ // pointers to those half-edges. Half-edges' direction goes out of the
+ // new voronoi vertex point. Returns a pointer to the new half-edge.
voronoi_edge_type *insert_new_edge(const site_event_type &site1,
const site_event_type &site3,
const circle_event_type &circle,
voronoi_edge_type *edge12,
voronoi_edge_type *edge23) {
- // Update bounding rectangle.
- //voronoi_rect_.update(circle.get_center());
-
- // Add new voronoi vertex with point at center of the circle.
- robust_vertices_.push_back(robust_voronoi_vertex<T>(
- circle.get_erc_x(), circle.get_erc_y(), circle.get_erc_denom()));
- vertex_records_.push_back(voronoi_vertex_type(&robust_vertices_.back(), edge12));
+ // Add a new voronoi vertex.
+ robust_vertices_.push_back(
+ robust_voronoi_vertex_type(circle.erc_x(),
+ circle.erc_y(),
+ circle.erc_denom()));
+ const robust_voronoi_vertex_type &robust_vertex =
+ robust_vertices_.back();
+ vertex_records_.push_back(voronoi_vertex_type(
+ make_point_2d(robust_vertex.x(), robust_vertex.y()), edge12));
voronoi_vertex_type &new_vertex = vertex_records_.back();
+ new_vertex.robust_voronoi_vertex(&robust_vertices_.back());
- // Update two input bisectors and their twin half-edges with
- // new voronoi vertex.
+ // Update vertex pointers of the old edges.
edge12->vertex0(&new_vertex);
edge23->vertex0(&new_vertex);
- // Add new half-edge.
+ // Add a new half-edge.
edge_records_.push_back(voronoi_edge_type());
voronoi_edge_type &new_edge1 = edge_records_.back();
- new_edge1.cell(&cell_records_[site1.get_site_index()]);
+ new_edge1.cell(&cell_records_[site1.index()]);
new_edge1.cell()->inc_num_incident_edges();
- // Add new half-edge.
+ // Add a new half-edge.
edge_records_.push_back(voronoi_edge_type());
voronoi_edge_type &new_edge2 = edge_records_.back();
- new_edge2.cell(&cell_records_[site3.get_site_index()]);
+ new_edge2.cell(&cell_records_[site3.index()]);
new_edge2.cell()->inc_num_incident_edges();
- // Update twin pointers of the new half-edges.
+ // Update twin pointers.
new_edge1.twin(&new_edge2);
new_edge2.twin(&new_edge1);
+
+ // Update vertex pointer.
new_edge2.vertex0(&new_vertex);
- // Update voronoi prev/next pointers of all half-edges incident
- // to the new voronoi vertex.
+ // Update voronoi prev/next pointers.
edge12->prev(&new_edge1);
new_edge1.next(edge12);
edge12->twin()->next(edge23);
@@ -775,20 +975,28 @@
edge23->twin()->next(&new_edge2);
new_edge2.prev(edge23->twin());
+ // Return a pointer to the new half-edge.
return &new_edge1;
}
+ // Remove zero-length edges from the voronoi output.
void simplify() {
voronoi_edge_iterator_type edge_it1;
voronoi_edge_iterator_type edge_it = edge_records_.begin();
num_cell_records_ = cell_records_.size();
- // Return in case of collinear sites input.
+ // All the initial sites are colinear.
if (vertex_records_.empty()) {
+ // Update edges counter.
num_edge_records_ = num_cell_records_ - 1;
+
+ // Return if there are no edges.
if (edge_it == edge_records_.end())
return;
+ // Update prev/next pointers of the edges. Those edges won't
+ // have any common endpoints, cause they are infinite.
+ // But still they follow each other using CCW ordering.
voronoi_edge_type *edge1 = &(*edge_it);
edge1->next(edge1);
edge1->prev(edge1);
@@ -799,7 +1007,7 @@
while (edge_it != edge_records_.end()) {
voronoi_edge_type *edge2 = &(*edge_it);
edge_it++;
-
+
edge1->next(edge2);
edge1->prev(edge2);
edge2->next(edge1);
@@ -814,61 +1022,75 @@
return;
}
- // Iterate over all edges and remove degeneracies.
+ // Iterate over all the edges and remove degeneracies
+ // (zero-length edges). Edge is considered to be zero-length
+ // if both its endpoints lie within some epsilon-rectangle.
while (edge_it != edge_records_.end()) {
edge_it1 = edge_it;
std::advance(edge_it, 2);
+ // Degenerate edges exist only among finite edges.
if (!edge_it1->vertex0() || !edge_it1->vertex1()) {
num_edge_records_++;
continue;
}
- const robust_voronoi_vertex<T> *v1 = edge_it1->vertex0()->robust_vertex();
- const robust_voronoi_vertex<T> *v2 = edge_it1->vertex1()->robust_vertex();
- detail::epsilon_robust_comparator<T> lhs1(v1->center_x * v2->denom);
- detail::epsilon_robust_comparator<T> rhs1(v1->denom * v2->center_x);
- detail::epsilon_robust_comparator<T> lhs2(v1->center_y * v2->denom);
- detail::epsilon_robust_comparator<T> rhs2(v1->denom * v2->center_y);
+ const voronoi_vertex_type *v1 = edge_it1->vertex0();
+ const voronoi_vertex_type *v2 = edge_it1->vertex1();
- if (lhs1.compares_undefined(rhs1, 64) && lhs2.compares_undefined(rhs2, 64)) {
- // Decrease number of cell incident edges.
+ // Make epsilon robust check.
+ if (v1->robust_vertex()->equals(v2->robust_vertex(), 64)) {
+ // Decrease number of cell's incident edges.
edge_it1->cell()->dec_num_incident_edges();
edge_it1->twin()->cell()->dec_num_incident_edges();
- // To guarantee N*lon(N) time we merge vertex with
- // less incident edges to the one with more.
+ // Corresponding cell incident edge pointer
+ // points to the degenerate edge.
if (edge_it1->cell()->incident_edge() == &(*edge_it1)) {
- if (edge_it1->cell()->incident_edge() == edge_it1->next()) {
+ // Update incident edge pointer.
+ if (edge_it1->cell()->incident_edge() ==
+ edge_it1->next()) {
edge_it1->cell()->incident_edge(NULL);
} else {
edge_it1->cell()->incident_edge(edge_it1->next());
}
}
- if (edge_it1->twin()->cell()->incident_edge() == edge_it1->twin()) {
- if (edge_it1->twin()->cell()->incident_edge() == edge_it1->twin()->next()) {
+
+ // Cell corresponding to the twin edge
+ // points to the degenerate edge.
+ if (edge_it1->twin()->cell()->incident_edge() ==
+ edge_it1->twin()) {
+ // Update incident edge pointer.
+ if (edge_it1->twin()->cell()->incident_edge() ==
+ edge_it1->twin()->next()) {
edge_it1->twin()->cell()->incident_edge(NULL);
} else {
- edge_it1->twin()->cell()->incident_edge(edge_it1->twin()->next());
+ edge_it1->twin()->cell()->incident_edge(
+ edge_it1->twin()->next());
}
}
- if (edge_it1->vertex0()->num_incident_edges() >=
- edge_it1->vertex1()->num_incident_edges()) {
- simplify_edge(&(*edge_it1));
+
+ // To guarantee N*lon(N) time 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 {
- simplify_edge(edge_it1->twin());
+ remove_edge(edge_it1->twin());
}
- // Remove zero length edges.
+ // Remove zero-length edge.
edge_records_.erase(edge_it1, edge_it);
} else {
+ // Count not degenerate edge.
num_edge_records_++;
}
}
robust_vertices_.clear();
- for (voronoi_vertex_iterator_type vertex_it = vertex_records_.begin();
- vertex_it != vertex_records_.end();) {
+ // Remove degenerate voronoi vertices with zero incident edges.
+ for (voronoi_vertex_iterator_type vertex_it =
+ vertex_records_.begin();
+ vertex_it != vertex_records_.end();) {
if (vertex_it->incident_edge() == NULL) {
vertex_it = vertex_records_.erase(vertex_it);
} else {
@@ -877,10 +1099,11 @@
}
}
- // Make some post processing.
+ // Update prev/next pointers for the ray-edges.
for (voronoi_cell_iterator_type cell_it = cell_records_.begin();
- cell_it != cell_records_.end(); cell_it++) {
- // Move to the previous edge while it is possible in the CW direction.
+ cell_it != cell_records_.end(); cell_it++) {
+ // Move to the previous edge while
+ // it is possible in the CW direction.
voronoi_edge_type *cur_edge = cell_it->incident_edge();
if (cur_edge) {
while (cur_edge->prev() != NULL) {
@@ -891,12 +1114,14 @@
break;
}
- // Rewind incident edge pointer to the leftmost edge for the boundary cells.
+ // Rewind incident edge pointer to the
+ // CW leftmost edge for the boundary cells.
cell_it->incident_edge(cur_edge);
- // Set up prev/next half-edge pointers for ray edges.
+ // Set up prev/next half-edge pointers for the ray-edges.
if (cur_edge->prev() == NULL) {
- voronoi_edge_type *last_edge = cell_it->incident_edge();
+ voronoi_edge_type *last_edge =
+ cell_it->incident_edge();
while (last_edge->next() != NULL)
last_edge = last_edge->next();
last_edge->next(cur_edge);
@@ -906,16 +1131,16 @@
}
}
- // Simplify degenerate edge.
- void simplify_edge(voronoi_edge_type *edge) {
+ // Remove degenerate edge.
+ void remove_edge(voronoi_edge_type *edge) {
voronoi_vertex_type *vertex1 = edge->vertex0();
voronoi_vertex_type *vertex2 = edge->vertex1();
// Update number of incident edges.
- vertex1->num_incident_edges(
- vertex1->num_incident_edges() + vertex2->num_incident_edges() - 2);
+ vertex1->num_incident_edges(vertex1->num_incident_edges() +
+ vertex2->num_incident_edges() - 2);
- // Update second vertex incident edges start and end points.
+ // Update the endpoints of the incident edges to the second vertex.
voronoi_edge_type *updated_edge = edge->twin()->rot_next();
while (updated_edge != edge->twin()) {
updated_edge->vertex0(vertex1);
@@ -931,32 +1156,32 @@
voronoi_edge_type *edge2_rot_prev = edge2->rot_prev();
voronoi_edge_type *edge2_rot_next = edge2->rot_next();
- // Update prev and next pointers of incident edges.
+ // Update prev/next pointers for the incident edges.
edge1_rot_next->twin()->next(edge2_rot_prev);
edge2_rot_prev->prev(edge1_rot_next->twin());
edge1_rot_prev->prev(edge2_rot_next->twin());
edge2_rot_next->twin()->next(edge1_rot_prev);
- // Change incident edge pointer of a vertex if it corresponds to the
+ // Change the pointer to the incident edge if it points to the
// degenerate edge.
if (vertex1->incident_edge() == edge) {
vertex1->incident_edge(edge->rot_prev());
}
- // Remove second vertex from the vertex records list.
+ // Set the incident edge point of the second vertex to NULL value.
if (vertex1 != vertex2) {
vertex2->incident_edge(NULL);
}
}
- std::list< robust_voronoi_vertex<T> > robust_vertices_;
+ std::list< robust_voronoi_vertex_type > robust_vertices_;
voronoi_cells_type cell_records_;
voronoi_vertices_type vertex_records_;
voronoi_edges_type edge_records_;
int num_cell_records_;
- int num_vertex_records_;
int num_edge_records_;
+ int num_vertex_records_;
BRect<coordinate_type> voronoi_rect_;
Modified: sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_sweepline.hpp
==============================================================================
--- sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_sweepline.hpp (original)
+++ sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_sweepline.hpp 2010-12-20 09:30:57 EST (Mon, 20 Dec 2010)
@@ -1,11 +1,11 @@
-// Boost sweepline/voronoi_sweepline.hpp header file
+// Boost sweepline/voronoi_sweepline.hpp header file
// Copyright Andrii Sydorchuk 2010.
// 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.
+// See http://www.boost.org for updates, documentation, and revision history.
#ifndef BOOST_SWEEPLINE_VORONOI_SWEEPLINE
#define BOOST_SWEEPLINE_VORONOI_SWEEPLINE
@@ -18,7 +18,7 @@
#include <queue>
#include <vector>
-#ifdef USE_MULTI_PRECISION_LIBRARY
+#ifdef USE_MULTI_PRECISION_LIBRARY
#pragma warning( disable : 4800 )
#include "gmpxx.h"
#endif
@@ -29,24 +29,37 @@
namespace boost {
namespace sweepline {
+ // Public methods to compute voronoi diagram.
+ // points - vector of input points.
+ // segments - vector of input segments.
+ // output - voronoi output data structure to hold voronoi diagram.
+ // The assumption is made that input doesn't contain segments that
+ // intersect or points lying on the segments. Also coordinates of
+ // the points and of the endpoints of the segments should be within
+ // signed integer range [-2^31, 2^31-1].
+ // Complexity - O(N*logN), memory usage - O(N), where N is the total
+ // number of points and segments.
+
template <typename T>
- static void build_voronoi(const std::vector< point_2d<T> > &points,
- voronoi_output<double> &output) {
+ static inline void build_voronoi(const std::vector< point_2d<T> > &points,
+ voronoi_output<double> &output) {
std::vector< std::pair< point_2d<T>, point_2d<T> > > segments_empty;
build_voronoi<T>(points, segments_empty, output);
}
template <typename T>
- static void build_voronoi(const std::vector< std::pair< point_2d<T>, point_2d<T> > > &segments,
- voronoi_output<double> &output) {
+ static inline void build_voronoi(
+ const std::vector< std::pair< point_2d<T>, point_2d<T> > > &segments,
+ voronoi_output<double> &output) {
std::vector< point_2d<T> > points_empty;
build_voronoi<T>(points_empty, segments, output);
}
template <typename T>
- static void build_voronoi(const std::vector< point_2d<T> > &points,
- const std::vector< std::pair< point_2d<T>, point_2d<T> > > &segments,
- voronoi_output<double> &output) {
+ static inline void build_voronoi(
+ const std::vector< point_2d<T> > &points,
+ const std::vector< std::pair< point_2d<T>, point_2d<T> > > &segments,
+ voronoi_output<double> &output) {
detail::voronoi_builder<double> builder(output);
builder.init(points, segments);
builder.run_sweepline();
Modified: sandbox/SOC/2010/sweepline/libs/sweepline/example/voronoi_visualizer.cpp
==============================================================================
--- sandbox/SOC/2010/sweepline/libs/sweepline/example/voronoi_visualizer.cpp (original)
+++ sandbox/SOC/2010/sweepline/libs/sweepline/example/voronoi_visualizer.cpp 2010-12-20 09:30:57 EST (Mon, 20 Dec 2010)
@@ -25,15 +25,15 @@
}
void build(QString file_path) {
- std::vector<iPoint2D> point_sites;
- std::vector<iSegment2D> segment_sites;
+ std::vector<ipoint_2d_type> point_sites;
+ std::vector<isegment_2d_type> segment_sites;
// Open file.
QFile data(file_path);
if (!data.open(QFile::ReadOnly))
QMessageBox::warning(this, tr("Voronoi Visualizer"),
tr("Disable to open file ") + file_path);
-
+
// Read points from the file.
QTextStream in_stream(&data);
int num_point_sites = 0;
@@ -75,7 +75,7 @@
glBegin(GL_POINTS);
for (it = cells.begin(); it != cells.end(); it++) {
if (!it->contains_segment())
- glVertex2f(it->get_point0().x(), it->get_point0().y());
+ glVertex2f(it->point0().x(), it->point0().y());
}
glEnd();
glPointSize(6);
@@ -83,8 +83,8 @@
glBegin(GL_LINES);
for (it = cells.begin(); it != cells.end(); it++) {
if (it->contains_segment()) {
- glVertex2f(it->get_point0().x(), it->get_point0().y());
- glVertex2f(it->get_point1().x(), it->get_point1().y());
+ glVertex2f(it->point0().x(), it->point0().y());
+ glVertex2f(it->point1().x(), it->point1().y());
}
}
glEnd();
@@ -109,7 +109,7 @@
glColor3f(0.0f, 1.0f, 0.0f);
glBegin(GL_LINES);
for (it = edges.begin(); it != edges.end(); it++) {
- std::vector<Point2D> temp_v =
+ std::vector<point_2d_type> temp_v =
voronoi_helper<coordinate_type>::get_intermediate_points(&(*it), brect_);
for (int i = 0; i < static_cast<int>(temp_v.size()) - 1; i++) {
glVertex2f(temp_v[i].x(), temp_v[i].y());
@@ -138,9 +138,9 @@
}
typedef double coordinate_type;
- typedef point_2d<coordinate_type> Point2D;
- typedef point_2d<int> iPoint2D;
- typedef std::pair<iPoint2D, iPoint2D> iSegment2D;
+ typedef point_2d<coordinate_type> point_2d_type;
+ typedef point_2d<int> ipoint_2d_type;
+ typedef std::pair<ipoint_2d_type, ipoint_2d_type> isegment_2d_type;
typedef voronoi_output<coordinate_type>::voronoi_cells_type voronoi_cells_type;
typedef voronoi_output<coordinate_type>::voronoi_vertices_type voronoi_vertices_type;
typedef voronoi_output<coordinate_type>::voronoi_edges_type voronoi_edges_type;
@@ -161,12 +161,12 @@
glWidget_ = new GLWidget();
file_dir_ = QDir(QDir::currentPath(), tr("*.txt"));
file_name_ = tr("");
-
+
QHBoxLayout *centralLayout = new QHBoxLayout;
centralLayout->addWidget(glWidget_);
centralLayout->addLayout(create_file_layout());
setLayout(centralLayout);
-
+
update_file_list();
setWindowTitle(tr("Voronoi Visualizer"));
layout()->setSizeConstraint( QLayout::SetFixedSize );
@@ -203,9 +203,9 @@
private:
QGridLayout *create_file_layout() {
QGridLayout *file_layout = new QGridLayout;
-
+
message_label_ = new QLabel("Double click item to build voronoi diagram:");
-
+
file_list_ = new QListWidget();
file_list_->connect(file_list_, SIGNAL(itemDoubleClicked(QListWidgetItem*)),
this, SLOT(build()));
Modified: sandbox/SOC/2010/sweepline/libs/sweepline/test/benchmark_test.cpp
==============================================================================
--- sandbox/SOC/2010/sweepline/libs/sweepline/test/benchmark_test.cpp (original)
+++ sandbox/SOC/2010/sweepline/libs/sweepline/test/benchmark_test.cpp 2010-12-20 09:30:57 EST (Mon, 20 Dec 2010)
@@ -5,7 +5,7 @@
// (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.
+// See http://www.boost.org for updates, documentation, and revision history.
#include <iostream>
#include <stdio.h>
Modified: sandbox/SOC/2010/sweepline/libs/sweepline/test/circle_event_creation_test.cpp
==============================================================================
--- sandbox/SOC/2010/sweepline/libs/sweepline/test/circle_event_creation_test.cpp (original)
+++ sandbox/SOC/2010/sweepline/libs/sweepline/test/circle_event_creation_test.cpp 2010-12-20 09:30:57 EST (Mon, 20 Dec 2010)
@@ -1,11 +1,11 @@
-// Boost sweepline library circle_event_creation_test.cpp file
+// Boost sweepline library circle_event_creation_test.cpp file
// Copyright Andrii Sydorchuk 2010.
// 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.
+// See http://www.boost.org for updates, documentation, and revision history.
#include "test_type_list.hpp"
#include "boost/sweepline/voronoi_sweepline.hpp"
@@ -41,7 +41,7 @@
site_event_type site2_1(static_cast<T>(min_int), static_cast<T>(max_int), 1);
site_event_type site2_2(static_cast<T>(max_int-1), static_cast<T>(max_int-1), 1);
site_event_type site3(static_cast<T>(max_int), static_cast<T>(max_int), 2);
- typename detail::circle_event<T> c_event;
+ typename detail::circle_event<T> c_event;
bool is_event = detail::create_circle_event_ppp(site1, site2_1, site3, c_event);
BOOST_CHECK_EQUAL(is_event, true);
is_event = detail::create_circle_event_ppp(site1, site2_2, site3, c_event);
@@ -53,7 +53,7 @@
typedef typename detail::site_event<T> site_event_type;
site_event_type segm_start(static_cast<T>(-4), static_cast<T>(0), 0);
site_event_type segm_end(static_cast<T>(0), static_cast<T>(4), 0);
- site_event_type segm_site(segm_start.get_point0(), segm_end.get_point0(), 0);
+ site_event_type segm_site(segm_start.point0(), segm_end.point0(), 0);
typename detail::circle_event<T> c_event;
bool is_event = detail::create_circle_event_pps(segm_start, segm_end, segm_site, 2, c_event);
BOOST_CHECK_EQUAL(is_event, false);
@@ -61,8 +61,8 @@
site_event_type site1_2(static_cast<T>(0), static_cast<T>(2), 0);
is_event = detail::create_circle_event_pps(site1_1, site1_2, segm_site, 2, c_event);
BOOST_CHECK_EQUAL(is_event, true);
- BOOST_CHECK_EQUAL(c_event.get_center().x() == static_cast<T>(-1), true);
- BOOST_CHECK_EQUAL(c_event.get_center().y() == static_cast<T>(1), true);
+ BOOST_CHECK_EQUAL(c_event.center().x() == static_cast<T>(-1), true);
+ BOOST_CHECK_EQUAL(c_event.center().y() == static_cast<T>(1), true);
is_event = detail::create_circle_event_pps(site1_1, site1_2, segm_site, 1, c_event);
BOOST_CHECK_EQUAL(is_event, false);
is_event = detail::create_circle_event_pps(site1_1, site1_2, segm_site, 3, c_event);
@@ -73,7 +73,7 @@
typedef typename detail::site_event<T> site_event_type;
site_event_type segm_start(static_cast<T>(-4), static_cast<T>(0), 0);
site_event_type segm_end(static_cast<T>(-4), static_cast<T>(20), 0);
- site_event_type segm_site(segm_start.get_point0(), segm_end.get_point0(), 0);
+ site_event_type segm_site(segm_start.point0(), segm_end.point0(), 0);
typename detail::circle_event<T> c_event;
site_event_type site1_1(static_cast<T>(-2), static_cast<T>(10), 0);
site_event_type site1_2(static_cast<T>(4), static_cast<T>(10), 0);
@@ -91,7 +91,7 @@
point_2d<T> segm1_start(static_cast<T>(1), static_cast<T>(0));
point_2d<T> segm1_end(static_cast<T>(7), static_cast<T>(0));
site_event_type segm_site1(segm1_start, segm1_end, 0);
- segm_site1.set_inverse();
+ segm_site1.inverse();
point_2d<T> segm2_start(static_cast<T>(-2), static_cast<T>(4));
point_2d<T> segm2_end(static_cast<T>(10), static_cast<T>(4));
site_event_type segm_site2(segm2_start, segm2_end, 1);
@@ -113,7 +113,7 @@
point_2d<T> segm1_start(static_cast<T>(-1), static_cast<T>(2));
point_2d<T> segm1_end(static_cast<T>(8), static_cast<T>(-10));
site_event_type segm_site1(segm1_start, segm1_end, 0);
- segm_site1.set_inverse();
+ segm_site1.inverse();
point_2d<T> segm2_start(static_cast<T>(-1), static_cast<T>(0));
point_2d<T> segm2_end(static_cast<T>(8), static_cast<T>(12));
site_event_type segm_site2(segm2_start, segm2_end, 1);
@@ -129,7 +129,7 @@
point_2d<T> segm1_start(static_cast<T>(1), static_cast<T>(0));
point_2d<T> segm1_end(static_cast<T>(6), static_cast<T>(0));
site_event_type segm_site1(segm1_start, segm1_end, 0);
- segm_site1.set_inverse();
+ segm_site1.inverse();
point_2d<T> segm2_start(static_cast<T>(-6), static_cast<T>(4));
point_2d<T> segm2_end(static_cast<T>(0), static_cast<T>(12));
site_event_type segm_site2(segm2_start, segm2_end, 1);
@@ -147,7 +147,7 @@
point_2d<T> point3(static_cast<T>(8), static_cast<T>(6));
point_2d<T> point4(static_cast<T>(0), static_cast<T>(12));
site_event_type segm_site1(point1, point2, 0);
- segm_site1.set_inverse();
+ segm_site1.inverse();
site_event_type segm_site2(point3, point4, 1);
site_event_type point_site(point1, 2);
typename detail::circle_event<T> c_event;
@@ -165,7 +165,7 @@
point_2d<T> point3(static_cast<T>(0), static_cast<T>(4));
point_2d<T> point4(static_cast<T>(4), static_cast<T>(4));
site_event_type segm_site1(point1, point2, 0);
- segm_site1.set_inverse();
+ segm_site1.inverse();
site_event_type segm_site2(point2, point3, 1);
site_event_type segm_site3(point3, point4, 2);
typename detail::circle_event<T> c_event;
@@ -181,7 +181,7 @@
point_2d<T> point3(static_cast<T>(-39), static_cast<T>(30));
point_2d<T> point4(static_cast<T>(1), static_cast<T>(60));
site_event_type segm_site1(point1, point2, 0);
- segm_site1.set_inverse();
+ segm_site1.inverse();
site_event_type segm_site2(point3, point4, 1);
site_event_type segm_site3(point4, point1, 2);
typename detail::circle_event<T> c_event;
Modified: sandbox/SOC/2010/sweepline/libs/sweepline/test/clipping_test.cpp
==============================================================================
--- sandbox/SOC/2010/sweepline/libs/sweepline/test/clipping_test.cpp (original)
+++ sandbox/SOC/2010/sweepline/libs/sweepline/test/clipping_test.cpp 2010-12-20 09:30:57 EST (Mon, 20 Dec 2010)
@@ -1,11 +1,11 @@
-// Boost sweepline library clipping_test.cpp file
+// Boost sweepline library clipping_test.cpp file
// Copyright Andrii Sydorchuk 2010.
// 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.
+// See http://www.boost.org for updates, documentation, and revision history.
#include <stdlib.h>
#include <time.h>
@@ -95,7 +95,7 @@
intersections.clear();
point_2d<T> test_direction(static_cast<T>(i), static_cast<T>(j));
voronoi_helper<T>::find_intersections(test_origin, test_direction,
- voronoi_helper<T>::SEGMENT, test_rect, intersections);
+ voronoi_helper<T>::SEGMENT, test_rect, intersections);
if (abs(i) >= 2 || abs(j) >= 2)
BOOST_CHECK_EQUAL(static_cast<int>(intersections.size()), 1);
else
@@ -117,7 +117,7 @@
T y = static_cast<T>(j) / static_cast<T>(26);
point_2d<T> test_direction(x, y);
voronoi_helper<T>::find_intersections(test_origin, test_direction,
- voronoi_helper<T>::SEGMENT, test_rect, intersections);
+ voronoi_helper<T>::SEGMENT, test_rect, intersections);
BOOST_CHECK_EQUAL(static_cast<int>(intersections.size()), 0);
}
}
@@ -135,7 +135,7 @@
std::vector< point_2d<T> > intersections;
voronoi_helper<T>::find_intersections(test_origin, test_direction1,
- voronoi_helper<T>::RAY, test_rect, intersections);
+ voronoi_helper<T>::RAY, test_rect, intersections);
BOOST_CHECK_EQUAL(static_cast<int>(intersections.size()), 2);
BOOST_CHECK_EQUAL(intersections[0].x() == static_cast<T>(0) &&
intersections[0].y() == static_cast<T>(2), true);
@@ -144,7 +144,7 @@
intersections.clear();
voronoi_helper<T>::find_intersections(test_origin, test_direction2,
- voronoi_helper<T>::RAY, test_rect, intersections);
+ voronoi_helper<T>::RAY, test_rect, intersections);
BOOST_CHECK_EQUAL(static_cast<int>(intersections.size()), 2);
BOOST_CHECK_EQUAL(intersections[0].x() == static_cast<T>(0) &&
intersections[0].y() == static_cast<T>(1), true);
@@ -153,7 +153,7 @@
intersections.clear();
voronoi_helper<T>::find_intersections(test_origin, test_direction3,
- voronoi_helper<T>::RAY, test_rect, intersections);
+ voronoi_helper<T>::RAY, test_rect, intersections);
BOOST_CHECK_EQUAL(static_cast<int>(intersections.size()), 2);
BOOST_CHECK_EQUAL(intersections[0].x() == static_cast<T>(1) &&
intersections[0].y() == static_cast<T>(2), true);
@@ -162,14 +162,14 @@
intersections.clear();
voronoi_helper<T>::find_intersections(test_origin, test_direction4,
- voronoi_helper<T>::RAY, test_rect, intersections);
+ voronoi_helper<T>::RAY, test_rect, intersections);
BOOST_CHECK_EQUAL(static_cast<int>(intersections.size()), 1);
BOOST_CHECK_EQUAL(intersections[0].x() == static_cast<T>(0) &&
intersections[0].y() == static_cast<T>(-1), true);
intersections.clear();
voronoi_helper<T>::find_intersections(test_origin, test_direction5,
- voronoi_helper<T>::RAY, test_rect, intersections);
+ voronoi_helper<T>::RAY, test_rect, intersections);
BOOST_CHECK_EQUAL(static_cast<int>(intersections.size()), 1);
BOOST_CHECK_EQUAL(intersections[0].x() == static_cast<T>(4) &&
intersections[0].y() == static_cast<T>(2), true);
@@ -189,7 +189,7 @@
T y = static_cast<T>(j) / static_cast<T>(26);
point_2d<T> test_direction(x, y);
voronoi_helper<T>::find_intersections(test_origin, test_direction,
- voronoi_helper<T>::RAY, test_rect, intersections);
+ voronoi_helper<T>::RAY, test_rect, intersections);
if (i && j)
BOOST_CHECK_EQUAL(static_cast<int>(intersections.size()), 1);
}
@@ -208,7 +208,7 @@
std::vector< point_2d<T> > intersections;
voronoi_helper<T>::find_intersections(test_origin, test_direction1,
- voronoi_helper<T>::LINE, test_rect, intersections);
+ voronoi_helper<T>::LINE, test_rect, intersections);
BOOST_CHECK_EQUAL(static_cast<int>(intersections.size()), 2);
BOOST_CHECK_EQUAL(intersections[0].x() == static_cast<T>(3) &&
intersections[0].y() == static_cast<T>(-1), true);
@@ -217,7 +217,7 @@
intersections.clear();
voronoi_helper<T>::find_intersections(test_origin, test_direction2,
- voronoi_helper<T>::LINE, test_rect, intersections);
+ voronoi_helper<T>::LINE, test_rect, intersections);
BOOST_CHECK_EQUAL(static_cast<int>(intersections.size()), 2);
BOOST_CHECK_EQUAL(intersections[0].x() == static_cast<T>(1) &&
intersections[0].y() == static_cast<T>(-1), true);
@@ -226,7 +226,7 @@
intersections.clear();
voronoi_helper<T>::find_intersections(test_origin, test_direction3,
- voronoi_helper<T>::LINE, test_rect, intersections);
+ voronoi_helper<T>::LINE, test_rect, intersections);
BOOST_CHECK_EQUAL(static_cast<int>(intersections.size()), 2);
BOOST_CHECK_EQUAL(intersections[0].x() == static_cast<T>(4) &&
intersections[0].y() == static_cast<T>(0.5), true);
@@ -235,14 +235,14 @@
intersections.clear();
voronoi_helper<T>::find_intersections(test_origin, test_direction4,
- voronoi_helper<T>::LINE, test_rect, intersections);
+ voronoi_helper<T>::LINE, test_rect, intersections);
BOOST_CHECK_EQUAL(static_cast<int>(intersections.size()), 1);
BOOST_CHECK_EQUAL(intersections[0].x() == static_cast<T>(0) &&
intersections[0].y() == static_cast<T>(-1), true);
intersections.clear();
voronoi_helper<T>::find_intersections(test_origin, test_direction5,
- voronoi_helper<T>::LINE, test_rect, intersections);
+ voronoi_helper<T>::LINE, test_rect, intersections);
BOOST_CHECK_EQUAL(static_cast<int>(intersections.size()), 1);
BOOST_CHECK_EQUAL(intersections[0].x() == static_cast<T>(4) &&
intersections[0].y() == static_cast<T>(2), true);
@@ -262,7 +262,7 @@
T y = static_cast<T>(j) / static_cast<T>(26);
point_2d<T> test_direction(x, y);
voronoi_helper<T>::find_intersections(test_origin, test_direction,
- voronoi_helper<T>::LINE, test_rect, intersections);
+ voronoi_helper<T>::LINE, test_rect, intersections);
if (i && j)
BOOST_CHECK_EQUAL(static_cast<int>(intersections.size()), 2);
}
Modified: sandbox/SOC/2010/sweepline/libs/sweepline/test/event_queue_test.cpp
==============================================================================
--- sandbox/SOC/2010/sweepline/libs/sweepline/test/event_queue_test.cpp (original)
+++ sandbox/SOC/2010/sweepline/libs/sweepline/test/event_queue_test.cpp 2010-12-20 09:30:57 EST (Mon, 20 Dec 2010)
@@ -5,7 +5,7 @@
// (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.
+// See http://www.boost.org for updates, documentation, and revision history.
#include <cmath>
@@ -24,7 +24,7 @@
detail::circle_events_queue<T> event_q;
BOOST_CHECK_EQUAL(event_q.empty(), true);
- event_q.reset();
+ event_q.clear();
for (int i = 0; i < 10; i++) {
T x = static_cast<T>(-i);
T y = static_cast<T>(10-i);
Modified: sandbox/SOC/2010/sweepline/libs/sweepline/test/event_types_test.cpp
==============================================================================
--- sandbox/SOC/2010/sweepline/libs/sweepline/test/event_types_test.cpp (original)
+++ sandbox/SOC/2010/sweepline/libs/sweepline/test/event_types_test.cpp 2010-12-20 09:30:57 EST (Mon, 20 Dec 2010)
@@ -5,7 +5,7 @@
// (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.
+// See http://www.boost.org for updates, documentation, and revision history.
#include "test_type_list.hpp"
#include "boost/sweepline/voronoi_sweepline.hpp"
@@ -49,7 +49,7 @@
BOOST_CHECK_EQUAL(site1.x(), static_cast<T>(1));
BOOST_CHECK_EQUAL(site1.y(), static_cast<T>(2));
- BOOST_CHECK_EQUAL(site1.get_site_index(), 0);
+ BOOST_CHECK_EQUAL(site1.index(), 0);
site2 = make_site_event<T>(static_cast<T>(0), static_cast<T>(2), 1);
bool arr1[] = { false, true, false, true, false, true };
@@ -73,8 +73,8 @@
point_2d<T> point5 = make_point_2d<T>(static_cast<T>(1), static_cast<T>(1));
site_event<T> site1 = make_site_event<T>(point1, point2, 0);
- BOOST_CHECK_EQUAL(point1 == site1.get_point1(), true);
- BOOST_CHECK_EQUAL(point2 == site1.get_point0(), true);
+ BOOST_CHECK_EQUAL(point1 == site1.point1(), true);
+ BOOST_CHECK_EQUAL(point2 == site1.point0(), true);
BOOST_CHECK_EQUAL(site1.is_segment(), true);
BOOST_CHECK_EQUAL(site1.is_vertical(), true);
@@ -119,8 +119,8 @@
point_2d<T> point4 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(10));
site_event<T> site1 = make_site_event<T>(point1, point2, 0);
- BOOST_CHECK_EQUAL(point1 == site1.get_point1(), true);
- BOOST_CHECK_EQUAL(point2 == site1.get_point0(), true);
+ BOOST_CHECK_EQUAL(point1 == site1.point1(), true);
+ BOOST_CHECK_EQUAL(point2 == site1.point0(), true);
BOOST_CHECK_EQUAL(site1.is_segment(), true);
BOOST_CHECK_EQUAL(site1.is_vertical(), false);
Modified: sandbox/SOC/2010/sweepline/libs/sweepline/test/node_comparer_test.cpp
==============================================================================
--- sandbox/SOC/2010/sweepline/libs/sweepline/test/node_comparer_test.cpp (original)
+++ sandbox/SOC/2010/sweepline/libs/sweepline/test/node_comparer_test.cpp 2010-12-20 09:30:57 EST (Mon, 20 Dec 2010)
@@ -5,7 +5,7 @@
// (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.
+// See http://www.boost.org for updates, documentation, and revision history.
#include "test_type_list.hpp"
#include "boost/sweepline/voronoi_sweepline.hpp"
@@ -248,11 +248,3 @@
BOOST_CHECK_EQUAL(node_comparer_test(new_node3, initial_node), false);
BOOST_CHECK_EQUAL(node_comparer_test(new_node4, initial_node), false);
}
-
-BOOST_AUTO_TEST_CASE_TEMPLATE(node_comparer_test_ps1, T, test_types) {
- // TODO(asydorchuk): add more tests there.
-}
-
-BOOST_AUTO_TEST_CASE_TEMPLATE(node_comparer_test_ss1, T, test_types) {
- // TODO(asydorchuk): add more tests there.
-}
Modified: sandbox/SOC/2010/sweepline/libs/sweepline/test/output_verification.hpp
==============================================================================
--- sandbox/SOC/2010/sweepline/libs/sweepline/test/output_verification.hpp (original)
+++ sandbox/SOC/2010/sweepline/libs/sweepline/test/output_verification.hpp 2010-12-20 09:30:57 EST (Mon, 20 Dec 2010)
@@ -5,7 +5,7 @@
// (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.
+// See http://www.boost.org for updates, documentation, and revision history.
#ifndef VORONOI_OUTPUT_VERIFICATION
#define VORONOI_OUTPUT_VERIFICATION
@@ -13,32 +13,28 @@
#include <map>
#include <vector>
-enum kOrientation {
- RIGHT_ORIENTATION = -1,
- COLLINEAR = 0,
- LEFT_ORIENTATION = 1,
-};
-
template <typename Point2D>
-kOrientation get_orientation(const Point2D &point1, const Point2D &point2, const Point2D &point3) {
+detail::kOrientation get_orientation(const Point2D &point1,
+ const Point2D &point2,
+ const Point2D &point3) {
typename Point2D::coordinate_type a = (point2.x() - point1.x()) * (point3.y() - point2.y());
typename Point2D::coordinate_type b = (point2.y() - point1.y()) * (point3.x() - point2.x());
if (a == b)
- return COLLINEAR;
- return (a < b) ? RIGHT_ORIENTATION : LEFT_ORIENTATION;
+ return detail::COLLINEAR;
+ return (a < b) ? detail::RIGHT_ORIENTATION : detail::LEFT_ORIENTATION;
}
template <typename Output>
bool verify_half_edge_orientation(const Output &output) {
- typedef typename Output::Point2D Point2D;
+ typedef typename Output::point_2d_type point_2d_type;
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->vertex0() != NULL && edge_it->vertex1() != NULL) {
- const Point2D &site_point = edge_it->cell()->get_point0();
- const Point2D &start_point = edge_it->vertex0()->vertex();
- const Point2D &end_point = edge_it->vertex1()->vertex();
- if (get_orientation(start_point, end_point, site_point) != LEFT_ORIENTATION)
+ if (edge_it->is_bounded()) {
+ const point_2d_type &site_point = edge_it->cell()->point0();
+ const point_2d_type &start_point = edge_it->vertex0()->vertex();
+ const point_2d_type &end_point = edge_it->vertex1()->vertex();
+ if (get_orientation(start_point, end_point, site_point) != detail::LEFT_ORIENTATION)
return false;
}
}
@@ -47,7 +43,7 @@
template <typename Output>
bool verify_cell_convexity(const Output &output) {
- typedef typename Output::Point2D Point2D;
+ typedef typename Output::point_2d_type point_2d_type;
typename Output::voronoi_cell_const_iterator_type cell_it;
for (cell_it = output.cell_records().begin();
cell_it != output.cell_records().end(); cell_it++) {
@@ -62,10 +58,10 @@
edge->vertex1() != NULL &&
edge->vertex1() == edge->next()->vertex0() &&
edge->next()->vertex1() != NULL) {
- const Point2D &vertex1 = edge->vertex0()->vertex();
- const Point2D &vertex2 = edge->vertex1()->vertex();
- const Point2D &vertex3 = edge->next()->vertex1()->vertex();
- if (get_orientation(vertex1, vertex2, vertex3) != LEFT_ORIENTATION)
+ const point_2d_type &vertex1 = edge->vertex0()->vertex();
+ const point_2d_type &vertex2 = edge->vertex1()->vertex();
+ const point_2d_type &vertex3 = edge->next()->vertex1()->vertex();
+ if (get_orientation(vertex1, vertex2, vertex3) != detail::LEFT_ORIENTATION)
return false;
}
edge = edge->next();
@@ -76,7 +72,7 @@
template <typename Output>
bool verify_incident_edges_ccw_order(const Output &output) {
- typedef typename Output::Point2D Point2D;
+ typedef typename Output::point_2d_type point_2d_type;
typedef typename Output::voronoi_edge_type voronoi_edge_type;
typename Output::voronoi_vertex_const_iterator_type vertex_it;
for (vertex_it = output.vertex_records().begin();
@@ -87,11 +83,11 @@
do {
const voronoi_edge_type *edge_next1 = edge->rot_next();
const voronoi_edge_type *edge_next2 = edge_next1->rot_next();
- const Point2D &site1 = edge->cell()->get_point0();
- const Point2D &site2 = edge_next1->cell()->get_point0();
- const Point2D &site3 = edge_next2->cell()->get_point0();
+ const point_2d_type &site1 = edge->cell()->point0();
+ const point_2d_type &site2 = edge_next1->cell()->point0();
+ const point_2d_type &site3 = edge_next2->cell()->point0();
- if (get_orientation(site1, site2, site3) != LEFT_ORIENTATION)
+ if (get_orientation(site1, site2, site3) != detail::LEFT_ORIENTATION)
return false;
edge = edge->rot_next();
@@ -101,22 +97,20 @@
}
template <typename Output>
-bool verfiy_no_half_edge_intersections(const Output &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::Point2D Point2D;
- std::map< Point2D, std::vector<Point2D> > edge_map;
- typedef typename std::map< Point2D, std::vector<Point2D> >::const_iterator
+ typedef typename Output::point_2d_type point_2d_type;
+ std::map< point_2d_type, std::vector<point_2d_type> > edge_map;
+ typedef typename std::map< point_2d_type, std::vector<point_2d_type> >::const_iterator
edge_map_iterator;
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->cell()->contains_segment() && edge_it->twin()->cell()->contains_segment())
- continue;
- if (edge_it->vertex0() != NULL && edge_it->vertex1() != NULL) {
- const Point2D &vertex0 = edge_it->vertex0()->vertex();
- const Point2D &vertex1 = edge_it->vertex1()->vertex();
+ if (edge_it->is_bounded()) {
+ const point_2d_type &vertex0 = edge_it->vertex0()->vertex();
+ const point_2d_type &vertex1 = edge_it->vertex1()->vertex();
if (vertex0 < vertex1)
edge_map[vertex0].push_back(vertex1);
}
@@ -127,13 +121,13 @@
// left to right checking for intersections between all pairs of edges
// that overlap in the x dimension.
// Complexity. Approximately N*sqrt(N). Worst case N^2.
- typedef typename std::vector<Point2D>::size_type size_type;
+ typedef typename std::vector<point_2d_type>::size_type size_type;
edge_map_iterator edge_map_it1, edge_map_it2, edge_map_it_bound;
for (edge_map_it1 = edge_map.begin();
edge_map_it1 != edge_map.end(); edge_map_it1++) {
- const Point2D &point1 = edge_map_it1->first;
+ const point_2d_type &point1 = edge_map_it1->first;
for (size_type i = 0; i < edge_map_it1->second.size(); i++) {
- const Point2D &point2 = edge_map_it1->second[i];
+ const point_2d_type &point2 = edge_map_it1->second[i];
typename Output::coordinate_type min_y1 = std::min(point1.y(), point2.y());
typename Output::coordinate_type max_y1 = std::max(point1.y(), point2.y());
@@ -143,9 +137,9 @@
edge_map_it2 = edge_map_it1;
edge_map_it2++;
for (; edge_map_it2 != edge_map_it_bound; edge_map_it2++) {
- const Point2D &point3 = edge_map_it2->first;
+ const point_2d_type &point3 = edge_map_it2->first;
for (size_type j = 0; j < edge_map_it2->second.size(); j++) {
- const Point2D &point4 = edge_map_it2->second[j];
+ const point_2d_type &point4 = edge_map_it2->second[j];
typename Output::coordinate_type min_y2 = std::min(point3.y(), point4.y());
typename Output::coordinate_type max_y2 = std::max(point3.y(), point4.y());
@@ -156,9 +150,9 @@
// Intersection check.
if (get_orientation(point1, point2, point3) *
- get_orientation(point1, point2, point4) == RIGHT_ORIENTATION &&
+ get_orientation(point1, point2, point4) == detail::RIGHT_ORIENTATION &&
get_orientation(point3, point4, point1) *
- get_orientation(point3, point4, point2) == RIGHT_ORIENTATION)
+ get_orientation(point3, point4, point2) == detail::RIGHT_ORIENTATION)
return false;
}
}
@@ -167,36 +161,32 @@
return true;
}
-template<typename P>
-struct segm_comparison {
- bool operator()(const std::pair<P, P>& segm1, const std::pair<P, P>& segm2) const {
- return segm1.first.x() < segm2.first.x();
- }
-};
-
-template<typename P>
-bool intersect(const std::pair<P, P> &segm1, const std::pair<P, P> &segm2) {
- if ((std::max)(segm1.first.y(), segm1.second.y()) <
- (std::min)(segm2.first.y(), segm2.second.y()) ||
- (std::max)(segm2.first.y(), segm2.second.y()) <
- (std::min)(segm1.first.y(), segm1.second.y()))
+template <typename P>
+bool do_intersect(const P &p1, const P &p2, const P &p3, const P &p4) {
+ if ((std::max)(p1.y(), p2.y()) < (std::min)(p3.y(), p4.y()) ||
+ (std::max)(p3.y(), p4.y()) < (std::min)(p1.y(), p2.y()))
+ return false;
+ if (p1 == p3)
+ return detail::orientation_test(p1, p2, p4) == detail::COLLINEAR;
+ if (p2 == p4)
+ return detail::orientation_test(p1, p2, p3) == detail::COLLINEAR;
+ if (p2 == p3)
return false;
- if (segm1.first == segm2.first)
- return detail::orientation_test(segm1.first, segm1.second, segm2.second)
- == detail::COLINEAR;
- if (segm1.second == segm2.second)
- return detail::orientation_test(segm1.first, segm1.second, segm2.first)
- == detail::COLINEAR;
- if ((get_orientation(segm1.first, segm1.second, segm2.first) *
- get_orientation(segm1.first, segm1.second, segm2.second) != LEFT_ORIENTATION) &&
- (get_orientation(segm2.first, segm2.second, segm1.first) *
- get_orientation(segm2.first, segm2.second, segm1.second) != LEFT_ORIENTATION))
+ if ((get_orientation(p1, p2, p3) * get_orientation(p1, p2, p4) != detail::LEFT_ORIENTATION) &&
+ (get_orientation(p3, p4, p1) * get_orientation(p3, p4, p2) != detail::LEFT_ORIENTATION))
return true;
return false;
}
-template<typename P>
-void remove_intersections(std::vector< std::pair<P, P> >& segm_vec) {
+template <typename P>
+struct segm_comparison {
+ bool operator()(const std::pair<P, P> &segm1, const std::pair<P, P> &segm2) const {
+ return segm1.first.x() < segm2.first.x();
+ }
+};
+
+template <typename P>
+void remove_intersections(std::vector< std::pair<P, P> > &segm_vec) {
for (int i = 0; i < static_cast<int>(segm_vec.size()); i++) {
if (segm_vec[i].first > segm_vec[i].second) {
(std::swap(segm_vec[i].first, segm_vec[i].second));
@@ -210,7 +200,7 @@
l_bound = upper_bound(segm_vec.begin(), segm_vec.end(), bound, segm_comparison<P>());
bool add = true;
for (b_it = it + 1; b_it != l_bound; b_it++) {
- if (intersect(*it, *b_it)) {
+ if (do_intersect(it->first, it->second, b_it->first, b_it->second)) {
add = false;
break;
}
@@ -241,7 +231,7 @@
if (mask & INCIDENT_EDGES_CCW_ORDER)
result &= verify_incident_edges_ccw_order(output);
if (mask & NO_HALF_EDGE_INTERSECTIONS)
- result &= verfiy_no_half_edge_intersections(output);
+ result &= verfiy_no_line_edge_intersections(output);
return result;
}
Modified: sandbox/SOC/2010/sweepline/libs/sweepline/test/predicates_test.cpp
==============================================================================
--- sandbox/SOC/2010/sweepline/libs/sweepline/test/predicates_test.cpp (original)
+++ sandbox/SOC/2010/sweepline/libs/sweepline/test/predicates_test.cpp 2010-12-20 09:30:57 EST (Mon, 20 Dec 2010)
@@ -1,11 +1,11 @@
-// Boost sweepline library predicates_test.cpp file
+// Boost sweepline library predicates_test.cpp file
// Copyright Andrii Sydorchuk 2010.
// 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.
+// See http://www.boost.org for updates, documentation, and revision history.
#include "test_type_list.hpp"
#include "boost/sweepline/voronoi_sweepline.hpp"
@@ -17,8 +17,8 @@
#define CHECK_ORIENTATION_EQUAL(p1, p2, p3, exp) \
BOOST_CHECK_EQUAL(detail::orientation_test(p1, p2, p3) == exp, true)
-#define CHECK_FAST_LESS_PREDICATE_PP(lp, rp, np, exp) \
- BOOST_CHECK_EQUAL(detail::fast_less_predicate(lp, rp, np) == exp, true)
+#define CHECK_EPS_LESS_PREDICATE_PP(lp, rp, np, exp) \
+ BOOST_CHECK_EQUAL(detail::eps_less_predicate(lp, rp, np) == exp, true)
#define CHECK_LESS_PREDICATE_PP(lp, rp, np, exp) \
BOOST_CHECK_EQUAL(detail::less_predicate(lp, rp, np) == exp, true)
@@ -47,12 +47,12 @@
point_2d<T> point4 = make_point_2d<T>(min_int, max_int);
point_2d<T> point5 = make_point_2d<T>(max_int - 1, max_int);
- CHECK_ORIENTATION_EQUAL(point1, point2, point3, detail::COLINEAR);
- CHECK_ORIENTATION_EQUAL(point1, point3, point2, detail::COLINEAR);
- CHECK_ORIENTATION_EQUAL(point2, point3, point1, detail::COLINEAR);
- CHECK_ORIENTATION_EQUAL(point2, point1, point3, detail::COLINEAR);
- CHECK_ORIENTATION_EQUAL(point3, point1, point2, detail::COLINEAR);
- CHECK_ORIENTATION_EQUAL(point3, point2, point1, detail::COLINEAR);
+ CHECK_ORIENTATION_EQUAL(point1, point2, point3, detail::COLLINEAR);
+ CHECK_ORIENTATION_EQUAL(point1, point3, point2, detail::COLLINEAR);
+ CHECK_ORIENTATION_EQUAL(point2, point3, point1, detail::COLLINEAR);
+ CHECK_ORIENTATION_EQUAL(point2, point1, point3, detail::COLLINEAR);
+ CHECK_ORIENTATION_EQUAL(point3, point1, point2, detail::COLLINEAR);
+ CHECK_ORIENTATION_EQUAL(point3, point2, point1, detail::COLLINEAR);
CHECK_ORIENTATION_EQUAL(point1, point4, point3, detail::RIGHT_ORIENTATION);
CHECK_ORIENTATION_EQUAL(point1, point3, point4, detail::LEFT_ORIENTATION);
@@ -76,16 +76,16 @@
point_2d<T> point3 = make_point_2d<T>(static_cast<T>(-2), static_cast<T>(1));
point_2d<T> site1 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(5));
- CHECK_FAST_LESS_PREDICATE_PP(point1, point2, site1, detail::UNDEFINED);
- CHECK_FAST_LESS_PREDICATE_PP(point3, point1, site1, detail::UNDEFINED);
+ CHECK_EPS_LESS_PREDICATE_PP(point1, point2, site1, detail::UNDEFINED);
+ CHECK_EPS_LESS_PREDICATE_PP(point3, point1, site1, detail::UNDEFINED);
point_2d<T> site2 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(4));
- CHECK_FAST_LESS_PREDICATE_PP(point1, point2, site2, detail::MORE);
- CHECK_FAST_LESS_PREDICATE_PP(point3, point1, site2, detail::MORE);
+ CHECK_EPS_LESS_PREDICATE_PP(point1, point2, site2, detail::MORE);
+ CHECK_EPS_LESS_PREDICATE_PP(point3, point1, site2, detail::MORE);
point_2d<T> site3 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(6));
- CHECK_FAST_LESS_PREDICATE_PP(point1, point2, site3, detail::LESS);
- CHECK_FAST_LESS_PREDICATE_PP(point3, point1, site3, detail::LESS);
+ CHECK_EPS_LESS_PREDICATE_PP(point1, point2, site3, detail::LESS);
+ CHECK_EPS_LESS_PREDICATE_PP(point3, point1, site3, detail::LESS);
}
// Test main point-point predicate.
@@ -128,13 +128,13 @@
point_2d<T> segm_start = make_point_2d<T>(static_cast<T>(-5), static_cast<T>(5));
point_2d<T> segm_end = make_point_2d<T>(static_cast<T>(2), static_cast<T>(-2));
typename detail::site_event<T> segm_site(segm_start, segm_end, 0);
-
+
point_2d<T> site_p1 = make_point_2d<T>(static_cast<T>(-2), static_cast<T>(4));
point_2d<T> new_p1 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(-1));
- segm_site.set_inverse();
+ segm_site.inverse();
CHECK_FAST_LESS_PREDICATE_PS(site_p1, segm_site, new_p1, false, detail::MORE);
CHECK_FAST_LESS_PREDICATE_PS(site_p1, segm_site, new_p1, true, detail::MORE);
-
+
point_2d<T> new_p2 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(1));
CHECK_FAST_LESS_PREDICATE_PS(site_p1, segm_site, new_p2, false, detail::MORE);
CHECK_FAST_LESS_PREDICATE_PS(site_p1, segm_site, new_p2, true, detail::UNDEFINED);
@@ -153,7 +153,7 @@
point_2d<T> segm_start = make_point_2d<T>(static_cast<T>(-5), static_cast<T>(5));
point_2d<T> segm_end = make_point_2d<T>(static_cast<T>(2), static_cast<T>(-2));
typename detail::site_event<T> segm_site(segm_start, segm_end, 0);
-
+
point_2d<T> site_p1 = make_point_2d<T>(static_cast<T>(-2), static_cast<T>(-4));
point_2d<T> site_p2 = make_point_2d<T>(static_cast<int>(-4), static_cast<int>(1));
@@ -162,7 +162,7 @@
CHECK_FAST_LESS_PREDICATE_PS(site_p1, segm_site, new_p1, true, detail::LESS);
CHECK_FAST_LESS_PREDICATE_PS(site_p2, segm_site, new_p1, false, detail::LESS);
CHECK_FAST_LESS_PREDICATE_PS(site_p2, segm_site, new_p1, true, detail::LESS);
-
+
point_2d<T> new_p2 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(-2));
CHECK_FAST_LESS_PREDICATE_PS(site_p1, segm_site, new_p2, false, detail::UNDEFINED);
CHECK_FAST_LESS_PREDICATE_PS(site_p1, segm_site, new_p2, true, detail::LESS);
@@ -188,8 +188,8 @@
point_2d<T> segm_end = make_point_2d<T>(static_cast<T>(2), static_cast<T>(-2));
typename detail::site_event<T> segm_site1(segm_start, segm_end, 0);
typename detail::site_event<T> segm_site2(segm_start, segm_end, 0);
- segm_site2.set_inverse();
-
+ segm_site2.inverse();
+
point_2d<T> site_p1 = make_point_2d<T>(static_cast<T>(-2), static_cast<T>(4));
point_2d<T> site_p2 = make_point_2d<T>(static_cast<T>(-2), static_cast<T>(-4));
point_2d<T> site_p3 = make_point_2d<T>(static_cast<int>(-4), static_cast<int>(1));
@@ -233,7 +233,7 @@
point_2d<T> segm_end1 = make_point_2d<T>(static_cast<T>(2), static_cast<T>(7));
typename detail::site_event<T> segm_site1_1(segm_start1, segm_end1, 0);
typename detail::site_event<T> segm_site1_2(segm_start1, segm_end1, 0);
- segm_site1_2.set_inverse();
+ segm_site1_2.inverse();
// New sites.
point_2d<T> new_site1 = make_point_2d<T>(static_cast<T>(2), static_cast<T>(7));
@@ -252,7 +252,7 @@
point_2d<T> segm_end1 = make_point_2d<T>(static_cast<T>(2), static_cast<T>(-2));
site_event_type segm_site1_1(segm_start1, segm_end1, 0);
site_event_type segm_site1_2(segm_start1, segm_end1, 0);
- segm_site1_2.set_inverse();
+ segm_site1_2.inverse();
// New sites.
point_2d<T> new_site1 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(2));
@@ -277,7 +277,7 @@
CHECK_LESS_PREDICATE_SS(segm_site1_2, segm_site3, new_site2, true);
CHECK_LESS_PREDICATE_SS(segm_site1_2, segm_site3, new_site3, true);
CHECK_LESS_PREDICATE_SS(segm_site1_2, segm_site3, new_site4, true);
- segm_site3.set_inverse();
+ segm_site3.inverse();
CHECK_LESS_PREDICATE_SS(segm_site3, segm_site1_2, new_site1, false);
CHECK_LESS_PREDICATE_SS(segm_site3, segm_site1_2, new_site2, false);
CHECK_LESS_PREDICATE_SS(segm_site3, segm_site1_2, new_site3, false);
@@ -291,7 +291,7 @@
point_2d<T> segm_start2 = make_point_2d<T>(static_cast<T>(-5), static_cast<T>(5));
point_2d<T> segm_end = make_point_2d<T>(static_cast<T>(-2), static_cast<T>(2));
site_event_type segm_site1(segm_start1, segm_end, 0);
- segm_site1.set_inverse();
+ segm_site1.inverse();
site_event_type segm_site2(segm_start2, segm_end, 1);
point_2d<T> point(-4, 2);
CHECK_LESS_PREDICATE_SS(segm_site1, segm_site2, point, false);
Modified: sandbox/SOC/2010/sweepline/libs/sweepline/test/sweepline_test.cpp
==============================================================================
--- sandbox/SOC/2010/sweepline/libs/sweepline/test/sweepline_test.cpp (original)
+++ sandbox/SOC/2010/sweepline/libs/sweepline/test/sweepline_test.cpp 2010-12-20 09:30:57 EST (Mon, 20 Dec 2010)
@@ -1,11 +1,11 @@
-// Boost sweepline library builder_test.cpp file
+// Boost sweepline library builder_test.cpp file
// Copyright Andrii Sydorchuk 2010.
// 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.
+// See http://www.boost.org for updates, documentation, and revision history.
#include <stdlib.h>
#include <time.h>
@@ -187,7 +187,7 @@
points.push_back(point1);
points.push_back(point2);
points.push_back(point3);
-
+
voronoi_output<coordinate_type> test_output;
build_voronoi(points, test_output);
VERIFY_VORONOI_OUTPUT(test_output, COMPLETE_VERIFICATION);
@@ -208,18 +208,18 @@
const voronoi_edge_type *edge1_1 = it->incident_edge();
const voronoi_edge_type *edge1_2 = edge1_1->twin();
- CHECK_EQUAL_POINTS(edge1_1->cell()->get_point0(), point3);
- CHECK_EQUAL_POINTS(edge1_2->cell()->get_point0(), point1);
+ CHECK_EQUAL_POINTS(edge1_1->cell()->point0(), point3);
+ CHECK_EQUAL_POINTS(edge1_2->cell()->point0(), point1);
const voronoi_edge_type *edge2_1 = edge1_1->rot_prev();
const voronoi_edge_type *edge2_2 = edge2_1->twin();
- CHECK_EQUAL_POINTS(edge2_1->cell()->get_point0(), point1);
- CHECK_EQUAL_POINTS(edge2_2->cell()->get_point0(), point2);
+ CHECK_EQUAL_POINTS(edge2_1->cell()->point0(), point1);
+ CHECK_EQUAL_POINTS(edge2_2->cell()->point0(), point2);
const voronoi_edge_type *edge3_1 = edge2_1->rot_prev();
const voronoi_edge_type *edge3_2 = edge3_1->twin();
- CHECK_EQUAL_POINTS(edge3_1->cell()->get_point0(), point2);
- CHECK_EQUAL_POINTS(edge3_2->cell()->get_point0(), point3);
+ CHECK_EQUAL_POINTS(edge3_1->cell()->point0(), point2);
+ CHECK_EQUAL_POINTS(edge3_2->cell()->point0(), point3);
BOOST_CHECK_EQUAL(edge1_2->twin() == edge1_1, true);
BOOST_CHECK_EQUAL(edge2_2->twin() == edge2_1, true);
@@ -256,7 +256,7 @@
points.push_back(point1);
points.push_back(point2);
points.push_back(point3);
-
+
voronoi_output<coordinate_type> test_output;
build_voronoi(points, test_output);
VERIFY_VORONOI_OUTPUT(test_output, COMPLETE_VERIFICATION);
@@ -271,18 +271,18 @@
const voronoi_edge_type *edge1_1 = it->incident_edge();
const voronoi_edge_type *edge1_2 = edge1_1->twin();
- CHECK_EQUAL_POINTS(edge1_1->cell()->get_point0(), point2);
- CHECK_EQUAL_POINTS(edge1_2->cell()->get_point0(), point1);
+ CHECK_EQUAL_POINTS(edge1_1->cell()->point0(), point2);
+ CHECK_EQUAL_POINTS(edge1_2->cell()->point0(), point1);
const voronoi_edge_type *edge2_1 = edge1_1->rot_prev();
- const voronoi_edge_type *edge2_2 = edge2_1->twin();
- CHECK_EQUAL_POINTS(edge2_1->cell()->get_point0(), point1);
- CHECK_EQUAL_POINTS(edge2_2->cell()->get_point0(), point3);
+ const voronoi_edge_type *edge2_2 = edge2_1->twin();
+ CHECK_EQUAL_POINTS(edge2_1->cell()->point0(), point1);
+ CHECK_EQUAL_POINTS(edge2_2->cell()->point0(), point3);
const voronoi_edge_type *edge3_1 = edge2_1->rot_prev();
const voronoi_edge_type *edge3_2 = edge3_1->twin();
- CHECK_EQUAL_POINTS(edge3_1->cell()->get_point0(), point3);
- CHECK_EQUAL_POINTS(edge3_2->cell()->get_point0(), point2);
+ CHECK_EQUAL_POINTS(edge3_1->cell()->point0(), point3);
+ CHECK_EQUAL_POINTS(edge3_2->cell()->point0(), point2);
BOOST_CHECK_EQUAL(edge1_2->twin() == edge1_1, true);
BOOST_CHECK_EQUAL(edge2_2->twin() == edge2_1, true);
@@ -317,7 +317,7 @@
static_cast<coordinate_type>(0)));
points.push_back(make_point_2d<coordinate_type>(static_cast<coordinate_type>(1),
static_cast<coordinate_type>(1)));
-
+
voronoi_output<coordinate_type> test_output;
build_voronoi(points, test_output);
VERIFY_VORONOI_OUTPUT(test_output, COMPLETE_VERIFICATION);
@@ -335,23 +335,23 @@
// Check voronoi edges.
const voronoi_edge_type *edge1_1 = it->incident_edge();
const voronoi_edge_type *edge1_2 = edge1_1->twin();
- CHECK_EQUAL_POINTS(edge1_1->cell()->get_point0(), points[1]);
- CHECK_EQUAL_POINTS(edge1_2->cell()->get_point0(), points[3]);
+ CHECK_EQUAL_POINTS(edge1_1->cell()->point0(), points[1]);
+ CHECK_EQUAL_POINTS(edge1_2->cell()->point0(), points[3]);
const voronoi_edge_type *edge2_1 = edge1_1->rot_prev();
- const voronoi_edge_type *edge2_2 = edge2_1->twin();
- CHECK_EQUAL_POINTS(edge2_1->cell()->get_point0(), points[3]);
- CHECK_EQUAL_POINTS(edge2_2->cell()->get_point0(), points[2]);
+ const voronoi_edge_type *edge2_2 = edge2_1->twin();
+ CHECK_EQUAL_POINTS(edge2_1->cell()->point0(), points[3]);
+ CHECK_EQUAL_POINTS(edge2_2->cell()->point0(), points[2]);
const voronoi_edge_type *edge3_1 = edge2_1->rot_prev();
const voronoi_edge_type *edge3_2 = edge3_1->twin();
- CHECK_EQUAL_POINTS(edge3_1->cell()->get_point0(), points[2]);
- CHECK_EQUAL_POINTS(edge3_2->cell()->get_point0(), points[0]);
+ CHECK_EQUAL_POINTS(edge3_1->cell()->point0(), points[2]);
+ CHECK_EQUAL_POINTS(edge3_2->cell()->point0(), points[0]);
const voronoi_edge_type *edge4_1 = edge3_1->rot_prev();
const voronoi_edge_type *edge4_2 = edge4_1->twin();
- CHECK_EQUAL_POINTS(edge4_1->cell()->get_point0(), points[0]);
- CHECK_EQUAL_POINTS(edge4_2->cell()->get_point0(), points[1]);
+ CHECK_EQUAL_POINTS(edge4_1->cell()->point0(), points[0]);
+ CHECK_EQUAL_POINTS(edge4_2->cell()->point0(), points[1]);
BOOST_CHECK_EQUAL(edge1_2->twin() == edge1_1, true);
BOOST_CHECK_EQUAL(edge2_2->twin() == edge2_1, true);
@@ -574,7 +574,7 @@
BOOST_CHECK_EQUAL(test_output.num_vertex_records(), 6);
BOOST_CHECK_EQUAL(test_output.num_edge_records(), 12);
VERIFY_VORONOI_OUTPUT(test_output, NO_HALF_EDGE_INTERSECTIONS);
-
+
}
BOOST_AUTO_TEST_CASE_TEMPLATE(segment_site_test8, T, test_types) {
@@ -638,7 +638,40 @@
#endif
#ifdef NDEBUG
-BOOST_AUTO_TEST_CASE_TEMPLATE(segment_random_test, T, test_types) {
+BOOST_AUTO_TEST_CASE_TEMPLATE(segment_random_test1, T, test_types) {
+ srand(static_cast<unsigned int>(time(NULL)));
+ voronoi_output<T> test_output;
+ std::vector< point_2d<T> > point_vec;
+ std::vector< std::pair< point_2d<T>, point_2d<T> > > segm_vec;
+ int num_runs = 100000;
+ int num_segments = 3;
+ point_vec.push_back(make_point_2d<T>(-100, -100));
+ point_vec.push_back(make_point_2d<T>(-100, 100));
+ point_vec.push_back(make_point_2d<T>(100, -100));
+ point_vec.push_back(make_point_2d<T>(100, 100));
+ for (int i = 0; i < num_runs; i++) {
+ for (int j = 0; j < num_segments; j++) {
+ T x1 = 0, y1 = 0, x2 = 0, y2 = 0;
+ while (x1 == x2 && y1 == y2) {
+ x1 = (rand() % 100) - 50;
+ y1 = (rand() % 100) - 50;
+ x2 = (rand() % 100) - 50;
+ y2 = (rand() % 100) - 50;
+ }
+ point_2d<T> point1(x1, y1);
+ point_2d<T> point2(x2, y2);
+ segm_vec.push_back(std::make_pair(point1, point2));
+ }
+ remove_intersections(segm_vec);
+ build_voronoi(point_vec, segm_vec, test_output);
+ VERIFY_VORONOI_OUTPUT(test_output, NO_HALF_EDGE_INTERSECTIONS);
+ segm_vec.clear();
+ }
+}
+#endif
+
+#ifdef NDEBUG
+BOOST_AUTO_TEST_CASE_TEMPLATE(segment_random_test2, T, test_types) {
srand(static_cast<unsigned int>(time(NULL)));
voronoi_output<T> test_output_small, test_output_large;
std::vector< std::pair< point_2d<T>, point_2d<T> > > segm_vec;
@@ -648,12 +681,12 @@
int mod_koef2[] = {100, 200, 300, 400};
int max_value[] = {100, 600, 5150, 50200};
int array_length = sizeof(num_segments) / sizeof(int);
- for (int k = 0; k < array_length; k++) {
+ for (int k = 3; k < array_length; k++) {
int koef = std::numeric_limits<int>::max() / max_value[k];
for (int i = 0; i < num_runs[k]; i++) {
for (int j = 0; j < num_segments[k]; j++) {
- T x1 = (rand() % mod_koef1[k]) - mod_koef1[k] / 2;
- T y1 = (rand() % mod_koef1[k]) - mod_koef1[k] / 2;
+ T x1 = (rand() % (mod_koef1[k] / 100)) - mod_koef1[k] / 2;
+ T y1 = (rand() % (mod_koef1[k] / 100)) - mod_koef1[k] / 2;
T dx = 0, dy = 0;
while (dx == 0 && dy == 0) {
dx = (rand() % mod_koef2[k]) - mod_koef2[k] / 2;
Modified: sandbox/SOC/2010/sweepline/libs/sweepline/test/test_type_list.hpp
==============================================================================
--- sandbox/SOC/2010/sweepline/libs/sweepline/test/test_type_list.hpp (original)
+++ sandbox/SOC/2010/sweepline/libs/sweepline/test/test_type_list.hpp 2010-12-20 09:30:57 EST (Mon, 20 Dec 2010)
@@ -5,7 +5,7 @@
// (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.
+// See http://www.boost.org for updates, documentation, and revision history.
#ifndef BOOST_SWEEPLINE_TEST_TYPE_LIST
#define BOOST_SWEEPLINE_TEST_TYPE_LIST
Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk