|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r65651 - in sandbox/SOC/2010/sweepline: boost/sweepline boost/sweepline/detail libs/sweepline/test libs/sweepline/test/voronoi_segment
From: sydorchuk.andriy_at_[hidden]
Date: 2010-09-28 15:28:13
Author: asydorchuk
Date: 2010-09-28 15:28:09 EDT (Tue, 28 Sep 2010)
New Revision: 65651
URL: http://svn.boost.org/trac/boost/changeset/65651
Log:
Finished all predicates ((point, segment, segment) and (segment, segment, segment) circle creation).
Added circle event creation tests.
Changed segment sites logic (two oriented vectors).
Fixed previous predicates according to the new logic.
Basic tests of voronoi of segments.
Added:
sandbox/SOC/2010/sweepline/libs/sweepline/test/voronoi_segment/segment_circle_event_creation_test.cpp (contents, props changed)
Text files modified:
sandbox/SOC/2010/sweepline/boost/sweepline/detail/voronoi_segment_formation.hpp | 1700 +++++++++++++++++++++++----------------
sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_segment_builder.hpp | 293 ++----
sandbox/SOC/2010/sweepline/libs/sweepline/test/Jamfile.v2 | 1
sandbox/SOC/2010/sweepline/libs/sweepline/test/voronoi_segment/segment_event_types_test.cpp | 10
sandbox/SOC/2010/sweepline/libs/sweepline/test/voronoi_segment/segment_predicates_test.cpp | 242 ++--
sandbox/SOC/2010/sweepline/libs/sweepline/test/voronoi_segment/segment_voronoi_builder_test.cpp | 50 +
6 files changed, 1244 insertions(+), 1052 deletions(-)
Modified: sandbox/SOC/2010/sweepline/boost/sweepline/detail/voronoi_segment_formation.hpp
==============================================================================
--- sandbox/SOC/2010/sweepline/boost/sweepline/detail/voronoi_segment_formation.hpp (original)
+++ sandbox/SOC/2010/sweepline/boost/sweepline/detail/voronoi_segment_formation.hpp 2010-09-28 15:28:09 EDT (Tue, 28 Sep 2010)
@@ -27,33 +27,17 @@
namespace detail {
///////////////////////////////////////////////////////////////////////////
- // GEOMETRY PREDICATES ////////////////////////////////////////////////////
+ // VORONOI EVENT TYPES ////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////
- // 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.
- bool almost_equal(double a, double b, long long maxUlps) {
- long long ll_a, ll_b;
- // Reinterpret double bits as long long.
- memcpy(&ll_a, &a, sizeof(double));
- memcpy(&ll_b, &b, sizeof(double));
+ template <typename T>
+ struct beach_line_node;
- // 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.
- if (ll_a < 0)
- ll_a = 0x8000000000000000LL - ll_a;
- if (ll_b < 0)
- ll_b = 0x8000000000000000LL - ll_b;
+ template <typename T>
+ struct beach_line_node_data;
- // Compare long long 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);
- }
+ template <typename BeachLineNode>
+ struct node_comparer;
enum kOrientation {
RIGHT_ORIENTATION = -1,
@@ -61,825 +45,1082 @@
LEFT_ORIENTATION = 1,
};
- // TODO(asydorchuk): Make templates specification for integer coordinate type,
- // as it is actually working with integer input.
+ // Site event type.
+ // Occurs when sweepline sweeps over one of the initial sites.
+ // Contains index of a site below the other sorted sites.
template <typename T>
- kOrientation orientation_test(const point_2d<T> &point1,
- const point_2d<T> &point2,
- const point_2d<T> &point3) {
- typedef long long ll;
- typedef unsigned long long ull;
- ull dif_x1, dif_x2, dif_y1, 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);
- 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;
+ struct site_event {
+ public:
+ typedef T coordinate_type;
+ typedef point_2d<T> Point2D;
- if (expr_l == 0)
- expr_l_plus = true;
- if (expr_r == 0)
- expr_r_plus = true;
+ site_event() : is_segment_(false),
+ is_vertical_(true),
+ is_inverse_(false) {}
- if ((expr_l_plus == expr_r_plus) && (expr_l == expr_r))
- return COLINEAR;
+ site_event(coordinate_type x, coordinate_type y, int index) :
+ point0_(x, y),
+ site_index_(index),
+ is_segment_(false),
+ is_vertical_(true),
+ is_inverse_(false) {}
- 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;
- }
- }
+ site_event(const Point2D &point, int index) :
+ point0_(point),
+ site_index_(index),
+ is_segment_(false),
+ is_vertical_(true),
+ is_inverse_(false) {}
- kOrientation orientation_test(long long dif_x1_, long long dif_y1_,
- long long dif_x2_, long long 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;
+ site_event(const Point2D &point1, const Point2D &point2, int index) :
+ point0_(point1),
+ point1_(point2),
+ site_index_(index),
+ is_segment_(true),
+ is_inverse_(false) {
+ if (point0_ > point1_)
+ (std::swap)(point0_, point1_);
+ is_vertical_ = (point0_.x() == point1_.x());
+ }
- 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;
+ 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_)));
+ }
- 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;
+ bool operator!=(const site_event &s_event) const {
+ return !((*this) == s_event);
}
- }
- enum kPredicateResult {
- LESS = -1,
- UNDEFINED = 0,
- MORE = 1,
- };
+ // 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();
- // 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),
- // 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
- // 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>
- 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());
- double fast_left_expr = fast_a1 * fast_b2 * fast_b2;
- double fast_right_expr = fast_a2 * fast_b1 * fast_b1;
-
- // Avoid cancellation.
- INT_PREDICATE_AVOID_CANCELLATION(fast_a1 * fast_a2 * fast_c,
- fast_left_expr, fast_right_expr);
- if (!almost_equal(fast_left_expr, fast_right_expr, 5))
- return (fast_left_expr < fast_right_expr) ? LESS : MORE;
- return UNDEFINED;
- }
+ if (!(this->is_segment_)) {
+ if (!s_event.is_segment_) {
+ return this->point0_.y() < s_event.point0_.y();
+ }
+ if (s_event.is_vertical_) {
+ return this->point0_.y() <= s_event.point0_.y();
+ }
+ return true;
+ } else {
+ if (!s_event.is_segment_ || s_event.is_vertical_) {
+ if (this->is_vertical_) {
+ return this->point0_.y() < s_event.point0_.y();
+ }
+ return false;
+ }
+ if (this->is_vertical_)
+ 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;
+ }
+ }
- template <typename T>
- bool less_predicate(const point_2d<T> &left_point,
- const point_2d<T> &right_point,
- const point_2d<T> &new_point) {
- kPredicateResult fast_res = fast_less_predicate(left_point, right_point, new_point);
- if (fast_res != UNDEFINED)
- return (fast_res == LESS);
+ bool operator<=(const site_event &s_event) const {
+ return !(s_event < (*this));
+ }
- 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;
+ bool operator>(const site_event &s_event) const {
+ return s_event < (*this);
+ }
- // a1 and a2 are greater than zero.
- 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()));
+ bool operator>=(const site_event &s_event) const {
+ return !((*this) < s_event);
+ }
- // 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);
- b1_sqr = b1 * b1;
- b2_sqr = b2 * b2;
- ull b1_sqr_mod = b1_sqr % a1;
- ull b2_sqr_mod = b2_sqr % a2;
+ coordinate_type x(bool oriented = false) const {
+ if (!oriented)
+ return point0_.x();
+ return is_inverse_ ? point1_.x() : point0_.x();
+ }
- // Compute left expression.
- INT_PREDICATE_COMPUTE_DIFFERENCE(static_cast<ll>(left_point.x()),
- static_cast<ll>(right_point.x()),
- l_expr, l_expr_plus);
- if (b2_sqr_mod * a1 < b1_sqr_mod * a2) {
- if (!l_expr_plus)
- l_expr++;
- else if (l_expr != 0)
- l_expr--;
- else {
- l_expr++;
- l_expr_plus = false;
- }
+ coordinate_type y(bool oriented = false) const {
+ if (!oriented)
+ return point0_.y();
+ return is_inverse_ ? point1_.y() : point0_.y();
}
- // Compute right expression.
- INT_PREDICATE_COMPUTE_DIFFERENCE(b1_sqr / a1, b2_sqr / a2, r_expr, r_expr_plus);
+ coordinate_type x0(bool oriented = false) const {
+ if (!oriented)
+ return point0_.x();
+ return is_inverse_ ? point1_.x() : point0_.x();
+ }
- // Compare left and right expressions.
- if (!l_expr_plus && r_expr_plus)
- return true;
- if (l_expr_plus && !r_expr_plus)
- return false;
- if (l_expr_plus && r_expr_plus)
- return l_expr < r_expr;
- return l_expr > r_expr;
- }
+ coordinate_type y0(bool oriented = false) const {
+ if (!oriented)
+ return point0_.y();
+ return is_inverse_ ? point1_.y() : point0_.y();
+ }
- template <typename T>
- kPredicateResult less_predicate_check(point_2d<T> segment_start, point_2d<T> segment_end,
- point_2d<T> site_point, point_2d<T> new_point,
- bool reverse_order) {
- typedef long long ll;
- typedef unsigned long long ull;
- kOrientation orientation1 = orientation_test(segment_start, segment_end, site_point);
- kOrientation orientation2 = orientation_test(segment_start, segment_end, new_point);
- if (orientation1 == COLINEAR)
- orientation1 = reverse_order ? LEFT_ORIENTATION : RIGHT_ORIENTATION;
- if (orientation1 != orientation2)
- return (orientation1 == RIGHT_ORIENTATION) ? LESS : MORE;
-
- 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());
- ll b = static_cast<ll>(segment_end.y()) - static_cast<ll>(segment_start.y());
-
- if (a == 0) {
- if (new_point.y() < site_point.y() && !reverse_order)
- return MORE;
- else if (new_point.y() > site_point.y() && reverse_order)
- return LESS;
- return UNDEFINED;
- } else {
- kOrientation orientation = orientation_test(a, b, dif_x, dif_y);
- if ((orientation == COLINEAR) ||
- ((orientation1 == RIGHT_ORIENTATION) ^ (orientation == RIGHT_ORIENTATION))) {
- if (orientation1 == RIGHT_ORIENTATION)
- return reverse_order ? LESS : UNDEFINED;
- return reverse_order ? UNDEFINED : MORE;
- }
+ coordinate_type x1(bool oriented = false) const {
+ if (!oriented)
+ return point1_.x();
+ return is_inverse_ ? point0_.x() : point1_.x();
}
- // 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);
- double fast_right_expr = static_cast<double>(b<<1) *
- static_cast<double>(dif_x) *
- static_cast<double>(dif_y);
- if (!(almost_equal(fast_left_expr, fast_right_expr, 4))) {
- if ((orientation1 == RIGHT_ORIENTATION) ^
- (fast_left_expr > fast_right_expr) ^ !reverse_order)
- return reverse_order ? LESS : MORE;
- return UNDEFINED;
+ coordinate_type y1(bool oriented = false) const {
+ if (!oriented)
+ return point1_.y();
+ return is_inverse_ ? point0_.y() : point1_.y();
}
- 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);
-
- 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;
- ull left_expr_mod = dif_y_sqr % dif_x_sqr;
-
- ull right_expr_denom = a_rob * dif_x_rob;
- ull right_expr = b_rob * dif_y_rob;
- ull right_expr_div = right_expr / right_expr_denom;
- ull right_expr_mod = right_expr % right_expr_denom;
-
- bool comparison_result;
- if ((b_sign != dif_y_sign) && right_expr_div)
- comparison_result = true;
- else {
- if (b_sign != dif_y_sign && right_expr_mod)
- right_expr_mod = right_expr_denom - right_expr_mod;
- else
- right_expr_div++;
- bool left_expr_odd = (left_expr_div % 2 == 1);
- left_expr_div >>= 1;
- if (left_expr_div != right_expr_div) {
- comparison_result = (left_expr_div > right_expr_div);
- } else {
- ull temp_right = right_expr_denom - right_expr_mod;
- if (temp_right > right_expr_mod) {
- if (left_expr_odd)
- comparison_result = true;
- else
- right_expr_mod <<= 1;
- } else {
- if (!left_expr_odd)
- comparison_result = false;
- else
- right_expr_mod = right_expr_mod - temp_right;
- }
- left_expr_div = left_expr_mod / dif_x_rob;
- right_expr_div = right_expr_mod / a_rob;
- if (left_expr_div != right_expr_div)
- comparison_result = (left_expr_div > right_expr_div);
- else {
- left_expr_mod = left_expr_mod % dif_x_rob;
- right_expr_mod = right_expr_mod % a_rob;
- comparison_result = (left_expr_mod * a_rob > right_expr_mod * dif_x_rob);
- }
- }
+ const Point2D &get_point0(bool oriented = false) const {
+ if (!oriented)
+ return point0_;
+ return is_inverse_ ? point1_ : point0_;
}
-
- if ((orientation1 == RIGHT_ORIENTATION) ^ comparison_result ^ !reverse_order)
- return reverse_order ? LESS : MORE;
- return UNDEFINED;
- }
-
-#ifdef USE_MULTI_PRECISION_LIBRARY
- template<typename T>
- bool mpz_less_predicate(point_2d<T> segment_start, point_2d<T> segment_end,
- point_2d<T> site_point, point_2d<T> new_point, bool reverse_order) {
- mpz_class segment_start_x, segment_start_y, segment_end_x, segment_end_y,
- site_point_x, site_point_y, new_point_x, new_point_y;
- segment_start_x = static_cast<int>(segment_start.x());
- segment_start_y = static_cast<int>(segment_start.y());
- segment_end_x = static_cast<int>(segment_end.x());
- segment_end_y = static_cast<int>(segment_end.y());
- site_point_x = static_cast<int>(site_point.x());
- site_point_y = static_cast<int>(site_point.y());
- new_point_x = static_cast<int>(new_point.x());
- new_point_y = static_cast<int>(new_point.y());
-
- mpz_class dif_x, dif_y, a, b, mul1, mul2, temp, left_expr, right_expr;
- dif_x = new_point_x - site_point_x;
- dif_y = new_point_y - site_point_y;
- a = segment_end_x - segment_start_x;
- b = segment_end_y - segment_start_y;
- mul1 = new_point_x - segment_end_x;
- mul2 = new_point_y - segment_end_y;
- temp = dif_x * dif_x + dif_y * dif_y;
- left_expr = (a * a + b * b) * temp * temp;
- right_expr = (2.0 * dif_x * (b * mul1 - a * mul2) - b * temp);
- right_expr = right_expr * right_expr;
-
- return (!reverse_order) ? (left_expr > right_expr) : (left_expr < right_expr);
- }
-#endif
-
- // 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.
- // 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.
- template <typename T>
- bool less_predicate(point_2d<T> segment_start, point_2d<T> segment_end,
- point_2d<T> site_point, point_2d<T> new_point, bool reverse_order) {
- kPredicateResult fast_res = less_predicate_check(segment_start, segment_end,
- site_point, new_point, reverse_order);
- if (fast_res != UNDEFINED)
- return (fast_res == LESS);
-
- double dif_x = static_cast<double>(new_point.x()) -
- static_cast<double>(site_point.x());
- double dif_y = static_cast<double>(new_point.y()) -
- static_cast<double>(site_point.y());
- double a = static_cast<double>(segment_end.x()) -
- static_cast<double>(segment_start.x());
- double b = static_cast<double>(segment_end.y()) -
- static_cast<double>(segment_start.y());
- double mul1 = static_cast<double>(new_point.x()) - static_cast<double>(segment_end.x());
- double mul2 = static_cast<double>(new_point.y()) - static_cast<double>(segment_end.y());
- double temp = dif_x * dif_x + dif_y * dif_y;
- double a_sqr = a * a;
- double b_sqr = b * b;
- double left_expr = (a_sqr * temp + (4.0 * dif_x * mul1) * b_sqr) * temp;
- double right_expr = (4.0 * dif_x * dif_x) * ((mul1 * mul1) * b_sqr + (mul2 * mul2) * a_sqr);
- double common_expr = (4.0 * dif_x * a) * (b * mul2);
- INT_PREDICATE_AVOID_CANCELLATION(common_expr * temp, right_expr, left_expr);
- INT_PREDICATE_AVOID_CANCELLATION(common_expr * (2.0 * dif_x * mul1), left_expr, right_expr);
+ const Point2D &get_point1(bool oriented = false) const {
+ if (!oriented)
+ return point1_;
+ return is_inverse_ ? point0_ : point1_;
+ }
- if (!almost_equal(left_expr, right_expr, 18)) {
- if (!reverse_order)
- return left_expr > right_expr;
- return left_expr < right_expr;
+ void set_site_index(int index) {
+ site_index_ = index;
}
-#ifndef USE_MULTI_PRECISION_LIBRARY
- return (!reverse_order) ? (left_expr > right_expr) : (left_expr < right_expr);
-#else
- return mpz_less_predicate(segment_start, segment_end, site_point,
- new_point, reverse_order);
-#endif
- }
+ void set_inverse() {
+ is_inverse_ ^= true;
+ }
- template <typename T>
- bool less_predicate(point_2d<T> segment1_start, point_2d<T> segment1_end,
- point_2d<T> segment2_start, point_2d<T> segment2_end,
- point_2d<T> new_point, bool reverse_order) {
- typedef unsigned long long ull;
- kOrientation orientation1 = orientation_test(segment1_start, segment1_end, segment2_start);
- kOrientation orientation2 = orientation_test(segment1_start, segment1_end, segment2_end);
- assert(static_cast<int>(orientation1) || static_cast<int>(orientation2));
- if (static_cast<int>(orientation1) * static_cast<int>(orientation2) == -1) {
- assert(!reverse_order);
- return less_predicate(segment2_start, segment2_end,
- segment1_start, segment1_end,
- new_point, true);
+ int get_site_index() const {
+ return site_index_;
}
- double intersection_x1 = 0.0;
- double intersection_x2 = 0.0;
- double a1 = static_cast<double>(segment1_start.x()) -
- static_cast<double>(segment1_end.x());
- if (a1 == 0.0) {
- // Avoid cancellation.
- intersection_x2 += static_cast<double>(new_point.x()) -
- static_cast<double>(segment1_end.x());
- } else {
- double b1 = static_cast<double>(segment1_start.y()) -
- static_cast<double>(segment1_end.y());
- double c1 = b1 * (static_cast<double>(new_point.x()) -
- static_cast<double>(segment1_end.x())) +
- a1 * segment1_end.y();
- double mul1 = sqrt(a1 * a1 + b1 * b1);
- if ((orientation1 == LEFT_ORIENTATION) || (orientation2 == LEFT_ORIENTATION)) {
- if (b1 >= 0.0) {
- mul1 = 1.0 / (b1 + mul1);
- } else {
- mul1 = (-b1 + mul1) / (a1 * a1);
- }
- } else {
- if (b1 >= 0.0) {
- mul1 = (-b1 - mul1) / (a1 * a1);
- } else {
- mul1 = 1.0 / (b1 - mul1);
- }
- }
- INT_PREDICATE_AVOID_CANCELLATION(a1 * mul1 * static_cast<double>(new_point.y()),
- intersection_x1, intersection_x2);
- INT_PREDICATE_AVOID_CANCELLATION(-c1 * mul1, intersection_x1, intersection_x2);
+ bool is_segment() const {
+ return is_segment_;
}
- double a2 = static_cast<double>(segment2_start.x()) -
- static_cast<double>(segment2_end.x());
- if (a2 == 0.0) {
- // Avoid cancellation.
- intersection_x1 += static_cast<double>(new_point.x()) -
- static_cast<double>(segment2_end.x());
- } else {
- double b2 = static_cast<double>(segment2_start.y()) -
- static_cast<double>(segment2_end.y());
- double c2 = b2 * (static_cast<double>(new_point.x()) -
- static_cast<double>(segment2_end.x())) +
- a2 * segment2_end.y();
- double mul2 = sqrt(a2 * a2 + b2 * b2);
- if (((orientation1 == LEFT_ORIENTATION) || (orientation2 == LEFT_ORIENTATION)) ^
- !reverse_order) {
- 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);
- }
- }
- INT_PREDICATE_AVOID_CANCELLATION(a2 * static_cast<double>(new_point.y()) * mul2,
- intersection_x2, intersection_x1);
- INT_PREDICATE_AVOID_CANCELLATION(-c2 * mul2, intersection_x2, intersection_x1);
+ bool is_vertical() const {
+ return is_vertical_;
}
- if (!almost_equal(intersection_x1, intersection_x2, 20)) {
- return ((orientation1 == LEFT_ORIENTATION) || (orientation2 == LEFT_ORIENTATION)) ^
- reverse_order ^ (intersection_x1 > intersection_x2);
+ bool is_inverse() const {
+ return is_inverse_;
}
-
- // TODO(asydorchuk): Add mpl support there.
- return ((orientation1 == LEFT_ORIENTATION) || (orientation2 == LEFT_ORIENTATION)) ^
- reverse_order ^ (intersection_x1 > intersection_x2);
- }
- ///////////////////////////////////////////////////////////////////////////
- // VORONOI EVENT TYPES ////////////////////////////////////////////////////
- ///////////////////////////////////////////////////////////////////////////
+ private:
+ Point2D point0_;
+ Point2D point1_;
+ int site_index_;
+ bool is_segment_;
+ bool is_vertical_;
+ bool is_inverse_;
+ };
template <typename T>
- struct beach_line_node;
+ site_event<T> make_site_event(T x, T y, int index) {
+ return site_event<T>(x, y, index);
+ }
template <typename T>
- struct beach_line_node_data;
+ site_event<T> make_site_event(const point_2d<T> &point, int index) {
+ return site_event<T>(point, index);
+ }
- template <typename BeachLineNode>
- struct node_comparer;
+ template <typename T>
+ 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);
+ }
- // Site event type.
- // Occurs when sweepline sweeps over one of the initial sites.
- // Contains index of a site below the other sorted sites.
+ // 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.
template <typename T>
- struct site_event {
+ struct circle_event {
public:
typedef T coordinate_type;
typedef point_2d<T> Point2D;
+ 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;
- site_event() : is_segment_(false), is_vertical_(true) {}
-
- site_event(coordinate_type x, coordinate_type y, int index) :
- point0_(x, y), site_index_(index), is_segment_(false), is_vertical_(true) {}
+ circle_event() : is_active_(true) {}
- site_event(const Point2D &point, int index) :
- point0_(point), site_index_(index), is_segment_(false), is_vertical_(true) {}
+ circle_event(coordinate_type c_x, coordinate_type c_y, coordinate_type lower_x) :
+ center_(c_x, c_y), lower_x_(lower_x), is_active_(true) {}
- site_event(const Point2D &point1, const Point2D &point2, int index) :
- point0_(point1), point1_(point2), site_index_(index), is_segment_(true) {
- if (point0_ > point1_)
- (std::swap)(point0_, point1_);
- is_vertical_ = (point0_.x() == point1_.x());
+ circle_event(const Point2D ¢er, coordinate_type lower_x) :
+ center_(center), lower_x_(lower_x), is_active_(true) {}
+
+ circle_event(const circle_event& c_event) {
+ center_ = c_event.center_;
+ lower_x_ = c_event.lower_x_;
+ bisector_node_ = c_event.bisector_node_;
+ is_active_ = c_event.is_active_;
}
- 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_)));
+ void operator=(const circle_event& c_event) {
+ center_ = c_event.center_;
+ lower_x_ = c_event.lower_x_;
+ bisector_node_ = c_event.bisector_node_;
+ is_active_ = c_event.is_active_;
}
- bool operator!=(const site_event &s_event) const {
- return !((*this) == s_event);
+ bool operator==(const circle_event &c_event) const {
+ return (center_.y() == c_event.y()) &&
+ (lower_x_ == c_event.lower_x_);
}
- //// All the sites are sorted due to coordinates of the first point.
- //// Point sites preceed segment sites with the same first point.
- //// Segments with the same first point use counterclockwise comparison.
- //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();
-
- // // If first points have different y coordinates, compare them.
- // if (this->point0_.y() != s_event.point0_.y() ||
- // !(this->is_segment_ || s_event.is_segment_))
- // return this->point0_.y() < s_event.point0_.y();
-
- // // Point site events go before segment site events with the same first point.
- // if (this->is_segment_ ^ s_event.is_segment_)
- // return !this->is_segment;
-
- // // Use angle comparison for segment sites with the same first point.
- // return orientation_test(this->point1_, this->point0_, s_event.point1_) ==
- // RIGHT_ORIENTATION;
- //}
+ bool operator!=(const circle_event &c_event) const {
+ return !((*this) == c_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 circle_event &c_event) const {
+ if (lower_x_ != c_event.lower_x_)
+ return lower_x_ < c_event.lower_x_;
+ return center_.y() < c_event.y();
+ }
- if (!(this->is_segment_)) {
- if (!s_event.is_segment_) {
- return this->point0_.y() < s_event.point0_.y();
- }
- if (s_event.is_vertical_) {
- return this->point0_.y() <= s_event.point0_.y();
- }
- return true;
- } else {
- if (!s_event.is_segment_ || s_event.is_vertical_) {
- if (this->is_vertical_) {
- return this->point0_.y() < s_event.point0_.y();
- }
- return false;
- }
- if (this->is_vertical_)
- 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_) ==
- RIGHT_ORIENTATION;
- }
+ bool operator<=(const circle_event &c_event) const {
+ return !(c_event < (*this));
}
- bool operator<=(const site_event &s_event) const {
- return !(s_event < (*this));
+ bool operator>(const circle_event &c_event) const {
+ return c_event < (*this);
}
- bool operator>(const site_event &s_event) const {
- return s_event < (*this);
+ bool operator>=(const circle_event &c_event) const {
+ return !((*this) < c_event);
}
- bool operator>=(const site_event &s_event) const {
- return !((*this) < s_event);
+ // Compares bottom voronoi circle point with site event point.
+ // 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.
+ int compare(const site_event<coordinate_type> &s_event) const {
+ if (s_event.x() != lower_x_) {
+ return (lower_x_ < s_event.x()) ? -1 : 1;
+ }
+ if (s_event.y() != center_.y())
+ return (center_.y() < s_event.y()) ? -1 : 1;
+ return 0;
}
coordinate_type x() const {
- return point0_.x();
+ return center_.x();
}
coordinate_type y() const {
- return point0_.y();
+ return center_.y();
}
- const Point2D &get_point0() const {
- return point0_;
+ const Point2D &get_center() const {
+ return center_;
}
- const Point2D &get_point1() const {
- return point1_;
+ const coordinate_type &get_lower_x() const {
+ return lower_x_;
}
- void set_site_index(int index) {
- site_index_ = index;
+ void set_bisector(beach_line_iterator iterator) {
+ bisector_node_ = iterator;
}
- int get_site_index() const {
- return site_index_;
+ void deactivate() {
+ is_active_ = false;
}
- bool is_segment() const {
- return is_segment_;
+ const beach_line_iterator &get_bisector() const {
+ return bisector_node_;
}
- bool is_vertical() const {
- return is_vertical_;
+ bool is_active() const {
+ return is_active_;
}
private:
- Point2D point0_;
- Point2D point1_;
- int site_index_;
- bool is_segment_;
- bool is_vertical_;
+ Point2D center_;
+ coordinate_type lower_x_;
+ beach_line_iterator bisector_node_;
+ bool is_active_;
};
template <typename T>
- site_event<T> make_site_event(T x, T y, int index) {
- return site_event<T>(x, y, index);
+ 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);
}
template <typename T>
- 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(const point_2d<T> &point1,
- const point_2d<T> &point2, int index) {
- return site_event<T>(point1, point2, index);
+ circle_event<T> make_circle_event(point_2d<T> center, T lower_x) {
+ return circle_event<T>(center, lower_x);
}
- // 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.
+ ///////////////////////////////////////////////////////////////////////////
+ // VORONOI CIRCLE EVENTS QUEUE ////////////////////////////////////////////
+ ///////////////////////////////////////////////////////////////////////////
+
+ // Event queue data structure, processes circle events.
template <typename T>
- struct circle_event {
+ class circle_events_queue {
public:
typedef T coordinate_type;
typedef point_2d<T> Point2D;
- 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 circle_event<T> circle_event_type;
+ typedef typename std::list<circle_event_type>::iterator circle_event_type_it;
- circle_event() : is_active_(true) {}
+ circle_events_queue() {}
- circle_event(coordinate_type c_x, coordinate_type c_y, coordinate_type lower_x) :
- center_(c_x, c_y), lower_x_(lower_x), is_active_(true) {}
+ void reset() {
+ while (!circle_events_.empty())
+ circle_events_.pop();
+ circle_events_list_.clear();
+ }
- circle_event(const Point2D ¢er, coordinate_type lower_x) :
- center_(center), lower_x_(lower_x), is_active_(true) {}
+ bool empty() {
+ remove_not_active_events();
+ return circle_events_.empty();
+ }
- circle_event(const circle_event& c_event) {
- center_ = c_event.center_;
- lower_x_ = c_event.lower_x_;
- bisector_node_ = c_event.bisector_node_;
- is_active_ = c_event.is_active_;
+ const circle_event_type &top() {
+ remove_not_active_events();
+ return (*circle_events_.top());
}
- void operator=(const circle_event& c_event) {
- center_ = c_event.center_;
- lower_x_ = c_event.lower_x_;
- bisector_node_ = c_event.bisector_node_;
- is_active_ = c_event.is_active_;
+ void pop() {
+ circle_event_type_it temp_it = circle_events_.top();
+ circle_events_.pop();
+ circle_events_list_.erase(temp_it);
}
- bool operator==(const circle_event &c_event) const {
- return (center_.y() == c_event.y()) &&
- (lower_x_ == c_event.lower_x_);
+ circle_event_type_it push(const circle_event_type &c_event) {
+ circle_events_list_.push_front(c_event);
+ circle_events_.push(circle_events_list_.begin());
+ return circle_events_list_.begin();
}
- bool operator!=(const circle_event &c_event) const {
- return !((*this) == c_event);
+ private:
+ struct comparison {
+ bool operator() (const circle_event_type_it &node1,
+ const circle_event_type_it &node2) const {
+ return (*node1) > (*node2);
+ }
+ };
+
+ void remove_not_active_events() {
+ while (!circle_events_.empty() && !circle_events_.top()->is_active())
+ pop();
}
- bool operator<(const circle_event &c_event) const {
- if (lower_x_ != c_event.lower_x_)
- return lower_x_ < c_event.lower_x_;
- return center_.y() < c_event.y();
+ std::priority_queue< circle_event_type_it,
+ std::vector<circle_event_type_it>,
+ comparison > circle_events_;
+ std::list<circle_event_type> circle_events_list_;
+
+ //Disallow copy constructor and operator=
+ circle_events_queue(const circle_events_queue&);
+ void operator=(const circle_events_queue&);
+ };
+
+ ///////////////////////////////////////////////////////////////////////////
+ // GEOMETRY PREDICATES ////////////////////////////////////////////////////
+ ///////////////////////////////////////////////////////////////////////////
+
+ // 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.
+ bool almost_equal(double a, double b, long long maxUlps) {
+ long long ll_a, ll_b;
+ // Reinterpret double bits as long long.
+ 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.
+ 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.
+ // 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.
+ template <typename T>
+ kOrientation orientation_test(const point_2d<T> &point1,
+ const point_2d<T> &point2,
+ const point_2d<T> &point3) {
+ typedef long long ll;
+ typedef unsigned long long ull;
+ ull dif_x1, dif_x2, dif_y1, 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);
+ 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;
+ }
+ }
+
+ kOrientation orientation_test(long long dif_x1_, long long dif_y1_,
+ long long dif_x2_, long long 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;
+ }
+ }
+
+ enum kPredicateResult {
+ LESS = -1,
+ UNDEFINED = 0,
+ MORE = 1,
+ };
+
+ // 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),
+ // 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
+ // 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>
+ 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());
+ double fast_left_expr = fast_a1 * fast_b2 * fast_b2;
+ double fast_right_expr = fast_a2 * fast_b1 * fast_b1;
+
+ // Avoid cancellation.
+ INT_PREDICATE_AVOID_CANCELLATION(fast_a1 * fast_a2 * fast_c,
+ fast_left_expr, fast_right_expr);
+ if (!almost_equal(fast_left_expr, fast_right_expr, 5))
+ return (fast_left_expr < fast_right_expr) ? LESS : MORE;
+ return UNDEFINED;
+ }
+
+ template <typename T>
+ bool less_predicate(const point_2d<T> &left_point,
+ const point_2d<T> &right_point,
+ const point_2d<T> &new_point) {
+ kPredicateResult fast_res = fast_less_predicate(left_point, right_point, new_point);
+ if (fast_res != UNDEFINED)
+ return (fast_res == LESS);
+
+ 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.
+ 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);
+ 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);
+ if (b2_sqr_mod * a1 < b1_sqr_mod * a2) {
+ if (!l_expr_plus)
+ l_expr++;
+ else if (l_expr != 0)
+ l_expr--;
+ else {
+ l_expr++;
+ l_expr_plus = false;
+ }
+ }
+
+ // Compute right expression.
+ INT_PREDICATE_COMPUTE_DIFFERENCE(b1_sqr / a1, b2_sqr / a2, r_expr, r_expr_plus);
+
+ // Compare left and right expressions.
+ if (!l_expr_plus && r_expr_plus)
+ return true;
+ if (l_expr_plus && !r_expr_plus)
+ return false;
+ if (l_expr_plus && r_expr_plus)
+ return l_expr < r_expr;
+ return l_expr > r_expr;
+ }
+
+ template <typename T>
+ 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),
+ 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();
+ 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());
+ ll b = static_cast<ll>(segment_end.y()) - static_cast<ll>(segment_start.y());
+
+ if (a == 0) {
+ if (new_point.y() < site_point.y() && !reverse_order)
+ return MORE;
+ else if (new_point.y() > site_point.y() && reverse_order)
+ return LESS;
+ return UNDEFINED;
+ } else {
+ kOrientation orientation = orientation_test(a, b, dif_x, dif_y);
+ if ((orientation == COLINEAR) ||
+ (!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);
+ double fast_right_expr = static_cast<double>(b<<1) *
+ static_cast<double>(dif_x) *
+ static_cast<double>(dif_y);
+ if (!(almost_equal(fast_left_expr, fast_right_expr, 4))) {
+ if (segment_site.is_inverse() ^ (fast_left_expr > fast_right_expr) ^ reverse_order)
+ return reverse_order ? LESS : MORE;
+ return UNDEFINED;
+ }
+
+ 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);
+
+ 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;
+ ull left_expr_mod = dif_y_sqr % dif_x_sqr;
+
+ ull right_expr_denom = a_rob * dif_x_rob;
+ ull right_expr = b_rob * dif_y_rob;
+ ull right_expr_div = right_expr / right_expr_denom;
+ ull right_expr_mod = right_expr % right_expr_denom;
+
+ bool comparison_result;
+ if ((b_sign != dif_y_sign) && right_expr_div)
+ comparison_result = true;
+ else {
+ if (b_sign != dif_y_sign && right_expr_mod)
+ right_expr_mod = right_expr_denom - right_expr_mod;
+ else
+ right_expr_div++;
+ bool left_expr_odd = (left_expr_div % 2 == 1);
+ left_expr_div >>= 1;
+ if (left_expr_div != right_expr_div) {
+ comparison_result = (left_expr_div > right_expr_div);
+ } else {
+ ull temp_right = right_expr_denom - right_expr_mod;
+ if (temp_right > right_expr_mod) {
+ if (left_expr_odd)
+ comparison_result = true;
+ else
+ right_expr_mod <<= 1;
+ } else {
+ if (!left_expr_odd)
+ comparison_result = false;
+ else
+ right_expr_mod = right_expr_mod - temp_right;
+ }
+ left_expr_div = left_expr_mod / dif_x_rob;
+ right_expr_div = right_expr_mod / a_rob;
+ if (left_expr_div != right_expr_div)
+ comparison_result = (left_expr_div > right_expr_div);
+ else {
+ left_expr_mod = left_expr_mod % dif_x_rob;
+ right_expr_mod = right_expr_mod % a_rob;
+ comparison_result = (left_expr_mod * a_rob > right_expr_mod * dif_x_rob);
+ }
+ }
}
+
+ if (segment_site.is_inverse() ^ comparison_result ^ reverse_order)
+ return reverse_order ? LESS : MORE;
+ return UNDEFINED;
+ }
- bool operator<=(const circle_event &c_event) const {
- return !(c_event < (*this));
- }
+#ifdef USE_MULTI_PRECISION_LIBRARY
+ template<typename T>
+ bool mpz_less_predicate(point_2d<T> segment_start, point_2d<T> segment_end,
+ point_2d<T> site_point, point_2d<T> new_point, bool reverse_order) {
+ mpz_class segment_start_x, segment_start_y, segment_end_x, segment_end_y,
+ site_point_x, site_point_y, new_point_x, new_point_y;
+ segment_start_x = static_cast<int>(segment_start.x());
+ segment_start_y = static_cast<int>(segment_start.y());
+ segment_end_x = static_cast<int>(segment_end.x());
+ segment_end_y = static_cast<int>(segment_end.y());
+ site_point_x = static_cast<int>(site_point.x());
+ site_point_y = static_cast<int>(site_point.y());
+ new_point_x = static_cast<int>(new_point.x());
+ new_point_y = static_cast<int>(new_point.y());
- bool operator>(const circle_event &c_event) const {
- return c_event < (*this);
- }
+ mpz_class dif_x, dif_y, a, b, mul1, mul2, temp, left_expr, right_expr;
+ dif_x = new_point_x - site_point_x;
+ dif_y = new_point_y - site_point_y;
+ a = segment_end_x - segment_start_x;
+ b = segment_end_y - segment_start_y;
+ mul1 = new_point_x - segment_end_x;
+ mul2 = new_point_y - segment_end_y;
+ temp = dif_x * dif_x + dif_y * dif_y;
+ left_expr = (a * a + b * b) * temp * temp;
+ right_expr = (2.0 * dif_x * (b * mul1 - a * mul2) - b * temp);
+ right_expr = right_expr * right_expr;
- bool operator>=(const circle_event &c_event) const {
- return !((*this) < c_event);
- }
+ return (!reverse_order) ? (left_expr > right_expr) : (left_expr < right_expr);
+ }
+#endif
- // Compares bottom voronoi circle point with site event point.
- // 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.
- int compare(const site_event<coordinate_type> &s_event) const {
- if (s_event.x() != lower_x_) {
- return (lower_x_ < s_event.x()) ? -1 : 1;
- }
- if (s_event.y() != center_.y())
- return (center_.y() < s_event.y()) ? -1 : 1;
- return 0;
+ // 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.
+ // 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.
+ template <typename T>
+ bool less_predicate(point_2d<T> site_point, site_event<T> segment_site,
+ point_2d<T> new_point, bool reverse_order) {
+ kPredicateResult fast_res = fast_less_predicate(
+ site_point, segment_site, new_point, reverse_order);
+ if (fast_res != UNDEFINED) {
+ return (fast_res == LESS);
}
- coordinate_type x() const {
- return center_.x();
- }
+ const point_2d<T> &segment_start = segment_site.get_point0();
+ const point_2d<T> &segment_end = segment_site.get_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()) -
+ static_cast<double>(site_point.y());
+ double a = static_cast<double>(segment_end.x()) -
+ static_cast<double>(segment_start.x());
+ double b = static_cast<double>(segment_end.y()) -
+ static_cast<double>(segment_start.y());
- coordinate_type y() const {
- return center_.y();
- }
+ double mul1 = static_cast<double>(new_point.x()) - static_cast<double>(segment_end.x());
+ double mul2 = static_cast<double>(new_point.y()) - static_cast<double>(segment_end.y());
+ double temp = dif_x * dif_x + dif_y * dif_y;
+ double a_sqr = a * a;
+ double b_sqr = b * b;
+ double left_expr = (a_sqr * temp + (4.0 * dif_x * mul1) * b_sqr) * temp;
+ double right_expr = (4.0 * dif_x * dif_x) * ((mul1 * mul1) * b_sqr + (mul2 * mul2) * a_sqr);
+ double common_expr = (4.0 * dif_x * a) * (b * mul2);
+ INT_PREDICATE_AVOID_CANCELLATION(common_expr * temp, right_expr, left_expr);
+ INT_PREDICATE_AVOID_CANCELLATION(common_expr * (2.0 * dif_x * mul1), left_expr, right_expr);
- const Point2D &get_center() const {
- return center_;
+ if (!almost_equal(left_expr, right_expr, 18)) {
+ if (!reverse_order)
+ return left_expr > right_expr;
+ return left_expr < right_expr;
}
- const coordinate_type &get_lower_x() const {
- return lower_x_;
- }
+#ifndef USE_MULTI_PRECISION_LIBRARY
+ return (!reverse_order) ? (left_expr > right_expr) : (left_expr < right_expr);
+#else
+ return mpz_less_predicate(segment_start, segment_end, site_point,
+ new_point, reverse_order);
+#endif
+ }
- void set_bisector(beach_line_iterator iterator) {
- bisector_node_ = iterator;
+ template <typename T>
+ 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) == RIGHT_ORIENTATION;
}
- void deactivate() {
- is_active_ = 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 = static_cast<double>(segment1_end.x()) -
+ static_cast<double>(segment1_start.x());
+ if (a1 == 0.0) {
+ // Avoid cancellation.
+ intersection_x2 += (static_cast<double>(new_point.x()) -
+ static_cast<double>(segment1_end.x())) * 0.5;
+ } else {
+ double b1 = static_cast<double>(segment1_end.y()) -
+ static_cast<double>(segment1_start.y());
+ double c1 = b1 * (static_cast<double>(new_point.x()) -
+ static_cast<double>(segment1_start.x())) +
+ a1 * segment1_start.y();
+ double mul1 = sqrt(a1 * a1 + b1 * b1);
+ if (left_site.is_inverse()) {
+ if (b1 >= 0.0) {
+ mul1 = 1.0 / (b1 + mul1);
+ } else {
+ mul1 = (-b1 + mul1) / (a1 * a1);
+ }
+ } else {
+ if (b1 >= 0.0) {
+ mul1 = (-b1 - mul1) / (a1 * a1);
+ } else {
+ mul1 = 1.0 / (b1 - mul1);
+ }
+ }
+ INT_PREDICATE_AVOID_CANCELLATION(a1 * mul1 * static_cast<double>(new_point.y()),
+ intersection_x1, intersection_x2);
+ INT_PREDICATE_AVOID_CANCELLATION(-c1 * mul1, intersection_x1, intersection_x2);
}
- const beach_line_iterator &get_bisector() const {
- return bisector_node_;
+ double a2 = static_cast<double>(segment2_end.x()) -
+ static_cast<double>(segment2_start.x());
+ if (a2 == 0.0) {
+ // Avoid cancellation.
+ intersection_x1 += (static_cast<double>(new_point.x()) -
+ static_cast<double>(segment2_end.x())) * 0.5;
+ } else {
+ double b2 = static_cast<double>(segment2_end.y()) -
+ static_cast<double>(segment2_start.y());
+ double c2 = b2 * (static_cast<double>(new_point.x()) -
+ static_cast<double>(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);
+ }
+ }
+ INT_PREDICATE_AVOID_CANCELLATION(a2 * static_cast<double>(new_point.y()) * mul2,
+ intersection_x2, intersection_x1);
+ INT_PREDICATE_AVOID_CANCELLATION(-c2 * mul2, intersection_x2, intersection_x1);
}
- bool is_active() const {
- return is_active_;
+ if (!almost_equal(intersection_x1, intersection_x2, 20)) {
+ return intersection_x1 < intersection_x2;
}
-
- private:
- Point2D center_;
- coordinate_type lower_x_;
- beach_line_iterator bisector_node_;
- bool is_active_;
- };
-
- template <typename T>
- circle_event<T> make_circle_event(T c_x, T c_y, T sqr_radius) {
- return circle_event<T>(c_x, c_y, sqr_radius);
- }
-
- template <typename T>
- circle_event<T> make_circle_event(point_2d<T> center, T sqr_radius) {
- return circle_event<T>(center, sqr_radius);
+
+ // TODO(asydorchuk): Add mpl support there.
+ return intersection_x1 < intersection_x2;
}
///////////////////////////////////////////////////////////////////////////
- // VORONOI CIRCLE EVENTS QUEUE ////////////////////////////////////////////
+ // CIRCLE EVENTS CREATION /////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////
-
- // Event queue data structure, processes circle events.
- template <typename T>
- class circle_events_queue {
- public:
- typedef T coordinate_type;
- typedef point_2d<T> Point2D;
- typedef circle_event<T> circle_event_type;
- typedef typename std::list<circle_event_type>::iterator circle_event_type_it;
- circle_events_queue() {}
+ // Create circle event from three point sites.
+ // TODO (asydorchuk): make precision optimizations.
+ template <typename T>
+ static bool create_circle_event_ppp(const site_event<T> &site1,
+ const site_event<T> &site2,
+ const site_event<T> &site3,
+ circle_event<T> &c_event) {
+ double dif_x1 = static_cast<double>(site1.x()) -
+ static_cast<double>(site2.x());
+ double dif_x2 = static_cast<double>(site2.x()) -
+ static_cast<double>(site3.x());
+ double dif_y1 = static_cast<double>(site1.y()) -
+ static_cast<double>(site2.y());
+ double dif_y2 = static_cast<double>(site2.y()) -
+ static_cast<double>(site3.y());
+ kOrientation valid_orientation = orientation_test(
+ static_cast<long long>(dif_x1), static_cast<long long>(dif_y1),
+ static_cast<long long>(dif_x2), static_cast<long long>(dif_y2));
+ if (valid_orientation != RIGHT_ORIENTATION)
+ return false;
+ double a = 2.0 * (dif_x1 * dif_y2 - dif_x2 * dif_y1);
+ double b1 = dif_x1 * (site1.x() + site2.x()) + dif_y1 * (site1.y() + site2.y());
+ double b2 = dif_x2 * (site2.x() + site3.x()) + dif_y2 * (site2.y() + site3.y());
+
+ // Create new circle event.
+ double c_x = (b1*dif_y2 - b2*dif_y1) / a;
+ double c_y = (b2*dif_x1 - b1*dif_x2) / a;
+ double radius = sqrt((c_x-site2.x())*(c_x-site2.x()) +
+ (c_y-site2.y())*(c_y-site2.y()));
+ c_event = make_circle_event<double>(c_x, c_y, c_x + radius);
+ return true;
+ }
- void reset() {
- while (!circle_events_.empty())
- circle_events_.pop();
- circle_events_list_.clear();
+ // Create circle event from two point sites and one segment site.
+ // TODO (asydorchuk): make_precision optimizations.
+ template <typename T>
+ static bool create_circle_event_pps(const site_event<T> &site1,
+ const site_event<T> &site2,
+ const site_event<T> &site3,
+ int segment_index,
+ circle_event<T> &c_event) {
+ if (segment_index != 2) {
+ if (orientation_test(site3.get_point0(), site1.get_point0(),
+ site2.get_point0()) != RIGHT_ORIENTATION &&
+ detail::orientation_test(site3.get_point1(), site1.get_point0(),
+ site2.get_point0()) != RIGHT_ORIENTATION)
+ return false;
+ } else {
+ if (orientation_test(site1.get_point0(), site3.get_point0(),
+ site2.get_point0()) != RIGHT_ORIENTATION &&
+ detail::orientation_test(site1.get_point0(), site3.get_point1(),
+ site2.get_point0()) != RIGHT_ORIENTATION)
+ return false;
}
- bool empty() {
- remove_not_active_events();
- return circle_events_.empty();
- }
+ double line_a = site3.get_point1().y() - site3.get_point0().y();
+ double line_b = site3.get_point0().x() - site3.get_point1().x();
+ double vec_x = site2.y() - site1.y();
+ double vec_y = site1.x() - site2.x();
+ double teta = line_b * vec_y + line_a * vec_x;
+ double A = line_a * (site1.x() - site3.get_point1().x()) +
+ line_b * (site1.y() - site3.get_point1().y());
+ double B = line_a * (site2.x() - site3.get_point1().x()) +
+ line_b * (site2.y() - site3.get_point1().y());
+ double denom = line_b * vec_x - line_a * vec_y;
+ double det = sqrt((teta * teta + denom * denom) * A * B);
+ double t;
+
+ if (orientation_test(static_cast<long long>(line_a),
+ static_cast<long long>(line_b),
+ static_cast<long long>(vec_x),
+ static_cast<long long>(vec_y)) == COLINEAR) {
+ t = (teta * teta - 4.0 * A * B) / (4.0 * teta * (A + B));
+ } else {
+ if (segment_index == 2)
+ det = -det;
+ t = (teta * (A + B) + 2.0 * det) / (2.0 * denom * denom);
+ }
+ double c_x = 0.5 * (site1.x() + site2.x()) + t * vec_x;
+ double c_y = 0.5 * (site1.y() + site2.y()) + t * vec_y;
+ double radius = sqrt((c_x-site2.x())*(c_x-site2.x()) +
+ (c_y-site2.y())*(c_y-site2.y()));
+ c_event = make_circle_event<double>(c_x, c_y, c_x + radius);
+ return true;
+ }
- const circle_event_type &top() {
- remove_not_active_events();
- return (*circle_events_.top());
+ // Create circle event from one point site and two segment sites.
+ // TODO (asydorchuk): make precision optimizations.
+ 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 (site1.get_site_index() == site2.get_site_index()) {
+ return false;
}
-
- void pop() {
- circle_event_type_it temp_it = circle_events_.top();
- circle_events_.pop();
- circle_events_list_.erase(temp_it);
+ 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);
+
+ kOrientation valid_orient1 = orientation_test(segm_end1, segm_start1, segm_start2);
+ kOrientation valid_orient2 = orientation_test(segm_end1, segm_start1, segm_end2);
+ if (valid_orient1 != RIGHT_ORIENTATION && valid_orient2 != RIGHT_ORIENTATION) {
+ return false;
}
- circle_event_type_it push(const circle_event_type &c_event) {
- circle_events_list_.push_front(c_event);
- circle_events_.push(circle_events_list_.begin());
- return circle_events_list_.begin();
+ double a1 = static_cast<double>(segm_end1.x() - segm_start1.x());
+ double b1 = static_cast<double>(segm_end1.y() - segm_start1.y());
+ double a2 = static_cast<double>(segm_end2.x() - segm_start2.x());
+ double b2 = static_cast<double>(segm_end2.y() - segm_start2.y());
+
+ kOrientation segm_orient = orientation_test(
+ static_cast<long long>(a1), static_cast<long long>(b1),
+ static_cast<long long>(a2), static_cast<long long>(b2));
+
+ if (segm_orient == COLINEAR) {
+ double a = a1 * a1 + b1 * b1;
+ double b = a1 * ((static_cast<double>(segm_start1.x()) +
+ static_cast<double>(segm_start2.x())) * 0.5 -
+ static_cast<double>(site1.x())) +
+ b1 * ((static_cast<double>(segm_start1.y()) +
+ static_cast<double>(segm_start2.y())) * 0.5 -
+ static_cast<double>(site1.y()));
+ double c = b1 * (static_cast<double>(segm_start2.x()) -
+ static_cast<double>(segm_start1.x())) -
+ a1 * (static_cast<double>(segm_start2.y()) -
+ static_cast<double>(segm_start1.y()));
+ double det = -(b1 * (static_cast<double>(site1.x()) -
+ static_cast<double>(segm_start1.x())) -
+ a1 * (static_cast<double>(site1.y()) -
+ static_cast<double>(segm_start1.y()))) *
+ (b1 * (static_cast<double>(site1.x()) -
+ static_cast<double>(segm_start2.x())) -
+ a1 * (static_cast<double>(site1.y()) -
+ static_cast<double>(segm_start2.y())));
+ double t = ((point_index == 2) ? (-b + sqrt(det)) : (-b - sqrt(det))) / a;
+ double c_x = 0.5 * (static_cast<double>(segm_start1.x() + segm_start2.x())) + a1 * t;
+ double c_y = 0.5 * (static_cast<double>(segm_start1.y() + segm_start2.y())) + b1 * t;
+ double radius = 0.5 * fabs(c) / sqrt(a);
+ c_event = make_circle_event<double>(c_x, c_y, c_x + radius);
+ return true;
+ } else {
+ double c1 = b1 * static_cast<double>(segm_end1.x()) -
+ a1 * static_cast<double>(segm_end1.y());
+ double c2 = a2 * static_cast<double>(segm_end2.y()) -
+ b2 * static_cast<double>(segm_end2.x());
+ double denom = (b1 * a2 - b2 * a1);
+ double intersection_x = (c1 * a2 + a1 * c2) / denom;
+ double intersection_y = (b1 * c2 + b2 * c1) / denom;
+ double sqr_sum1 = sqrt(a1 * a1 + b1 * b1);
+ double sqr_sum2 = sqrt(a2 * a2 + b2 * b2);
+ double vec_x = a1 * sqr_sum2 + a2 * sqr_sum1;
+ double vec_y = b1 * sqr_sum2 + b2 * sqr_sum1;
+ double dx = intersection_x - site1.get_point0().x();
+ double dy = intersection_y - site1.get_point0().y();
+ double a = a1 * a2 + b1 * b2 + sqr_sum1 * sqr_sum2;
+ double b = dx * vec_x + dy * vec_y;
+ double det = -2.0 * a * (b1 * dx - a1 * dy) * (b2 * dx - a2 * dy);
+ double t = ((point_index == 2) ? (-b + sqrt(det)) : (-b - sqrt(det))) / (a * a);
+ double c_x = intersection_x + vec_x * t;
+ double c_y = intersection_y + vec_y * t;
+ double radius = fabs(denom * t);
+ c_event = make_circle_event<double>(c_x, c_y, c_x + radius);
+ return true;
}
+ }
- private:
- struct comparison {
- bool operator() (const circle_event_type_it &node1,
- const circle_event_type_it &node2) const {
- return (*node1) > (*node2);
- }
- };
-
- void remove_not_active_events() {
- while (!circle_events_.empty() && !circle_events_.top()->is_active())
- pop();
+ // Create circle event from 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()) {
+ return false;
}
-
- std::priority_queue< circle_event_type_it,
- std::vector<circle_event_type_it>,
- comparison > circle_events_;
- std::list<circle_event_type> circle_events_list_;
-
- //Disallow copy constructor and operator=
- circle_events_queue(const circle_events_queue&);
- void operator=(const circle_events_queue&);
- };
+ double a1 = static_cast<double>(site1.x1(true) - site1.x0(true));
+ double b1 = static_cast<double>(site1.y1(true) - site1.y0(true));
+ double c1 = b1 * static_cast<double>(site1.x0(true)) -
+ a1 * static_cast<double>(site1.y0(true));
+ double a2 = static_cast<double>(site2.x1(true) - site2.x0(true));
+ double b2 = static_cast<double>(site2.y1(true) - site2.y0(true));
+ double c2 = b2 * static_cast<double>(site2.x0(true)) -
+ a2 * static_cast<double>(site2.y0(true));
+ double a3 = static_cast<double>(site3.x1(true) - site3.x0(true));
+ double b3 = static_cast<double>(site3.y1(true) - site3.y0(true));
+ double c3 = b3 * static_cast<double>(site3.x0(true)) -
+ a3 * static_cast<double>(site3.y0(true));
+ double sqr_sum1 = sqrt(a1 * a1 + b1 * b1);
+ double sqr_sum2 = sqrt(a2 * a2 + b2 * b2);
+ double sqr_sum3 = sqrt(a3 * a3 + b3 * b3);
+ double vec_a1 = -a1 * sqr_sum2 + a2 * sqr_sum1;
+ double vec_b1 = -b1 * sqr_sum2 + b2 * sqr_sum1;
+ double vec_c1 = -c1 * sqr_sum2 + c2 * sqr_sum1;
+ double vec_a2 = -a2 * sqr_sum3 + a3 * sqr_sum2;
+ double vec_b2 = -b2 * sqr_sum3 + b3 * sqr_sum2;
+ double vec_c2 = -c2 * sqr_sum3 + c3 * sqr_sum2;
+ double det = vec_a1 * vec_b2 - vec_b1 * vec_a2;
+ double c_x = (vec_a1 * vec_c2 - vec_c1 * vec_a2) / det;
+ double c_y = (vec_b1 * vec_c2 - vec_c1 * vec_b2) / det;
+ double radius = fabs(-b2 * c_x + a2 * c_y + c2) / sqr_sum2;
+ c_event = make_circle_event<double>(c_x, c_y, c_x + radius);
+ return true;
+ }
///////////////////////////////////////////////////////////////////////////
// VORONOI BEACH LINE TYPES ///////////////////////////////////////////////
@@ -936,19 +1177,27 @@
right_site_ = site;
}
+ void set_left_site_inverse() {
+ left_site_.set_inverse();
+ }
+
+ void set_right_site_inverse() {
+ right_site_.set_inverse();
+ }
+
coordinate_type get_comparison_x() const {
return (left_site_.get_site_index() >= right_site_.get_site_index()) ?
left_site_.x() : right_site_.x();
}
- std::pair<coordinate_type, int> get_comparison_y() const {
- // If node is lookup node.
- if (left_site_.get_site_index() == right_site_.get_site_index())
+ 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()) {
return std::make_pair(left_site_.y(), 0);
+ }
if (left_site_.get_site_index() > right_site_.get_site_index()) {
- if (left_site_.x() != right_site_.x() &&
- left_site_.is_segment() && left_site_.is_vertical())
- return std::make_pair(left_site_.get_point1().y(), 1);
+ if (!is_new_node && left_site_.is_segment() && left_site_.is_vertical()) {
+ return std::make_pair(left_site_.y1(), 1);
+ }
return std::make_pair(left_site_.y(), 1);
}
return std::make_pair(right_site_.y(), -1);
@@ -959,9 +1208,9 @@
left_site_.get_site_index() : right_site_.get_site_index();
}
- const Point2D &get_comparison_point() const {
+ const site_event_type &get_comparison_site() const {
return (left_site_.get_site_index() >= right_site_.get_site_index()) ?
- left_site_.get_point0() : right_site_.get_point0();
+ left_site_ : right_site_;
}
bool less(const Point2D &new_site) const {
@@ -988,19 +1237,15 @@
}
bool less_ps(const Point2D &new_site) const {
- return less_predicate(right_site_.get_point0(), right_site_.get_point1(),
- left_site_.get_point0(), new_site, false);
+ return less_predicate(left_site_.get_point0(), right_site_, new_site, false);
}
bool less_sp(const Point2D &new_site) const {
- return less_predicate(left_site_.get_point0(), left_site_.get_point1(),
- right_site_.get_point0(), new_site, true);
+ return less_predicate(right_site_.get_point0(), left_site_, new_site, true);
}
bool less_ss(const Point2D &new_site) const {
- return less_predicate(left_site_.get_point0(), left_site_.get_point1(),
- right_site_.get_point0(), right_site_.get_point1(),
- new_site, false);
+ return less_predicate(left_site_, right_site_, new_site);
}
private:
@@ -1053,16 +1298,27 @@
coordinate_type node2_x = node2.get_comparison_x();
if (node1_x < node2_x) {
- return node1.less(node2.get_comparison_point());
+ return node1.less(node2.get_comparison_site().get_point0());
} else if (node1_x > node2_x) {
- return !node2.less(node1.get_comparison_point());
+ return !node2.less(node1.get_comparison_site().get_point0());
} else {
- std::pair<coordinate_type, int> y1 = node1.get_comparison_y();
- std::pair<coordinate_type, int> y2 = node2.get_comparison_y();
- if (y1 != y2)
- return y1 < y2;
- return (node1.get_comparison_index() > node2.get_comparison_index()) ?
- (y1.second > 0) : (y2.second < 0);
+ 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();
+ if (y1.first != y2.first) {
+ return y1.first < y2.first;
+ }
+ return (!node1.get_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);
+ if (y1.first != y2.first) {
+ return y1.first < y2.first;
+ }
+ return (!node2.get_comparison_site().is_segment()) ? (y2.second > 0) : true;
+ }
}
}
};
Modified: sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_segment_builder.hpp
==============================================================================
--- sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_segment_builder.hpp (original)
+++ sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_segment_builder.hpp 2010-09-28 15:28:09 EDT (Tue, 28 Sep 2010)
@@ -47,12 +47,6 @@
// Clear all data structures.
reset();
- // Return if there are no sites.
- if (points.empty() && segments.empty()) {
- site_events_iterator_ = site_events_.begin();
- return;
- }
-
// TODO(asydorchuk): Add segments intersection check.
// Avoid additional memory reallocations.
@@ -90,14 +84,26 @@
// 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.
output_.init(site_events_.size());
+ }
- int skip = 0;
- if (site_events_.size() == 1) {
+ void reset() {
+ output_.reset();
+ circle_events_.reset();
+ beach_line_.clear();
+ site_events_.clear();
+ site_events_iterator_ = site_events_.begin();
+ }
+
+ void run_sweepline() {
+ // Init beach line.
+ if (site_events_.empty()) {
+ return;
+ } else if (site_events_.size() == 1) {
// Handle one input site case.
output_.process_single_site(site_events_[0]);
- skip = 1;
site_events_iterator_++;
} else {
+ int skip = 0;
// Init beach line.
while(site_events_iterator_ != site_events_.end() &&
site_events_iterator_->x() == site_events_.begin()->x() &&
@@ -114,29 +120,20 @@
init_beach_line_collinear_sites();
}
}
- }
-
- void reset() {
- output_.reset();
- circle_events_.reset();
- beach_line_.clear();
- site_events_.clear();
- site_events_iterator_ = site_events_.begin();
- }
- void run_sweepline() {
// Algorithm stops when there are no events to process.
while (!circle_events_.empty() ||
!(site_events_iterator_ == site_events_.end())) {
- if (circle_events_.empty())
+ if (circle_events_.empty()) {
process_site_event();
- else if (site_events_iterator_ == site_events_.end())
+ } else if (site_events_iterator_ == site_events_.end()) {
process_circle_event();
- else {
- if (circle_events_.top().compare(*site_events_iterator_) > 0)
+ } else {
+ if (circle_events_.top().compare(*site_events_iterator_) > 0) {
process_site_event();
- else
+ } else {
process_circle_event();
+ }
}
}
@@ -173,6 +170,7 @@
typedef typename detail::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;
void init_beach_line() {
// The first site is always a point.
@@ -181,17 +179,8 @@
site_events_iterator_type it_second = site_events_.begin();
it_second++;
- // The second site might be a point or a segment.
- // Create two new beach line nodes.
- Key new_left_node(*it_first, *it_second);
- Key new_right_node(*it_second, *it_first);
-
- // Update output.
- edge_type *edge = output_.insert_new_edge(*it_first, *it_second);
-
- // Insert two new nodes into the binary search tree.
- beach_line_.insert(std::pair<Key, Value>(new_left_node, Value(edge)));
- beach_line_.insert(std::pair<Key, Value>(new_right_node, Value(edge->twin)));
+ // Insert new nodes.
+ insert_new_arc(*it_first, *it_first, *it_second, beach_line_.begin());
// The second site has been already processed.
site_events_iterator_++;
@@ -222,61 +211,69 @@
// Uses special comparison function for the lower bound and insertion
// operations.
void process_site_event() {
- const site_event_type &site_event = *site_events_iterator_;
+ site_event_type site_event = *site_events_iterator_;
site_events_iterator_++;
+ // If new site is end point of some segment, remove temporary nodes from
+ // beach line data structure.
+ if (!site_event.is_segment()) {
+ while (!end_points_.empty() && end_points_.top().first == site_event.get_point0()) {
+ beach_line_.erase(end_points_.top().second);
+ end_points_.pop();
+ }
+ }
+
// 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);
- beach_line_iterator position = it;
+ int it_dist = site_event.is_segment() ? 2 : 1;
- site_event_type site_arc;
if (it == beach_line_.end()) {
it--;
- site_arc = it->first.get_right_site();
+ 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_event, position);
- new_node_it--;
+ beach_line_iterator new_node_it = insert_new_arc(site_arc, site_arc, site_event, it);
+ std::advance(new_node_it, -it_dist);
// Add candidate circle to the event queue.
- activate_circle_event(it->first.get_left_site(),
- it->first.get_right_site(),
- site_event,
- new_node_it);
+ activate_circle_event(it->first.get_left_site(), it->first.get_right_site(),
+ site_event, new_node_it);
} else if (it == beach_line_.begin()) {
- site_arc = it->first.get_left_site();
+ const site_event_type &site_arc = it->first.get_left_site();
// Insert new arc into the sweepline.
- beach_line_iterator new_node_it = insert_new_arc(site_arc, site_event, position);
- new_node_it++;
+ insert_new_arc(site_arc, site_arc, site_event, it);
// Add candidate circle to the event queue.
- activate_circle_event(site_event,
- it->first.get_left_site(),
- it->first.get_right_site(),
- new_node_it);
+ if (site_event.is_segment()) {
+ site_event.set_inverse();
+ }
+ activate_circle_event(site_event, it->first.get_left_site(),
+ it->first.get_right_site(), it);
} else {
- site_arc = it->first.get_left_site();
- const site_event_type &site2 = it->first.get_left_site();
+ const site_event_type &site_arc2 = it->first.get_left_site();
const site_event_type &site3 = it->first.get_right_site();
// Remove 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();
-
// Insert new arc into the sweepline.
- beach_line_iterator new_node_it = insert_new_arc(site_arc, site_event, position);
+ beach_line_iterator new_node_it =
+ insert_new_arc(site_arc1, site_arc2, site_event, it);
// Add candidate circles to the event queue.
- new_node_it--;
- activate_circle_event(site1, site2, site_event, new_node_it);
- new_node_it++;
- new_node_it++;
- activate_circle_event(site_event, site2, site3, new_node_it);
+ std::advance(new_node_it, -it_dist);
+ activate_circle_event(site1, site_arc1, site_event, new_node_it);
+ if (site_event.is_segment()) {
+ site_event.set_inverse();
+ }
+ std::advance(new_node_it, it_dist + 1);
+ activate_circle_event(site_event, site_arc2, site3, new_node_it);
}
}
@@ -312,6 +309,12 @@
// why we use const_cast there and take all the responsibility that
// map data structure keeps correct ordering.
const_cast<Key &>(it_first->first).set_right_site(site3);
+ if (site1.is_segment() && site1.get_point0(true) == site3.get_point0(true)) {
+ const_cast<Key &>(it_first->first).set_left_site_inverse();
+ }
+ if (site3.is_segment() && site3.get_point1(true) == site1.get_point0(true)) {
+ const_cast<Key &>(it_first->first).set_right_site_inverse();
+ }
it_first->second.edge = output_.insert_new_edge(site1, site2, site3, circle_event,
bisector1, bisector2);
beach_line_.erase(it_last);
@@ -342,145 +345,33 @@
}
// Insert new arc below site arc into the beach line.
- beach_line_iterator insert_new_arc(const site_event_type &site_arc,
+ 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_arc, site_event);
- Key new_right_node(site_event, site_arc);
+ Key new_left_node(site_arc1, site_event);
+ Key new_right_node(site_event, site_arc2);
+ if (site_event.is_segment()) {
+ new_right_node.set_left_site_inverse();
+ }
// Insert two new nodes into the binary search tree.
// Update output.
- edge_type *edge = output_.insert_new_edge(site_arc, site_event);
+ edge_type *edge = output_.insert_new_edge(site_arc2, site_event);
beach_line_iterator it = beach_line_.insert(position,
std::pair<Key, Value>(new_right_node, Value(edge->twin)));
- beach_line_.insert(it, std::pair<Key, Value>(new_left_node, Value(edge)));
+ if (site_event.is_segment()) {
+ Key new_node(site_event, site_event);
+ new_node.set_right_site_inverse();
+ 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));
+ }
+ beach_line_.insert(position, std::pair<Key, Value>(new_left_node, Value(edge)));
return it;
}
- // Create circle event from three point sites.
- bool create_circle_event_ppp(const site_event_type &site1,
- const site_event_type &site2,
- const site_event_type &site3,
- circle_event_type &c_event) const {
- // Check if bisectors intersect.
- if (detail::orientation_test(site1.get_point0(),site2.get_point0(),
- site3.get_point0()) != detail::RIGHT_ORIENTATION)
- return false;
-
- coordinate_type a = ((site1.x() - site2.x()) * (site2.y() - site3.y()) -
- (site1.y() - site2.y()) * (site2.x() - site3.x())) *
- static_cast<coordinate_type>(2.0);
-
- coordinate_type b1 = (site1.x() - site2.x()) * (site1.x() + site2.x()) +
- (site1.y() - site2.y()) * (site1.y() + site2.y());
- coordinate_type b2 = (site2.x() - site3.x()) * (site2.x() + site3.x()) +
- (site2.y() - site3.y()) * (site2.y() + site3.y());
-
- // Create new circle event.
- coordinate_type c_x = (b1*(site2.y() - site3.y()) - b2*(site1.y() - site2.y())) / a;
- coordinate_type c_y = (b2*(site1.x() - site2.x()) - b1*(site2.x() - site3.x())) / a;
- coordinate_type radius = sqrt((c_x-site2.x())*(c_x-site2.x()) +
- (c_y-site2.y())*(c_y-site2.y()));
- c_event = detail::make_circle_event<coordinate_type>(c_x, c_y, c_x + radius);
- return true;
- }
-
- // Create circle event from two point sites and one segment site.
- bool create_circle_event_pps(const site_event_type &site1,
- const site_event_type &site2,
- const site_event_type &site3,
- int segment_index,
- circle_event_type &c_event) const {
- // Check if bisectors intersect.
- detail::kOrientation orientation1 = detail::orientation_test(
- site3.get_point0(), site3.get_point1(), site1.get_point0());
- detail::kOrientation orientation2 = detail::orientation_test(
- site3.get_point0(), site3.get_point1(), site2.get_point0());
-
- // If point sites are situated on different sides of segment return false.
- if (static_cast<int>(orientation1) * static_cast<int>(orientation2) == -1)
- return false;
-
- if ((orientation1 == detail::COLINEAR) && (orientation2 == detail::COLINEAR))
- return false;
-
- bool right_oriented = (orientation1 == detail::RIGHT_ORIENTATION) ||
- (orientation2 == detail::RIGHT_ORIENTATION);
-
- // Additional orientation test if segment index is equal to 1 or 3.
- if (segment_index != 2) {
- const Point2D &test_point = ((segment_index == 1) ^ right_oriented) ?
- site3.get_point1() : site3.get_point0();
- if ((segment_index == 1 &&
- detail::orientation_test(test_point, site1.get_point0(), site2.get_point0()) !=
- detail::RIGHT_ORIENTATION) ||
- (segment_index == 3 &&
- detail::orientation_test(site1.get_point0(), site2.get_point0(), test_point) !=
- detail::RIGHT_ORIENTATION))
- return false;
- } else {
- if (detail::orientation_test(site1.get_point0(), site3.get_point0(),
- site2.get_point0()) != detail::RIGHT_ORIENTATION &&
- detail::orientation_test(site1.get_point0(), site3.get_point1(),
- site2.get_point0()) != detail::RIGHT_ORIENTATION)
- 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 vec_x = site2.y() - site1.y();
- double vec_y = site1.x() - site2.x();
- double teta = line_b * vec_y + line_a * vec_x;
- double A = line_a * (site1.x() - site3.get_point1().x()) +
- line_b * (site1.y() - site3.get_point1().y());
- double B = line_a * (site2.x() - site3.get_point1().x()) +
- line_b * (site2.y() - site3.get_point1().y());
- double denom = line_b * vec_x - line_a * vec_y;
- double det = 0.0;
- if (orientation1 != detail::COLINEAR && orientation2 != detail::COLINEAR)
- det = sqrt((teta * teta + denom * denom) * A * B);
- double t;
-
- if (detail::orientation_test(static_cast<long long>(line_a),
- static_cast<long long>(line_b),
- static_cast<long long>(vec_x),
- static_cast<long long>(vec_y)) == detail::COLINEAR) {
- t = (teta * teta - 4.0 * A * B) / (4.0 * teta * (A + B));
- } else {
- if (segment_index == 2)
- det = -det;
- t = (teta * (A + B) + 2.0 * det) / (2.0 * denom * denom);
- }
- double c_x = 0.5 * (site1.x() + site2.x()) + t * vec_x;
- double c_y = 0.5 * (site1.y() + site2.y()) + t * vec_y;
- double radius = sqrt((c_x-site2.x())*(c_x-site2.x()) +
- (c_y-site2.y())*(c_y-site2.y()));
- c_event = detail::make_circle_event<coordinate_type>(c_x, c_y, c_x + radius);
- return true;
- }
-
- // Create circle event from one point site and two segment sites.
- bool create_circle_event_pss(const site_event_type &site1,
- const site_event_type &site2,
- const site_event_type &site3,
- int point_index,
- circle_event_type &c_event) const {
- // Check if bisectors intersect.
- // Not implemented yet.
- return false;
- }
-
- // Create circle event from three segment sites.
- bool create_circle_event_sss(const site_event_type &site1,
- const site_event_type &site2,
- const site_event_type &site3,
- circle_event_type &c_event) const {
- // Check if bisectors intersecto.
- // Not implemented yet.
- return false;
- }
-
// Create circle event from the given three points.
bool create_circle_event(const site_event_type &site1,
const site_event_type &site2,
@@ -489,29 +380,29 @@
if (!site1.is_segment()) {
if (!site2.is_segment()) {
if (!site3.is_segment()) {
- return create_circle_event_ppp(site1, site2, site3, c_event);
+ return detail::create_circle_event_ppp(site1, site2, site3, c_event);
} else {
- return create_circle_event_pps(site1, site2, site3, 3, c_event);
+ return detail::create_circle_event_pps(site1, site2, site3, 3, c_event);
}
} else {
if (!site3.is_segment()) {
- return create_circle_event_pps(site1, site3, site2, 2, c_event);
+ return detail::create_circle_event_pps(site1, site3, site2, 2, c_event);
} else {
- return create_circle_event_pss(site1, site2, site3, 1, c_event);
+ return detail::create_circle_event_pss(site1, site2, site3, 1, c_event);
}
}
} else {
if (!site2.is_segment()) {
if (!site3.is_segment()) {
- return create_circle_event_pps(site2, site3, site1, 1, c_event);
+ return detail::create_circle_event_pps(site2, site3, site1, 1, c_event);
} else {
- return create_circle_event_pss(site2, site1, site3, 2, c_event);
+ return detail::create_circle_event_pss(site2, site1, site3, 2, c_event);
}
} else {
if (!site3.is_segment()) {
- return create_circle_event_pss(site3, site1, site2, 3, c_event);
+ return detail::create_circle_event_pss(site3, site1, site2, 3, c_event);
} else {
- return create_circle_event_sss(site1, site2, site3, c_event);
+ return detail::create_circle_event_sss(site1, site2, site3, c_event);
}
}
}
@@ -530,8 +421,16 @@
}
private:
+ struct end_point_comparison {
+ bool operator() (const end_point_type &end1, const end_point_type &end2) const {
+ return end1.first > end2.first;
+ }
+ };
+
std::vector<site_event_type> site_events_;
site_events_iterator_type site_events_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_;
Modified: sandbox/SOC/2010/sweepline/libs/sweepline/test/Jamfile.v2
==============================================================================
--- sandbox/SOC/2010/sweepline/libs/sweepline/test/Jamfile.v2 (original)
+++ sandbox/SOC/2010/sweepline/libs/sweepline/test/Jamfile.v2 2010-09-28 15:28:09 EDT (Tue, 28 Sep 2010)
@@ -27,6 +27,7 @@
test-suite "sweepline_segment"
:
+ [ run voronoi_segment/segment_circle_event_creation_test.cpp ]
[ run voronoi_segment/segment_event_queue_test.cpp ]
[ run voronoi_segment/segment_event_types_test.cpp ]
[ run voronoi_segment/segment_node_comparer_test.cpp ]
Added: sandbox/SOC/2010/sweepline/libs/sweepline/test/voronoi_segment/segment_circle_event_creation_test.cpp
==============================================================================
--- (empty file)
+++ sandbox/SOC/2010/sweepline/libs/sweepline/test/voronoi_segment/segment_circle_event_creation_test.cpp 2010-09-28 15:28:09 EDT (Tue, 28 Sep 2010)
@@ -0,0 +1,191 @@
+// Boost sweepline library segment_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.
+
+#include "../test_type_list.hpp"
+#include "boost/sweepline/voronoi_segment_sweepline.hpp"
+using namespace boost::sweepline;
+
+#define BOOST_TEST_MODULE circle_event_creation
+#include <boost/test/test_case_template.hpp>
+
+#define CHECK_CIRCLE_EVENT(c_e, c_x, c_y, l_x) \
+ BOOST_CHECK_EQUAL(c_e.get_center().x() == static_cast<T>(c_x), true); \
+ BOOST_CHECK_EQUAL(c_e.get_center().y() == static_cast<T>(c_y), true); \
+ BOOST_CHECK_EQUAL(c_e.get_lower_x() == static_cast<T>(l_x), true)
+
+// Test for (point, point, point) case.
+BOOST_AUTO_TEST_CASE_TEMPLATE(circle_event_creation_ppp_test1, T, test_types) {
+ typedef typename detail::site_event<T> site_event_type;
+ site_event_type site1(static_cast<T>(0), static_cast<T>(0), 0);
+ site_event_type site2(static_cast<T>(-8), static_cast<T>(0), 1);
+ site_event_type site3(static_cast<T>(0), static_cast<T>(6), 2);
+ typename detail::circle_event<T> c_event;
+ bool is_event = detail::create_circle_event_ppp(site1, site2, site3, c_event);
+ BOOST_CHECK_EQUAL(is_event, true);
+ CHECK_CIRCLE_EVENT(c_event, -4, 3, 1);
+ is_event = detail::create_circle_event_ppp(site2, site1, site3, c_event);
+ BOOST_CHECK_EQUAL(is_event, false);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(circle_event_creation_ppp_test2, T, test_types) {
+ typedef typename detail::site_event<T> site_event_type;
+ int min_int = std::numeric_limits<int>::min();
+ int max_int = std::numeric_limits<int>::max();
+ site_event_type site1(static_cast<T>(min_int), static_cast<T>(min_int), 0);
+ 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;
+ 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);
+ BOOST_CHECK_EQUAL(is_event, false);
+}
+
+// Test for (point, point, segment) case.
+BOOST_AUTO_TEST_CASE_TEMPLATE(circle_event_creation_pps_test1, T, test_types) {
+ 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);
+ 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);
+ site_event_type site1_1(static_cast<T>(-2), static_cast<T>(0), 0);
+ 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);
+ 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);
+ BOOST_CHECK_EQUAL(is_event, false);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(circle_event_creation_pps_test2, T, test_types) {
+ 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);
+ 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);
+ bool is_event = detail::create_circle_event_pps(site1_1, site1_2, segm_site, 1, c_event);
+ BOOST_CHECK_EQUAL(is_event, true);
+ CHECK_CIRCLE_EVENT(c_event, 1, 6, 6);
+ is_event = detail::create_circle_event_pps(site1_2, site1_1, segm_site, 3, c_event);
+ BOOST_CHECK_EQUAL(is_event, true);
+ CHECK_CIRCLE_EVENT(c_event, 1, 14, 6);
+}
+
+// Test for (point, segment, segment) case.
+BOOST_AUTO_TEST_CASE_TEMPLATE(circle_event_creation_pss_test1, T, test_types) {
+ typedef typename detail::site_event<T> site_event_type;
+ point_2d<T> segm1_start(static_cast<T>(0), 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();
+ point_2d<T> segm2_start(static_cast<T>(-3), static_cast<T>(4));
+ point_2d<T> segm2_end(static_cast<T>(9), static_cast<T>(4));
+ site_event_type segm_site2(segm2_start, segm2_end, 1);
+ typename detail::circle_event<T> c_event;
+
+ site_event_type site1(static_cast<T>(5), static_cast<T>(2), 2);
+ bool is_event = detail::create_circle_event_pss(site1, segm_site1, segm_site2, 3, c_event);
+ BOOST_CHECK_EQUAL(is_event, true);
+ CHECK_CIRCLE_EVENT(c_event, 3, 2, 5);
+
+ site_event_type site2(static_cast<T>(0), static_cast<T>(0), 2);
+ is_event = detail::create_circle_event_pss(site2, segm_site1, segm_site2, 2, c_event);
+ BOOST_CHECK_EQUAL(is_event, true);
+ CHECK_CIRCLE_EVENT(c_event, 0, 2, 2);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(circle_event_creation_pss_test2, T, test_types) {
+ typedef typename detail::site_event<T> site_event_type;
+ 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();
+ 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);
+ typename detail::circle_event<T> c_event;
+ site_event_type site1(static_cast<T>(1), static_cast<T>(1), 2);
+ bool is_event = detail::create_circle_event_pss(site1, segm_site1, segm_site2, 2, c_event);
+ BOOST_CHECK_EQUAL(is_event, true);
+ CHECK_CIRCLE_EVENT(c_event, 6, 1, 11);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(circle_event_creation_pss_test3, T, test_types) {
+ typedef typename detail::site_event<T> site_event_type;
+ point_2d<T> segm1_start(static_cast<T>(0), static_cast<T>(0));
+ point_2d<T> segm1_end(static_cast<T>(5), static_cast<T>(0));
+ site_event_type segm_site1(segm1_start, segm1_end, 0);
+ segm_site1.set_inverse();
+ point_2d<T> segm2_start(static_cast<T>(-7), static_cast<T>(4));
+ point_2d<T> segm2_end(static_cast<T>(-1), static_cast<T>(12));
+ site_event_type segm_site2(segm2_start, segm2_end, 1);
+ typename detail::circle_event<T> c_event;
+ site_event_type site1(static_cast<T>(0), static_cast<T>(0), 2);
+ bool is_event = detail::create_circle_event_pss(site1, segm_site1, segm_site2, 2, c_event);
+ BOOST_CHECK_EQUAL(is_event, true);
+ CHECK_CIRCLE_EVENT(c_event, 0, 5, 5);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(circle_event_creation_pss_test4, T, test_types) {
+ typedef typename detail::site_event<T> site_event_type;
+ point_2d<T> point1(static_cast<T>(0), static_cast<T>(0));
+ point_2d<T> point2(static_cast<T>(4), static_cast<T>(0));
+ point_2d<T> point3(static_cast<T>(7), static_cast<T>(6));
+ point_2d<T> point4(static_cast<T>(-1), static_cast<T>(12));
+ site_event_type segm_site1(point1, point2, 0);
+ segm_site1.set_inverse();
+ site_event_type segm_site2(point3, point4, 1);
+ site_event_type point_site(point1, 2);
+ typename detail::circle_event<T> c_event;
+ bool is_event = detail::create_circle_event_pss(
+ point_site, segm_site1, segm_site2, 2, c_event);
+ BOOST_CHECK_EQUAL(is_event, true);
+ CHECK_CIRCLE_EVENT(c_event, 0, 5, 5);
+}
+
+// Test for (segment, segment, segment) case.
+BOOST_AUTO_TEST_CASE_TEMPLATE(circle_event_creation_sss_test1, T, test_types) {
+ typedef typename detail::site_event<T> site_event_type;
+ point_2d<T> point1(static_cast<T>(4), static_cast<T>(0));
+ point_2d<T> point2(static_cast<T>(0), static_cast<T>(0));
+ 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();
+ site_event_type segm_site2(point2, point3, 1);
+ site_event_type segm_site3(point3, point4, 2);
+ typename detail::circle_event<T> c_event;
+ bool is_event = detail::create_circle_event_sss(segm_site1, segm_site2, segm_site3, c_event);
+ BOOST_CHECK_EQUAL(is_event, true);
+ CHECK_CIRCLE_EVENT(c_event, 2, 2, 4);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(circle_event_creation_sss_test2, T, test_types) {
+ typedef typename detail::site_event<T> site_event_type;
+ point_2d<T> point1(static_cast<T>(40), static_cast<T>(30));
+ point_2d<T> point2(static_cast<T>(0), static_cast<T>(0));
+ point_2d<T> point3(static_cast<T>(-40), static_cast<T>(30));
+ point_2d<T> point4(static_cast<T>(0), static_cast<T>(60));
+ site_event_type segm_site1(point1, point2, 0);
+ segm_site1.set_inverse();
+ site_event_type segm_site2(point3, point4, 1);
+ site_event_type segm_site3(point4, point1, 2);
+ typename detail::circle_event<T> c_event;
+ bool is_event = detail::create_circle_event_sss(segm_site1, segm_site2, segm_site3, c_event);
+ BOOST_CHECK_EQUAL(is_event, true);
+ CHECK_CIRCLE_EVENT(c_event, 0, 30, 24);
+}
Modified: sandbox/SOC/2010/sweepline/libs/sweepline/test/voronoi_segment/segment_event_types_test.cpp
==============================================================================
--- sandbox/SOC/2010/sweepline/libs/sweepline/test/voronoi_segment/segment_event_types_test.cpp (original)
+++ sandbox/SOC/2010/sweepline/libs/sweepline/test/voronoi_segment/segment_event_types_test.cpp 2010-09-28 15:28:09 EDT (Tue, 28 Sep 2010)
@@ -150,13 +150,13 @@
point4 = make_point_2d<T>(static_cast<T>(10), static_cast<T>(9));
site2 = make_site_event<T>(point2, point4, 0);
- EVENT_TYPES_CHECK_COMPARISON(site1, site2, arr1_1);
- EVENT_TYPES_CHECK_COMPARISON(site2, site1, arr1_2);
+ EVENT_TYPES_CHECK_COMPARISON(site1, site2, arr1_2);
+ EVENT_TYPES_CHECK_COMPARISON(site2, site1, arr1_1);
point4 = make_point_2d<T>(static_cast<T>(9), static_cast<T>(10));
site2 = make_site_event<T>(point2, point4, 0);
- EVENT_TYPES_CHECK_COMPARISON(site1, site2, arr1_2);
- EVENT_TYPES_CHECK_COMPARISON(site2, site1, arr1_1);
+ EVENT_TYPES_CHECK_COMPARISON(site1, site2, arr1_1);
+ EVENT_TYPES_CHECK_COMPARISON(site2, site1, arr1_2);
}
BOOST_AUTO_TEST_CASE_TEMPLATE(circle_event_test1, T, test_types) {
@@ -211,4 +211,4 @@
site = make_site_event<T>(static_cast<T>(4), static_cast<T>(2), 0);
BOOST_CHECK_EQUAL(circle.compare(site), -1);
-}
\ No newline at end of file
+}
Modified: sandbox/SOC/2010/sweepline/libs/sweepline/test/voronoi_segment/segment_predicates_test.cpp
==============================================================================
--- sandbox/SOC/2010/sweepline/libs/sweepline/test/voronoi_segment/segment_predicates_test.cpp (original)
+++ sandbox/SOC/2010/sweepline/libs/sweepline/test/voronoi_segment/segment_predicates_test.cpp 2010-09-28 15:28:09 EDT (Tue, 28 Sep 2010)
@@ -14,27 +14,27 @@
#define BOOST_TEST_MODULE predicates_test
#include <boost/test/test_case_template.hpp>
-#define CHECK_ORIENTATION_EQUAL(p1, p2, p3, orientation) \
- BOOST_CHECK_EQUAL(detail::orientation_test(p1, p2, p3) == orientation, true)
+#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(p1, p2, p3, exp) \
- BOOST_CHECK_EQUAL(detail::fast_less_predicate(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_LESS_PREDICATE_PP(p1, p2, p3, exp) \
- BOOST_CHECK_EQUAL(detail::less_predicate(p1, p2, p3) == exp, true)
+#define CHECK_LESS_PREDICATE_PP(lp, rp, np, exp) \
+ BOOST_CHECK_EQUAL(detail::less_predicate(lp, rp, np) == exp, true)
-#define CHECK_LESS_PREDICATE_CHECK_PS(p1, p2, p3, p4, r, exp) \
- BOOST_CHECK_EQUAL(detail::less_predicate_check(p1, p2, p3, p4, r) == exp, true); \
+#define CHECK_FAST_LESS_PREDICATE_PS(lp, rs, np, reverse, exp) \
+ BOOST_CHECK_EQUAL(detail::fast_less_predicate(lp, rs, np, reverse) == exp, true); \
if (exp != detail::UNDEFINED) { \
bool exp_res = (exp == detail::LESS) ? true : false; \
- BOOST_CHECK_EQUAL(detail::less_predicate(p1, p2, p3, p4, r), exp_res); \
+ BOOST_CHECK_EQUAL(detail::less_predicate(lp, rs, np, reverse), exp_res); \
}
-#define CHECK_LESS_PREDICATE_PS(p1, p2, p3, p4, r, exp) \
- BOOST_CHECK_EQUAL(detail::less_predicate(p1, p2, p3, p4, r) == exp, true)
+#define CHECK_LESS_PREDICATE_PS(lp, rs, np, reverse, exp) \
+ BOOST_CHECK_EQUAL(detail::less_predicate(lp, rs, np, reverse) == exp, true)
-#define CHECK_LESS_PREDICATE_SS(p1, p2, p3, p4, p5, r, exp) \
- BOOST_CHECK_EQUAL(detail::less_predicate(p1, p2, p3, p4, p5, r) == exp, true)
+#define CHECK_LESS_PREDICATE_SS(ls, rs, np, exp) \
+ BOOST_CHECK_EQUAL(detail::less_predicate(ls, rs, np) == exp, true)
// This test uses integer values in the range [-2^31, 2^31), to validate
// orientation predicate for the hole integer range input coordinates.
@@ -70,7 +70,7 @@
}
// Test fast point-point predicate.
-BOOST_AUTO_TEST_CASE_TEMPLATE(fast_less_predicate_test1, T, test_types) {
+BOOST_AUTO_TEST_CASE_TEMPLATE(fast_less_predicate_point_point_test1, T, test_types) {
point_2d<T> point1 = make_point_2d<T>(static_cast<T>(-5), static_cast<T>(0));
point_2d<T> point2 = make_point_2d<T>(static_cast<T>(-8), static_cast<T>(9));
point_2d<T> point3 = make_point_2d<T>(static_cast<T>(-2), static_cast<T>(1));
@@ -108,188 +108,178 @@
}
// Vertical segment case.
-BOOST_AUTO_TEST_CASE_TEMPLATE(less_predicate_check_point_segment_test1, T, test_types) {
+BOOST_AUTO_TEST_CASE_TEMPLATE(fast_less_predicate_point_segment_test1, T, test_types) {
point_2d<T> segm_start = make_point_2d<T>(static_cast<T>(-4), static_cast<T>(0));
point_2d<T> segm_end = make_point_2d<T>(static_cast<T>(-4), static_cast<T>(20));
+ typename detail::site_event<T> segm_site(segm_start, segm_end, 0);
point_2d<T> site_p = make_point_2d<T>(static_cast<T>(-2), static_cast<T>(10));
point_2d<T> new_p1 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(11));
point_2d<T> new_p2 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(9));
- CHECK_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p, new_p1,
- false, detail::UNDEFINED);
- CHECK_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p, new_p2,
- false, detail::MORE);
- CHECK_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p, new_p1,
- true, detail::LESS);
- CHECK_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p, new_p2,
- true, detail::UNDEFINED);
+ CHECK_FAST_LESS_PREDICATE_PS(site_p, segm_site, new_p1, false, detail::UNDEFINED);
+ CHECK_FAST_LESS_PREDICATE_PS(site_p, segm_site, new_p2, false, detail::MORE);
+ CHECK_FAST_LESS_PREDICATE_PS(site_p, segm_site, new_p1, true, detail::LESS);
+ CHECK_FAST_LESS_PREDICATE_PS(site_p, segm_site, new_p2, true, detail::UNDEFINED);
}
// Not vertical segment case. Site is to the left of the segment vector.
-BOOST_AUTO_TEST_CASE_TEMPLATE(less_predicate_check_point_segment_test2, T, test_types) {
+BOOST_AUTO_TEST_CASE_TEMPLATE(fast_less_predicate_point_segment_test2, T, test_types) {
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));
- CHECK_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p1, new_p1,
- false, detail::MORE);
- CHECK_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p1, new_p1,
- true, detail::MORE);
+ segm_site.set_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_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p1, new_p2,
- false, detail::MORE);
- CHECK_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p1, new_p2,
- true, detail::UNDEFINED);
+ 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);
point_2d<T> new_p3 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(4));
- CHECK_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p1, new_p3,
- false, detail::MORE);
- CHECK_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p1, new_p3,
- true, detail::UNDEFINED);
+ CHECK_FAST_LESS_PREDICATE_PS(site_p1, segm_site, new_p3, false, detail::MORE);
+ CHECK_FAST_LESS_PREDICATE_PS(site_p1, segm_site, new_p3, true, detail::UNDEFINED);
point_2d<T> new_p4 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(5));
- CHECK_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p1, new_p4,
- false, detail::UNDEFINED);
- CHECK_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p1, new_p4,
- true, detail::LESS);
+ CHECK_FAST_LESS_PREDICATE_PS(site_p1, segm_site, new_p4, false, detail::UNDEFINED);
+ CHECK_FAST_LESS_PREDICATE_PS(site_p1, segm_site, new_p4, true, detail::LESS);
}
// Not vertical segment case. Site is to the right of the segment vector.
-BOOST_AUTO_TEST_CASE_TEMPLATE(less_predicate_check_point_segment_test3, T, test_types) {
+BOOST_AUTO_TEST_CASE_TEMPLATE(fast_less_predicate_point_segment_test3, T, test_types) {
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));
point_2d<T> new_p1 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(1));
- CHECK_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p1, new_p1,
- false, detail::LESS);
- CHECK_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p1, new_p1,
- true, detail::LESS);
- CHECK_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p2, new_p1,
- false, detail::LESS);
- CHECK_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p2, new_p1,
- true, detail::LESS);
+ CHECK_FAST_LESS_PREDICATE_PS(site_p1, segm_site, new_p1, false, detail::LESS);
+ 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_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p1, new_p2,
- false, detail::UNDEFINED);
- CHECK_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p1, new_p2,
- true, detail::LESS);
- CHECK_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p2, new_p2,
- false, detail::UNDEFINED);
- CHECK_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p2, new_p2,
- true, detail::LESS);
+ 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);
+ CHECK_FAST_LESS_PREDICATE_PS(site_p2, segm_site, new_p2, false, detail::UNDEFINED);
+ CHECK_FAST_LESS_PREDICATE_PS(site_p2, segm_site, new_p2, true, detail::LESS);
point_2d<T> new_p3 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(-8));
- CHECK_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p1, new_p3,
- false, detail::UNDEFINED);
- CHECK_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p1, new_p3,
- true, detail::LESS);
- CHECK_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p2, new_p3,
- false, detail::UNDEFINED);
- CHECK_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p2, new_p3,
- true, detail::LESS);
+ CHECK_FAST_LESS_PREDICATE_PS(site_p1, segm_site, new_p3, false, detail::UNDEFINED);
+ CHECK_FAST_LESS_PREDICATE_PS(site_p1, segm_site, new_p3, true, detail::LESS);
+ CHECK_FAST_LESS_PREDICATE_PS(site_p2, segm_site, new_p3, false, detail::UNDEFINED);
+ CHECK_FAST_LESS_PREDICATE_PS(site_p2, segm_site, new_p3, true, detail::LESS);
point_2d<T> new_p4 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(-9));
- CHECK_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p1, new_p4,
- false, detail::MORE);
- CHECK_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p1, new_p4,
- true, detail::UNDEFINED);
- CHECK_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p2, new_p4,
- false, detail::MORE);
- CHECK_LESS_PREDICATE_CHECK_PS(segm_start, segm_end, site_p2, new_p4,
- true, detail::UNDEFINED);
+ CHECK_FAST_LESS_PREDICATE_PS(site_p1, segm_site, new_p4, false, detail::MORE);
+ CHECK_FAST_LESS_PREDICATE_PS(site_p1, segm_site, new_p4, true, detail::UNDEFINED);
+ CHECK_FAST_LESS_PREDICATE_PS(site_p2, segm_site, new_p4, false, detail::MORE);
+ CHECK_FAST_LESS_PREDICATE_PS(site_p2, segm_site, new_p4, true, detail::UNDEFINED);
}
// Test main point-segment predicate.
BOOST_AUTO_TEST_CASE_TEMPLATE(less_predicate_point_segment_test1, T, test_types) {
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_site1(segm_start, segm_end, 0);
+ typename detail::site_event<T> segm_site2(segm_start, segm_end, 0);
+ segm_site2.set_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));
point_2d<T> new_p1 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(1));
- CHECK_LESS_PREDICATE_PS(segm_start, segm_end, site_p1, new_p1, true, false);
+ CHECK_LESS_PREDICATE_PS(site_p1, segm_site2, new_p1, true, false);
point_2d<T> new_p2 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(4));
- CHECK_LESS_PREDICATE_PS(segm_start, segm_end, site_p1, new_p2, true, true);
+ CHECK_LESS_PREDICATE_PS(site_p1, segm_site2, new_p2, true, true);
point_2d<T> new_p3 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(5));
- CHECK_LESS_PREDICATE_PS(segm_start, segm_end, site_p1, new_p3, false, false);
+ CHECK_LESS_PREDICATE_PS(site_p1, segm_site2, new_p3, false, false);
point_2d<T> new_p4 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(7));
- CHECK_LESS_PREDICATE_PS(segm_start, segm_end, site_p1, new_p4, false, true);
+ CHECK_LESS_PREDICATE_PS(site_p1, segm_site2, new_p4, false, true);
point_2d<T> new_p5 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(-2));
- CHECK_LESS_PREDICATE_PS(segm_start, segm_end, site_p2, new_p5, false, false);
- CHECK_LESS_PREDICATE_PS(segm_start, segm_end, site_p3, new_p5, false, false);
+ CHECK_LESS_PREDICATE_PS(site_p2, segm_site1, new_p5, false, false);
+ CHECK_LESS_PREDICATE_PS(site_p3, segm_site1, new_p5, false, false);
point_2d<T> new_p6 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(-8));
- CHECK_LESS_PREDICATE_PS(segm_start, segm_end, site_p2, new_p6, false, false);
- CHECK_LESS_PREDICATE_PS(segm_start, segm_end, site_p3, new_p6, false, false);
+ CHECK_LESS_PREDICATE_PS(site_p2, segm_site1, new_p6, false, false);
+ CHECK_LESS_PREDICATE_PS(site_p3, segm_site1, new_p6, false, false);
point_2d<T> new_p7 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(-9));
- CHECK_LESS_PREDICATE_PS(segm_start, segm_end, site_p2, new_p7, true, true);
- CHECK_LESS_PREDICATE_PS(segm_start, segm_end, site_p3, new_p7, true, true);
+ CHECK_LESS_PREDICATE_PS(site_p2, segm_site1, new_p7, true, true);
+ CHECK_LESS_PREDICATE_PS(site_p3, segm_site1, new_p7, true, true);
point_2d<T> new_p8 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(-18));
- CHECK_LESS_PREDICATE_PS(segm_start, segm_end, site_p2, new_p8, true, false);
- CHECK_LESS_PREDICATE_PS(segm_start, segm_end, site_p3, new_p8, true, false);
+ CHECK_LESS_PREDICATE_PS(site_p2, segm_site1, new_p8, true, false);
+ CHECK_LESS_PREDICATE_PS(site_p3, segm_site1, new_p8, true, false);
point_2d<T> new_p9 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(-1));
- CHECK_LESS_PREDICATE_PS(segm_start, segm_end, site_p2, new_p9, false, true);
- CHECK_LESS_PREDICATE_PS(segm_start, segm_end, site_p3, new_p9, false, true);
+ CHECK_LESS_PREDICATE_PS(site_p2, segm_site1, new_p9, false, true);
+ CHECK_LESS_PREDICATE_PS(site_p3, segm_site1, new_p9, false, true);
}
// Test main segment-segment predicate.
BOOST_AUTO_TEST_CASE_TEMPLATE(less_predicate_segment_segment_test1, T, test_types) {
+ point_2d<T> segm_start1 = make_point_2d<T>(static_cast<T>(-5), static_cast<T>(0));
+ 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();
+
+ // New sites.
+ point_2d<T> new_site1 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(5));
+ point_2d<T> new_site2 = make_point_2d<T>(static_cast<T>(1), static_cast<T>(5));
+ point_2d<T> new_site3 = make_point_2d<T>(static_cast<T>(-1), static_cast<T>(5));
+
+ CHECK_LESS_PREDICATE_SS(segm_site1_1, segm_site1_2, new_site1, false);
+ CHECK_LESS_PREDICATE_SS(segm_site1_1, segm_site1_2, new_site2, true);
+ CHECK_LESS_PREDICATE_SS(segm_site1_1, segm_site1_2, new_site3, false);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(less_predicate_segment_segment_test2, T, test_types) {
+ typedef typename detail::site_event<T> site_event_type;
+
point_2d<T> segm_start1 = make_point_2d<T>(static_cast<T>(-5), static_cast<T>(5));
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();
- // Common start point case.
- point_2d<T> segm_start2_1 = make_point_2d<T>(static_cast<T>(-5), static_cast<T>(5));
- point_2d<T> segm_end2_1 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(6));
+ // New sites.
point_2d<T> new_site1 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(2));
- CHECK_LESS_PREDICATE_SS(segm_start1, segm_end1, segm_start2_1, segm_end2_1,
- new_site1, false, false);
-
- point_2d<T> new_site2 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(9));
- CHECK_LESS_PREDICATE_SS(segm_start1, segm_end1, segm_start2_1, segm_end2_1,
- new_site2, false, true);
-
- point_2d<T> new_site3 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(10));
- CHECK_LESS_PREDICATE_SS(segm_start1, segm_end1, segm_start2_1, segm_end2_1,
- new_site3, false, true);
-
- // No common segment end points.
- point_2d<T> segm_start2_2 = make_point_2d<T>(static_cast<T>(-2), static_cast<T>(4));
- point_2d<T> segm_end2_2 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(4));
- CHECK_LESS_PREDICATE_SS(segm_start1, segm_end1, segm_start2_2, segm_end2_2,
- new_site1, false, false);
- CHECK_LESS_PREDICATE_SS(segm_start1, segm_end1, segm_start2_2, segm_end2_2,
- new_site1, true, false);
- CHECK_LESS_PREDICATE_SS(segm_start2_2, segm_end2_2, segm_start1, segm_end1,
- new_site1, false, false);
-
- CHECK_LESS_PREDICATE_SS(segm_start1, segm_end1, segm_start2_2, segm_end2_2,
- new_site2, false, true);
- CHECK_LESS_PREDICATE_SS(segm_start1, segm_end1, segm_start2_2, segm_end2_2,
- new_site2, true, true);
- CHECK_LESS_PREDICATE_SS(segm_start2_2, segm_end2_2, segm_start1, segm_end1,
- new_site2, false, true);
-
- point_2d<T> new_site4 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(5));
- CHECK_LESS_PREDICATE_SS(segm_start1, segm_end1, segm_start2_2, segm_end2_2,
- new_site4, false, true);
- CHECK_LESS_PREDICATE_SS(segm_start1, segm_end1, segm_start2_2, segm_end2_2,
- new_site4, true, false);
- CHECK_LESS_PREDICATE_SS(segm_start2_2, segm_end2_2, segm_start1, segm_end1,
- new_site4, false, false);
-}
\ No newline at end of file
+ point_2d<T> new_site2 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(5));
+ point_2d<T> new_site3 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(6));
+ point_2d<T> new_site4 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(8));
+
+ // Common end points.
+ point_2d<T> segm_start2 = make_point_2d<T>(static_cast<T>(-5), static_cast<T>(5));
+ point_2d<T> segm_end2 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(6));
+ site_event_type segm_site2(segm_start2, segm_end2, 1);
+ CHECK_LESS_PREDICATE_SS(segm_site1_2, segm_site2, new_site1, false);
+ CHECK_LESS_PREDICATE_SS(segm_site1_2, segm_site2, new_site2, true);
+ CHECK_LESS_PREDICATE_SS(segm_site1_2, segm_site2, new_site3, true);
+ CHECK_LESS_PREDICATE_SS(segm_site1_2, segm_site2, new_site4, true);
+
+ // No common end points.
+ point_2d<T> segm_start3 = make_point_2d<T>(static_cast<T>(-2), static_cast<T>(4));
+ point_2d<T> segm_end3 = make_point_2d<T>(static_cast<T>(0), static_cast<T>(4));
+ site_event_type segm_site3(segm_start3, segm_end3, 2);
+ CHECK_LESS_PREDICATE_SS(segm_site1_2, segm_site3, new_site1, false);
+ 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();
+ 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);
+ CHECK_LESS_PREDICATE_SS(segm_site3, segm_site1_2, new_site4, true);
+}
Modified: sandbox/SOC/2010/sweepline/libs/sweepline/test/voronoi_segment/segment_voronoi_builder_test.cpp
==============================================================================
--- sandbox/SOC/2010/sweepline/libs/sweepline/test/voronoi_segment/segment_voronoi_builder_test.cpp (original)
+++ sandbox/SOC/2010/sweepline/libs/sweepline/test/voronoi_segment/segment_voronoi_builder_test.cpp 2010-09-28 15:28:09 EDT (Tue, 28 Sep 2010)
@@ -635,7 +635,6 @@
test_voronoi_builder.run_sweepline();
}
-
BOOST_AUTO_TEST_CASE_TEMPLATE(segment_sites_test2, T, test_types) {
typedef T coordinate_type;
voronoi_builder<coordinate_type> test_voronoi_builder;
@@ -708,4 +707,51 @@
test_voronoi_builder.run_sweepline();
// TODO(asydorchuk): Add output checks there.
-}
\ No newline at end of file
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(segment_site_test6, T, test_types) {
+ typedef T coordinate_type;
+ voronoi_builder<coordinate_type> test_voronoi_builder;
+ std::vector< point_2d<T> > point_vec;
+ std::vector< typename voronoi_builder<coordinate_type>::Segment2D > segm_vec;
+ point_2d<T> point1 = make_point_2d<T>(-1, 1);
+ point_2d<T> point2 = make_point_2d<T>(1, 0);
+ point_2d<T> point3 = make_point_2d<T>(1, 2);
+ segm_vec.push_back(std::make_pair< point_2d<T>, point_2d<T> >(point2, point3));
+ point_vec.push_back(point1);
+ test_voronoi_builder.init(point_vec, segm_vec);
+ test_voronoi_builder.run_sweepline();
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(segment_site_test7, T, test_types) {
+ typedef T coordinate_type;
+ voronoi_builder<coordinate_type> test_voronoi_builder;
+ std::vector< point_2d<T> > point_vec;
+ std::vector< typename voronoi_builder<coordinate_type>::Segment2D > segm_vec;
+ point_2d<T> point1 = make_point_2d<T>(0, 0);
+ point_2d<T> point2 = make_point_2d<T>(4, 0);
+ point_2d<T> point3 = make_point_2d<T>(0, 4);
+ point_2d<T> point4 = make_point_2d<T>(4, 4);
+ segm_vec.push_back(std::make_pair(point1, point2));
+ segm_vec.push_back(std::make_pair(point2, point3));
+ segm_vec.push_back(std::make_pair(point3, point4));
+ test_voronoi_builder.init(segm_vec);
+ test_voronoi_builder.run_sweepline();
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(segment_site_test8, T, test_types) {
+ typedef T coordinate_type;
+ voronoi_builder<coordinate_type> test_voronoi_builder;
+ std::vector< point_2d<T> > point_vec;
+ std::vector< typename voronoi_builder<coordinate_type>::Segment2D > segm_vec;
+ point_2d<T> point1 = make_point_2d<T>(0, 0);
+ point_2d<T> point2 = make_point_2d<T>(4, 0);
+ point_2d<T> point3 = make_point_2d<T>(4, 4);
+ point_2d<T> point4 = make_point_2d<T>(0, 4);
+ segm_vec.push_back(std::make_pair(point1, point2));
+ segm_vec.push_back(std::make_pair(point2, point3));
+ segm_vec.push_back(std::make_pair(point3, point4));
+ segm_vec.push_back(std::make_pair(point4, point1));
+ test_voronoi_builder.init(segm_vec);
+ test_voronoi_builder.run_sweepline();
+}
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