Subject: [Boost-commit] svn:boost r70099 - in sandbox/SOC/2010/sweepline: boost/sweepline boost/sweepline/detail libs/sweepline/example libs/sweepline/test
From: sydorchuk.andriy_at_[hidden]
Date: 2011-03-17 21:35:21
Author: asydorchuk
Date: 2011-03-17 21:35:19 EDT (Thu, 17 Mar 2011)
New Revision: 70099
URL: http://svn.boost.org/trac/boost/changeset/70099
Log:
Added mpz based sqr evaluator class.
Added (segment, segment, segment) robust predicate implementation.
Added mpz sqr evaluator tests.
Updated build rules.
Moved to polygon_ulong_long_type instead of unsigned long long.
Added:
sandbox/SOC/2010/sweepline/boost/sweepline/detail/mpz_arithmetic.hpp (contents, props changed)
sandbox/SOC/2010/sweepline/libs/sweepline/test/mpz_arithmetic_test.cpp (contents, props changed)
Text files modified:
sandbox/SOC/2010/sweepline/boost/sweepline/detail/voronoi_formation.hpp | 263 ++++++++++++++++++++++++++++++---------
sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_output.hpp | 26 +-
sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_sweepline.hpp | 13 +
sandbox/SOC/2010/sweepline/libs/sweepline/example/voronoi_visualizer.cpp | 11 +
sandbox/SOC/2010/sweepline/libs/sweepline/test/Jamfile.v2 | 37 +++++
sandbox/SOC/2010/sweepline/libs/sweepline/test/sweepline_test.cpp | 166 ++++++++++++------------
6 files changed, 347 insertions(+), 169 deletions(-)
Added: sandbox/SOC/2010/sweepline/boost/sweepline/detail/mpz_arithmetic.hpp
==============================================================================
--- (empty file)
+++ sandbox/SOC/2010/sweepline/boost/sweepline/detail/mpz_arithmetic.hpp 2011-03-17 21:35:19 EDT (Thu, 17 Mar 2011)
@@ -0,0 +1,277 @@
+// Boost sweepline/mpz_arithmetic.hpp header file
+
+// Copyright Andrii Sydorchuk 2010.
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+// See http://www.boost.org for updates, documentation, and revision history.
+
+#ifndef BOOST_SWEEPLINE_MPZ_ARITHMETIC
+#define BOOST_SWEEPLINE_MPZ_ARITHMETIC
+
+#include <cmath>
+
+#pragma warning(disable:4800)
+#include <gmpxx.h>
+
+namespace boost {
+namespace sweepline {
+namespace detail {
+
+ // Allows to compute expression without allocation of additional memory
+ // for temporary mpz variables. Expression should contain less then SZ
+ // temporary values.
+ class mpz_wrapper {
+ public:
+ mpz_wrapper() {}
+
+ explicit mpz_wrapper(int input) : m_(input) {}
+
+ explicit mpz_wrapper(double input) : m_(input) {}
+
+ mpz_wrapper(const mpz_class& input) : m_(input) {}
+
+ mpz_wrapper(const mpz_wrapper& w) : m_(w.m_) {
+ cur_ = 0;
+ }
+
+ mpz_wrapper& operator=(int that) {
+ m_ = that;
+ return *this;
+ }
+
+ mpz_wrapper& operator=(double that) {
+ m_ = that;
+ return *this;
+ }
+
+ mpz_wrapper& operator=(const mpz_wrapper& that) {
+ m_ = that.m_;
+ return *this;
+ }
+
+ double get_d() const {
+ return m_.get_d();
+ }
+
+ std::string get_str() const {
+ return m_.get_str();
+ }
+
+ bool operator==(const mpz_wrapper& that) const {
+ return m_ == that.m_;
+ }
+
+ bool operator!=(const mpz_wrapper& that) const {
+ return m_ != that.m_;
+ }
+
+ bool operator<(const mpz_wrapper& that) const {
+ return m_ < that.m_;
+ }
+
+ bool operator<=(const mpz_wrapper& that) const {
+ return m_ <= that.m_;
+ }
+
+ bool operator>(const mpz_wrapper& that) const {
+ return m_ > that.m_;
+ }
+
+ bool operator>=(const mpz_wrapper& that) const {
+ return m_ >= that.m_;
+ }
+
+ mpz_wrapper& operator-() const {
+ temp_[cur_].m_ = -this->m_;
+ return temp_[next_cur()];
+ }
+
+ mpz_wrapper& operator+(const mpz_wrapper& that) const {
+ temp_[cur_].m_ = this->m_ + that.m_;
+ return temp_[next_cur()];
+ }
+
+ mpz_wrapper& operator-(const mpz_wrapper& that) const {
+ temp_[cur_].m_ = this->m_ - that.m_;
+ return temp_[next_cur()];
+ }
+
+ mpz_wrapper& operator*(const mpz_wrapper& that) const {
+ temp_[cur_].m_ = this->m_ * that.m_;
+ return temp_[next_cur()];
+ }
+
+ mpz_wrapper& operator+=(const mpz_wrapper& that) {
+ this->m_ += that.m_;
+ return *this;
+ }
+
+ mpz_wrapper& operator-=(const mpz_wrapper& that) {
+ this->m_ -= that.m_;
+ return *this;
+ }
+
+ mpz_wrapper& operator*=(const mpz_wrapper& that) {
+ this->m_ *= that.m_;
+ return *this;
+ }
+
+ bool is_positive() const {
+ return m_.__get_mp()->_mp_size > 0;
+ }
+
+ bool is_negative() const {
+ return m_.__get_mp()->_mp_size < 0;
+ }
+
+ private:
+ static int next_cur() {
+ int ret_val = cur_;
+ cur_ = ++cur_ % SZ;
+ return ret_val;
+ }
+
+ mpz_class m_;
+ static const int SZ = 8;
+ static int cur_;
+ static mpz_wrapper temp_[SZ];
+ };
+
+ int mpz_wrapper::cur_ = 0;
+ mpz_wrapper mpz_wrapper::temp_[mpz_wrapper::SZ];
+
+ static mpz_wrapper& mpz_cross(double a1, double b1, double a2, double b2) {
+ static mpz_wrapper a[2], b[2];
+ a[0] = a1;
+ a[1] = a2;
+ b[0] = b1;
+ b[1] = b2;
+ return a[0] * b[1] - a[1] * b[0];
+ }
+
+ static void swap(mpz_wrapper &a, mpz_wrapper &b) {
+ static mpz_wrapper c;
+ c = a;
+ a = b;
+ b = c;
+ }
+
+ static void sign_sort(mpz_wrapper *A, mpz_wrapper *B, int SZ) {
+ for (int i = SZ-1; i >= 0; --i)
+ for (int j = 0; j < i; ++j)
+ if (!A[j].is_positive() && !A[j+1].is_negative()) {
+ swap(A[j], A[j+1]);
+ swap(B[j], B[j+1]);
+ }
+ }
+
+ template <int N>
+ struct sqr_expr_evaluator {
+ static double eval(mpz_class *A, mpz_class *B);
+ };
+
+ // Evaluates expression:
+ // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) with
+ // 7 * EPS relative error in the worst case.
+ template <>
+ struct sqr_expr_evaluator<2> {
+ static double eval(mpz_wrapper *A, mpz_wrapper *B) {
+ sign_sort(A, B, 2);
+ double ret_val = fabs(A[0].get_d()) * sqrt(B[0].get_d()) +
+ fabs(A[1].get_d()) * sqrt(B[1].get_d());
+ if (!A[1].is_negative())
+ return ret_val;
+ if (!A[0].is_positive())
+ return -ret_val;
+ return (A[0] * A[0] * B[0] - A[1] * A[1] * B[1]).get_d() / ret_val;
+ }
+ };
+
+ // Evaluates expression:
+ // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) + A[2] * sqrt(B[2])
+ // with 14 * EPS relative error in the worst case.
+ template <>
+ struct sqr_expr_evaluator<3> {
+ static double eval(mpz_wrapper *A, mpz_wrapper *B) {
+ sign_sort(A, B, 3);
+ static mpz_wrapper cA[2], cB[2];
+ double res = fabs(A[0].get_d()) * sqrt(B[0].get_d()) +
+ fabs(A[1].get_d()) * sqrt(B[1].get_d()) +
+ fabs(A[2].get_d()) * sqrt(B[2].get_d());
+ if (!A[2].is_negative())
+ return res;
+ if (!A[0].is_positive())
+ return -res;
+ cA[0] = A[0] * A[0] * B[0];
+ cB[0] = 1;
+ if (!A[1].is_negative()) {
+ cA[0] += A[1] * A[1] * B[1] - A[2] * A[2] * B[2];
+ cA[1] = A[0] * A[1];
+ cA[1] = cA[1] + cA[1];
+ cB[1] = B[0] * B[1];
+ } else {
+ cA[0] -= A[1] * A[1] * B[1] + A[2] * A[2] * B[2];
+ cA[1] = A[1] * A[2];
+ cA[1] = -cA[1];
+ cA[1] = cA[1] + cA[1];
+ cB[1] = B[1] * B[2];
+ }
+ return sqr_expr_evaluator<2>::eval(cA, cB) / res;
+ }
+ };
+
+ // Evaluates expression:
+ // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) + A[2] * sqrt(B[2]) + A[3] * sqrt(B[3])
+ // with 23 * EPS relative error in the worst case.
+ template <>
+ struct sqr_expr_evaluator<4> {
+ static double eval(mpz_wrapper *A, mpz_wrapper *B) {
+ sign_sort(A, B, 4);
+ static mpz_wrapper cA[3], cB[3];
+ double res = fabs(A[0].get_d()) * sqrt(B[0].get_d()) +
+ fabs(A[1].get_d()) * sqrt(B[1].get_d()) +
+ fabs(A[2].get_d()) * sqrt(B[2].get_d()) +
+ fabs(A[3].get_d()) * sqrt(B[3].get_d());
+ if (!A[3].is_negative())
+ return res;
+ if (!A[0].is_positive())
+ return -res;
+ bool three_neg = !A[1].is_positive();
+ bool three_pos = !A[2].is_negative();
+ if (!A[1].is_positive() || !A[2].is_negative()) {
+ if ((A[1] * A[1] * B[1] <= A[2] * A[2] * B[2]) ^ A[2].is_negative()) {
+ swap(A[1], A[2]);
+ swap(B[1], B[2]);
+ }
+ }
+ cA[0] = 0;
+ for (int i = 0; i < 4; ++i)
+ if (i < 2)
+ cA[0] += A[i] * A[i] * B[i];
+ else
+ cA[0] -= A[i] * A[i] * B[i];
+ cA[1] = A[0] * A[1];
+ cA[1] += cA[1];
+ cA[2] = A[2] * A[3];
+ cA[2] += cA[2];
+ cA[2] = -cA[2];
+ cB[0] = 1;
+ cB[1] = B[0] * B[1];
+ cB[2] = B[2] * B[3];
+ double numerator = sqr_expr_evaluator<3>::eval(cA, cB);
+ if (!three_pos && !three_neg)
+ return numerator / res;
+ double denominator = fabs(A[0].get_d()) * sqrt(B[0].get_d()) +
+ fabs(A[3].get_d()) * sqrt(B[3].get_d());
+ A[2] = -A[2];
+ return numerator / (denominator + sqr_expr_evaluator<2>::eval(&A[1], &B[1]));
+ }
+ };
+
+} // detail
+} // sweepline
+} // boost
+
+#endif
Modified: sandbox/SOC/2010/sweepline/boost/sweepline/detail/voronoi_formation.hpp
==============================================================================
--- sandbox/SOC/2010/sweepline/boost/sweepline/detail/voronoi_formation.hpp (original)
+++ sandbox/SOC/2010/sweepline/boost/sweepline/detail/voronoi_formation.hpp 2011-03-17 21:35:19 EDT (Thu, 17 Mar 2011)
@@ -389,8 +389,8 @@
erc_type rhs1 = denom_ * that.lower_x_;
erc_type lhs2 = center_y_ * that.denom_;
erc_type rhs2 = denom_ * that.center_y_;
- return (lhs1.compare(rhs1, 64) == UNDEFINED &&
- lhs2.compare(rhs2, 64) == UNDEFINED);
+ return (lhs1.compare(rhs1, 128) == UNDEFINED &&
+ lhs2.compare(rhs2, 128) == UNDEFINED);
}
bool operator!=(const circle_event &that) const {
@@ -407,7 +407,7 @@
// Compare x-coordinates of the rightmost points. If the result
// is undefined we assume that x-coordinates are equal.
- kPredicateResult pres = lhs1.compare(rhs1, 64);
+ kPredicateResult pres = lhs1.compare(rhs1, 128);
if (pres != UNDEFINED)
return (pres == LESS);
@@ -416,7 +416,7 @@
erc_type rhs2 = denom_ * that.center_y_;
// Compare y-coordinates of the rightmost points.
- return (lhs2.compare(rhs2, 64) == LESS);
+ return (lhs2.compare(rhs2, 128) == LESS);
}
bool operator<=(const circle_event &that) const {
@@ -632,12 +632,12 @@
// Return true if the value is positive, else false.
template <typename T>
static inline bool convert_to_65_bit(
- T value, BOOST_POLYGON_UNSIGNED_LONG_LONG &res) {
+ T value, polygon_ulong_long_type &res) {
if (value >= 0) {
- res = static_cast<BOOST_POLYGON_UNSIGNED_LONG_LONG>(value);
+ res = static_cast<polygon_ulong_long_type>(value);
return true;
} else {
- res = static_cast<BOOST_POLYGON_UNSIGNED_LONG_LONG>(-value);
+ res = static_cast<polygon_ulong_long_type>(-value);
return false;
}
}
@@ -648,7 +648,7 @@
// their integer reinterpretatoins differ in not more than maxUlps units.
static inline bool almost_equal(double a, double b,
unsigned int maxUlps) {
- BOOST_POLYGON_UNSIGNED_LONG_LONG ll_a, ll_b;
+ polygon_ulong_long_type ll_a, ll_b;
// Reinterpret double bits as 64-bit signed integer.
memcpy(&ll_a, &a, sizeof(double));
@@ -680,15 +680,15 @@
template <typename T>
static kOrientation orientation_test(T dif_x1_, T dif_y1_,
T dif_x2_, T dif_y2_) {
- BOOST_POLYGON_UNSIGNED_LONG_LONG dif_x1, dif_y1, dif_x2, dif_y2;
+ polygon_ulong_long_type dif_x1, dif_y1, dif_x2, dif_y2;
bool dif_x1_plus, dif_x2_plus, dif_y1_plus, dif_y2_plus;
dif_x1_plus = convert_to_65_bit(dif_x1_, dif_x1);
dif_y1_plus = convert_to_65_bit(dif_y1_, dif_y1);
dif_x2_plus = convert_to_65_bit(dif_x2_, dif_x2);
dif_y2_plus = convert_to_65_bit(dif_y2_, dif_y2);
- BOOST_POLYGON_UNSIGNED_LONG_LONG expr_l = dif_x1 * dif_y2;
- BOOST_POLYGON_UNSIGNED_LONG_LONG expr_r = dif_x2 * dif_y1;
+ polygon_ulong_long_type expr_l = dif_x1 * dif_y2;
+ polygon_ulong_long_type expr_r = dif_x2 * dif_y1;
bool expr_l_plus = (dif_x1_plus == dif_y2_plus) ? true : false;
bool expr_r_plus = (dif_x2_plus == dif_y1_plus) ? true : false;
@@ -743,19 +743,19 @@
}
// Compute robust cross_product: a1 * b2 - b1 * a2.
- // The result is correct with epsilon relative error equal to 1EPS.
+ // The result is correct with epsilon relative error equal to 2EPS.
template <typename T>
static double robust_cross_product(T a1_, T b1_, T a2_, T b2_) {
- BOOST_POLYGON_UNSIGNED_LONG_LONG a1, b1, a2, b2;
+ polygon_ulong_long_type a1, b1, a2, b2;
bool a1_plus, a2_plus, b1_plus, b2_plus;
a1_plus = convert_to_65_bit(a1_, a1);
b1_plus = convert_to_65_bit(b1_, b1);
a2_plus = convert_to_65_bit(a2_, a2);
b2_plus = convert_to_65_bit(b2_, b2);
- BOOST_POLYGON_UNSIGNED_LONG_LONG expr_l = a1 * b2;
+ polygon_ulong_long_type expr_l = a1 * b2;
bool expr_l_plus = (a1_plus == b2_plus) ? true : false;
- BOOST_POLYGON_UNSIGNED_LONG_LONG expr_r = b1 * a2;
+ polygon_ulong_long_type expr_r = b1 * a2;
bool expr_r_plus = (a2_plus == b1_plus) ? true : false;
if (expr_l == 0)
@@ -1021,7 +1021,7 @@
// Find the x-coordinate (relative to the sweepline position) of the point
// of the intersection of the horizontal line going through the new site
// with the arc corresponding to the segment site.
- // The relative error is atmost 7EPS.
+ // The relative error is atmost 8EPS.
template <typename T>
static double find_distance_to_segment_arc(const site_event<T> &segment,
const point_2d<T> &new_point) {
@@ -1049,9 +1049,9 @@
k = 1.0 / (b1 - k);
}
}
- // Relative error of the robust cross product is 1EPS.
+ // Relative error of the robust cross product is 2EPS.
// Relative error of the k is atmost 5EPS.
- // The resulting relative error is atmost 7EPS.
+ // The resulting relative error is atmost 8EPS.
return robust_cross_product(a1, b1, a3, b3) * k;
}
}
@@ -1063,8 +1063,8 @@
static bool robust_65bit_less_predicate(const point_2d<T> &left_point,
const point_2d<T> &right_point,
const point_2d<T> &new_point) {
- BOOST_POLYGON_UNSIGNED_LONG_LONG a1, a2, b1, b2;
- BOOST_POLYGON_UNSIGNED_LONG_LONG b1_sqr, b2_sqr, l_expr, r_expr;
+ polygon_ulong_long_type a1, a2, b1, b2;
+ polygon_ulong_long_type b1_sqr, b2_sqr, l_expr, r_expr;
bool l_expr_plus, r_expr_plus;
// Compute a1 = (x0-x1), a2 = (x0-x2), b1 = (y0 - y1), b2 = (y0 - y2).
@@ -1076,8 +1076,8 @@
// Compute b1_sqr = (y0-y1)^2 and b2_sqr = (y0-y2)^2.
b1_sqr = b1 * b1;
b2_sqr = b2 * b2;
- BOOST_POLYGON_UNSIGNED_LONG_LONG b1_sqr_mod = b1_sqr % a1;
- BOOST_POLYGON_UNSIGNED_LONG_LONG b2_sqr_mod = b2_sqr % a2;
+ polygon_ulong_long_type b1_sqr_mod = b1_sqr % a1;
+ polygon_ulong_long_type b2_sqr_mod = b2_sqr % a2;
// Compute l_expr = (x1 - x2).
l_expr_plus = convert_to_65_bit(left_point.x() - right_point.x(), l_expr);
@@ -1085,18 +1085,18 @@
// In case (b2_sqr / a2) < (b1_sqr / a1), decrease left_expr by 1.
if (b2_sqr_mod * a1 < b1_sqr_mod * a2) {
if (!l_expr_plus)
- l_expr++;
+ ++l_expr;
else if (l_expr != 0)
- l_expr--;
+ --l_expr;
else {
- l_expr++;
+ ++l_expr;
l_expr_plus = false;
}
}
// Compute right expression.
- BOOST_POLYGON_UNSIGNED_LONG_LONG value1 = b1_sqr / a1;
- BOOST_POLYGON_UNSIGNED_LONG_LONG value2 = b2_sqr / a2;
+ polygon_ulong_long_type value1 = b1_sqr / a1;
+ polygon_ulong_long_type value2 = b2_sqr / a2;
if (value1 >= value2) {
r_expr = value1 - value2;
r_expr_plus = true;
@@ -1243,10 +1243,11 @@
double dist1 = find_distance_to_point_arc(site_point, new_point);
double dist2 = find_distance_to_segment_arc(segment_site, new_point);
- if (!almost_equal(dist1, dist2, 10)) {
+ if (!almost_equal(dist1, dist2, 11)) {
return reverse_order ^ (dist1 < dist2);
}
+ // TODO(asydorchuk): Add mpl support there.
return reverse_order ^ (dist1 < dist2);
}
@@ -1271,8 +1272,8 @@
double dist1 = find_distance_to_segment_arc(left_segment, new_point);
double dist2 = find_distance_to_segment_arc(right_segment, new_point);
- // The undefined ulp range is equal to 7EPS + 7EPS <= 14ULP.
- if (!almost_equal(dist1, dist2, 14)) {
+ // The undefined ulp range is equal to 8EPS + 8EPS <= 16ULP.
+ if (!almost_equal(dist1, dist2, 16)) {
return dist1 < dist2;
}
@@ -1294,7 +1295,79 @@
static inline T sqr_distance(T dif_x, T dif_y) {
return dif_x * dif_x + dif_y * dif_y;
}
-
+/*
+ // Recompute parameters of the circle event using high-precision library.
+ // In the worst case parameters of the circle have next relative errors:
+ // r(c_x) = 4 * EPS;
+ // r(c_y) = 4 * EPS;
+ // r(lower_x) = 8 * EPS.
+ template <typename T>
+ static bool create_circle_event_ppp_gmpxx(const site_event<T> &site1,
+ const site_event<T> &site2,
+ const site_event<T> &site3,
+ circle_event<T> &c_event) {
+ static mpz_wrapper mpz_dif_x[3], mpz_dif_y[3], mpz_sum_x[3], mpz_sum_y[3],
+ mpz_numerator[3], mpz_c_x, mpz_c_y, mpz_sqr_r;
+ mpz_dif_x[0] = site1.x() - site2.x();
+ mpz_dif_x[1] = site2.x() - site3.x();
+ mpz_dif_x[2] = site1.x() - site3.x();
+ mpz_dif_y[0] = site1.y() - site2.y();
+ mpz_dif_y[1] = site2.y() - site3.y();
+ mpz_dif_y[2] = site1.y() - site3.y();
+
+ // Evaluate orientation.
+ double orientation = static_cast<double>(
+ mpz_dif_x[0] * mpz_dif_y[1] - mpz_dif_x[1] * mpz_dif_y[0]);
+
+ // If use this function only to recompute parameters of the circle
+ // event, this check won't be required.
+ if (orientation_test(orientation) != RIGHT_ORIENTATION)
+ return false;
+
+ // Evaluate inverse orientation to avoid division in calculations.
+ // r(inv_orientation) = 2 * EPS.
+ double inv_orientation = 0.5 / orientation;
+
+ mpz_sum_x[0] = site1.x() + site2.x();
+ mpz_sum_x[1] = site2.x() + site3.x();
+ mpz_sum_y[0] = site1.y() + site2.y();
+ mpz_sum_y[1] = site2.y() + site3.y();
+
+ mpz_numerator[1] = mpz_dif_x[0] * mpz_sum_x[0] + mpz_dif_y[0] * mpz_sum_y[0];
+ mpz_numerator[2] = mpz_dif_x[1] * mpz_sum_x[1] + mpz_dif_y[1] * mpz_sum_y[1];
+
+ mpz_c_x = mpz_numerator[1] * mpz_dif_y[1] - mpz_numerator[2] * mpz_dif_y[0];
+ mpz_c_y = mpz_numerator[2] * mpz_dif_x[0] - mpz_numerator[1] * mpz_dif_x[1];
+
+ // Evaluate radius of the circle.
+ mpz_sqr_r = 1.0;
+ for (int i = 0; i < 3; ++i)
+ mpz_sqr_r *= mpz_dif_x[i] * mpz_dif_x[i] + mpz_dif_y[i] * mpz_dif_y[i];
+
+ // Evaluate coordinates of the center of the circle.
+ // r(c_x) = r(c_y) = 4 * EPS.
+ double c_x = static_cast<double>(mpz_c_x) * inv_orientation;
+ double c_y = static_cast<double>(mpz_c_y) * inv_orientation;
+
+ // r(r) = 1.5 * EPS < 2 * EPS.
+ double r = sqrt(static_cast<double>(mpz_sqr_r));
+
+ // If c_x >= 0 then lower_x = c_x + r,
+ // else lower_x = (c_x * c_x - r * r) / (c_x - r).
+ // To guarantee epsilon relative error.
+ if (c_x >= 0) {
+ // r(lower_x) = 5 * EPS.
+ c_event = circle_event<double>(c_x, c_y, c_x + r * fabs(inv_orientation));
+ } else {
+ mpz_numerator[0] = mpz_c_x * mpz_c_x - mpz_sqr_r;
+ // r(lower_x) = 8 * EPS.
+ double lower_x = static_cast<double>(mpz_numerator[0]) * fabs(inv_orientation);
+ lower_x /= (mpz_c_x >= 0) ? (-mpz_c_x.get_d() - r) : (mpz_c_x.get_d() - r);
+ c_event = circle_event<double>(c_x, c_y, lower_x);
+ }
+ return true;
+ }
+*/
// Find parameters of the inscribed circle that is tangent to three
// point sites.
template <typename T>
@@ -1482,6 +1555,7 @@
ix += a1 * c2 * inv_orientation;
iy += b1 * c2 * inv_orientation;
iy += b2 * c1 * inv_orientation;
+
b += ix * (a1 * sqr_sum2);
b += ix * (a2 * sqr_sum1);
b += iy * (b1 * sqr_sum2);
@@ -1507,7 +1581,60 @@
}
return true;
}
-
+/*
+ template <typename T>
+ static bool create_circle_event_sss_gmpxx(const site_event<T> &site1,
+ const site_event<T> &site2,
+ const site_event<T> &site3,
+ circle_event<T> &c_event) {
+ if (site1.index() == site2.index() || site2.index() == site3.index())
+ return false;
+ static mpz_wrapper a[3], b[3], c[3], sqr_len[4], cross[4];
+ a[0] = site1.x1(true) - site1.x0(true);
+ a[1] = site2.x1(true) - site2.x0(true);
+ a[2] = site3.x1(true) - site3.x0(true);
+
+ b[0] = site1.y1(true) - site1.y0(true);
+ b[1] = site2.y1(true) - site2.y0(true);
+ b[2] = site3.y1(true) - site3.y0(true);
+
+ c[0] = mpz_cross(site1.x0(true), site1.y0(true), site1.x1(true), site1.y1(true));
+ c[1] = mpz_cross(site2.x0(true), site2.y0(true), site2.x1(true), site2.y1(true));
+ c[2] = mpz_cross(site3.x0(true), site3.y0(true), site3.x1(true), site3.y1(true));
+
+ for (int i = 0; i < 3; ++i) {
+ int j = (i+1) % 3;
+ int k = (i+2) % 3;
+ cross[i] = a[j] * b[k] - a[k] * b[j];
+ sqr_len[i] = a[i] * a[i] + b[i] * b[i];
+ }
+ double denom = sqr_expr_evaluator<3>::eval(cross, sqr_len);
+
+ for (int i = 0; i < 3; ++i) {
+ int j = (i+1) % 3;
+ int k = (i+2) % 3;
+ cross[i] = b[j] * c[k] - b[k] * c[j];
+ sqr_len[i] = a[i] * a[i] + b[i] * b[i];
+ }
+ double c_y = sqr_expr_evaluator<3>::eval(cross, sqr_len) / denom;
+
+ cross[3] = 0;
+ for (int i = 0; i < 3; ++i) {
+ int j = (i+1) % 3;
+ int k = (i+2) % 3;
+ cross[i] = a[j] * c[k] - a[k] * c[j];
+ sqr_len[i] = a[i] * a[i] + b[i] * b[i];
+ cross[3] += cross[i] * b[i];
+ }
+ double c_x = sqr_expr_evaluator<3>::eval(cross, sqr_len) / denom;
+
+ sqr_len[3] = 1;
+ double lower_x = sqr_expr_evaluator<4>::eval(cross, sqr_len) / denom;
+
+ c_event = circle_event<double>(c_x, c_y, lower_x);
+ return true;
+ }
+*/
// Find parameters of the inscribed circle that is tangent to three
// segment sites.
template <typename T>
@@ -1515,28 +1642,39 @@
const site_event<T> &site2,
const site_event<T> &site3,
circle_event<T> &c_event) {
- if (site1.index() == site2.index() || site2.index() == site3.index()) {
+ // TODO(asydorchuk): Check if this one is required.
+ if (site1.index() == site2.index() || site2.index() == site3.index())
return false;
- }
double a1 = site1.x1(true) - site1.x0(true);
double b1 = site1.y1(true) - site1.y0(true);
- double c1 = robust_cross_product(b1, a1, site1.y0(true), site1.x0(true));
+ double c1 = robust_cross_product(site1.x0(true), site1.y0(true),
+ site1.x1(true), site1.y1(true));
double a2 = site2.x1(true) - site2.x0(true);
double b2 = site2.y1(true) - site2.y0(true);
- double c2 = robust_cross_product(b2, a2, site2.y0(true), site2.x0(true));
+ double c2 = robust_cross_product(site2.x0(true), site2.y0(true),
+ site2.x1(true), site2.y1(true));
double a3 = site3.x1(true) - site3.x0(true);
double b3 = site3.y1(true) - site3.y0(true);
- double c3 = robust_cross_product(b3, a3, site3.y0(true), site3.x0(true));
+ double c3 = robust_cross_product(site3.x0(true), site3.y0(true),
+ site3.x1(true), site3.y1(true));
double len1 = sqrt(a1 * a1 + b1 * b1);
double len2 = sqrt(a2 * a2 + b2 * b2);
double len3 = sqrt(a3 * a3 + b3 * b3);
double cross_12 = robust_cross_product(a1, b1, a2, b2);
double cross_23 = robust_cross_product(a2, b2, a3, b3);
double cross_31 = robust_cross_product(a3, b3, a1, b1);
- epsilon_robust_comparator<double> det, c_x, c_y, r;
- det += cross_12 * len3;
- det += cross_23 * len1;
- det += cross_31 * len2;
+ epsilon_robust_comparator<double> denom, c_x, c_y, r;
+
+ // denom = cross_12 * len3 + cross_23 * len1 + cross_31 * len2.
+ denom += cross_12 * len3;
+ denom += cross_23 * len1;
+ denom += cross_31 * len2;
+
+ // denom * r = (b2 * c_x - a2 * c_y - c2 * denom) / len2.
+ r -= cross_12 * c3;
+ r -= cross_23 * c1;
+ r -= cross_31 * c2;
+
c_x += a1 * c2 * len3;
c_x -= a2 * c1 * len3;
c_x += a2 * c3 * len1;
@@ -1549,11 +1687,7 @@
c_y -= b3 * c2 * len1;
c_y += b3 * c1 * len2;
c_y -= b1 * c3 * len2;
- r += b2 * c_x;
- r -= a2 * c_y;
- r -= c2 * det;
- r /= len2;
- c_event = circle_event<double>(c_x, c_y, c_x + r, det);
+ c_event = circle_event<double>(c_x, c_y, c_x + r, denom);
return true;
}
@@ -1711,6 +1845,7 @@
void edge(voronoi_edge<T> *new_edge) {
edge_ = new_edge;
}
+
private:
typename circle_events_queue<T>::circle_event_type_it circle_event_it_;
voronoi_edge<T> *edge_;
@@ -1850,7 +1985,7 @@
// Create a site event from each input point.
for (typename std::vector< point_2d<iType> >::const_iterator it = points.begin();
- it != points.end(); it++) {
+ it != points.end(); ++it) {
site_events_.push_back(make_site_event(static_cast<T>(it->x()),
static_cast<T>(it->y()),
0));
@@ -1861,7 +1996,7 @@
// 2) the endpoint of the segment;
// 3) the segment itself.
for (typename std::vector<iSegment2D>::const_iterator it = segments.begin();
- it != segments.end(); it++) {
+ it != segments.end(); ++it) {
T x1 = static_cast<T>(it->first.x());
T y1 = static_cast<T>(it->first.y());
T x2 = static_cast<T>(it->second.x());
@@ -1879,7 +2014,7 @@
site_events_.begin(), site_events_.end()), site_events_.end());
// Number the sites.
- for (size_t cur = 0; cur < site_events_.size(); cur++)
+ for (size_t cur = 0; cur < site_events_.size(); ++cur)
site_events_[cur].index(cur);
// Init the site's iterator.
@@ -1904,7 +2039,7 @@
} else if (site_events_.size() == 1) {
// Handle one input site case.
output_.process_single_site(site_events_[0]);
- site_event_iterator_++;
+ ++site_event_iterator_;
} else {
int skip = 0;
@@ -1912,8 +2047,8 @@
while(site_event_iterator_ != site_events_.end() &&
site_event_iterator_->x() == site_events_.begin()->x() &&
site_event_iterator_->is_vertical()) {
- site_event_iterator_++;
- skip++;
+ ++site_event_iterator_;
+ ++skip;
}
if (skip == 1) {
@@ -1969,21 +2104,21 @@
// Get the first and the second site events.
site_event_iterator_type it_first = site_events_.begin();
site_event_iterator_type it_second = site_events_.begin();
- it_second++;
+ ++it_second;
// Update the beach line.
insert_new_arc(*it_first, *it_first, *it_second, beach_line_.begin());
// The second site has been already processed.
// Move the site's iterator.
- site_event_iterator_++;
+ ++site_event_iterator_;
}
// Insert bisectors for all collinear initial sites.
void init_beach_line_collinear_sites() {
site_event_iterator_type it_first = site_events_.begin();
site_event_iterator_type it_second = site_events_.begin();
- it_second++;
+ ++it_second;
while (it_second != site_event_iterator_) {
// Create a new beach line node.
key_type new_node(*it_first, *it_second);
@@ -1996,8 +2131,8 @@
std::pair<key_type, value_type>(new_node, value_type(edge)));
// Update iterators.
- it_first++;
- it_second++;
+ ++it_first;
+ ++it_second;
}
}
@@ -2007,7 +2142,7 @@
site_event_type site_event = *site_event_iterator_;
// Move the site's iterator.
- site_event_iterator_++;
+ ++site_event_iterator_;
// If a new site is an end point of some segment,
// remove temporary nodes from the beach line data structure.
@@ -2033,7 +2168,7 @@
if (it == beach_line_.end()) {
// The above arc corresponds to the second arc of the last node.
// Move the iterator to the last node.
- it--;
+ --it;
// Get the second site of the last node
const site_event_type &site_arc = it->first.right_site();
@@ -2074,7 +2209,7 @@
// Remove the candidate circle from the event queue.
it->second.deactivate_circle_event();
- it--;
+ --it;
const site_event_type &site_arc1 = it->first.right_site();
const site_event_type &site1 = it->first.left_site();
@@ -2122,7 +2257,7 @@
edge_type *bisector2 = it_first->second.edge();
// Get the half-edge corresponding to the first bisector - (A, B).
- it_first--;
+ --it_first;
edge_type *bisector1 = it_first->second.edge();
// Get the A site.
@@ -2151,14 +2286,14 @@
// to the left for potential circle events.
if (it_first != beach_line_.begin()) {
it_first->second.deactivate_circle_event();
- it_first--;
+ --it_first;
const site_event_type &site_l1 = it_first->first.left_site();
activate_circle_event(site_l1, site1, site3, it_last);
}
// Check the new triplet formed by the neighboring arcs
// to the right for potential circle events.
- it_last++;
+ ++it_last;
if (it_last != beach_line_.end()) {
it_last->second.deactivate_circle_event();
const site_event_type &site_r1 = it_last->first.right_site();
@@ -2293,8 +2428,8 @@
void operator=(const voronoi_builder&);
};
+} // detail
} // sweepline
} // boost
-} // detail
#endif
Modified: sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_output.hpp
==============================================================================
--- sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_output.hpp (original)
+++ sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_output.hpp 2011-03-17 21:35:19 EDT (Thu, 17 Mar 2011)
@@ -404,7 +404,7 @@
if (point_site == segment_site_start ||
point_site == segment_site_end)
return;
-
+
// Apply the linear transformation to move start point of the
// segment to the point with coordinates (0, 0) and the direction
// of the segment to coincide the positive direction of the x-axis.
@@ -606,8 +606,8 @@
friend class voronoi_output<coordinate_type>;
void incident_edge(voronoi_edge_type *e) { incident_edge_ = e; }
- void inc_num_incident_edges() { num_incident_edges_++; }
- void dec_num_incident_edges() { num_incident_edges_--; }
+ void inc_num_incident_edges() { ++num_incident_edges_; }
+ void dec_num_incident_edges() { --num_incident_edges_; }
point_2d_type point0_;
point_2d_type point1_;
@@ -1040,13 +1040,13 @@
voronoi_edge_type *edge1 = &(*edge_it);
edge1->next(edge1);
edge1->prev(edge1);
- edge_it++;
+ ++edge_it;
edge1 = &(*edge_it);
- edge_it++;
+ ++edge_it;
while (edge_it != edge_records_.end()) {
voronoi_edge_type *edge2 = &(*edge_it);
- edge_it++;
+ ++edge_it;
edge1->next(edge2);
edge1->prev(edge2);
@@ -1054,7 +1054,7 @@
edge2->prev(edge1);
edge1 = &(*edge_it);
- edge_it++;
+ ++edge_it;
}
edge1->next(edge1);
@@ -1071,7 +1071,7 @@
// Degenerate edges exist only among finite edges.
if (!edge_it1->vertex0() || !edge_it1->vertex1()) {
- num_edge_records_++;
+ ++num_edge_records_;
continue;
}
@@ -1079,7 +1079,7 @@
const voronoi_vertex_type *v2 = edge_it1->vertex1();
// Make epsilon robust check.
- if (v1->robust_vertex()->equals(v2->robust_vertex(), 64)) {
+ if (v1->robust_vertex()->equals(v2->robust_vertex(), 128)) {
// Decrease number of cell's incident edges.
edge_it1->cell()->dec_num_incident_edges();
edge_it1->twin()->cell()->dec_num_incident_edges();
@@ -1122,7 +1122,7 @@
edge_records_.erase(edge_it1, edge_it);
} else {
// Count not degenerate edge.
- num_edge_records_++;
+ ++num_edge_records_;
}
}
robust_vertices_.clear();
@@ -1134,14 +1134,14 @@
if (vertex_it->incident_edge() == NULL) {
vertex_it = vertex_records_.erase(vertex_it);
} else {
- vertex_it++;
- num_vertex_records_++;
+ ++vertex_it;
+ ++num_vertex_records_;
}
}
// Update prev/next pointers for the ray-edges.
for (voronoi_cell_iterator_type cell_it = cell_records_.begin();
- cell_it != cell_records_.end(); cell_it++) {
+ cell_it != cell_records_.end(); ++cell_it) {
// Move to the previous edge while
// it is possible in the CW direction.
voronoi_edge_type *cur_edge = cell_it->incident_edge();
Modified: sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_sweepline.hpp
==============================================================================
--- sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_sweepline.hpp (original)
+++ sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_sweepline.hpp 2011-03-17 21:35:19 EDT (Thu, 17 Mar 2011)
@@ -19,14 +19,17 @@
#include <stack>
#include <vector>
-#ifdef USE_MULTI_PRECISION_LIBRARY
- #pragma warning( disable : 4800 )
- #include "gmpxx.h"
+#ifndef BOOST_POLYGON_USE_LONG_LONG
+#include <boost/config.hpp>
+#ifdef BOOST_HAS_LONG_LONG
+#define BOOST_POLYGON_USE_LONG_LONG
+typedef boost::long_long_type polygon_long_long_type;
+typedef boost::ulong_long_type polygon_ulong_long_type;
+#endif
#endif
-
-#define BOOST_POLYGON_UNSIGNED_LONG_LONG unsigned long long
#include "voronoi_output.hpp"
+//#include "detail/mpz_arithmetic.hpp"
#include "detail/voronoi_formation.hpp"
namespace boost {
Modified: sandbox/SOC/2010/sweepline/libs/sweepline/example/voronoi_visualizer.cpp
==============================================================================
--- sandbox/SOC/2010/sweepline/libs/sweepline/example/voronoi_visualizer.cpp (original)
+++ sandbox/SOC/2010/sweepline/libs/sweepline/example/voronoi_visualizer.cpp 2011-03-17 21:35:19 EDT (Thu, 17 Mar 2011)
@@ -7,6 +7,8 @@
// See http://www.boost.org for updates, documentation, and revision history.
+#include <iostream>
+
#include <QtOpenGL/QGLWidget>
#include <QtGui/QtGui>
@@ -68,8 +70,15 @@
}
protected:
+ void initializeGL() {
+ glHint(GL_POINT_SMOOTH_HINT, GL_NICEST);
+ glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
+ glEnable(GL_BLEND);
+ glEnable(GL_POINT_SMOOTH);
+ }
+
void paintGL() {
- qglClearColor(QColor::fromRgb(0, 0, 0));
+ qglClearColor(QColor::fromRgb(0, 0, 0));
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
// Draw voronoi sites.
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 2011-03-17 21:35:19 EDT (Thu, 17 Mar 2011)
@@ -11,7 +11,7 @@
;
-test-suite "sweepline"
+alias "sweepline"
:
[ run circle_event_creation_test.cpp ]
[ run clipping_test.cpp ]
@@ -22,7 +22,38 @@
[ run sweepline_test.cpp ]
;
-test-suite "sweepline_benchmark"
+alias "sweepline_benchmark"
:
[ run benchmark_test.cpp ]
- ;
\ No newline at end of file
+ ;
+
+path-constant GMP_PATH : "/usr/local" ;
+
+lib gmp
+ : # sources
+ : # requirements
+ <name>gmp <search>$(GMP_PATH)/lib
+ : #default-build
+ : # usage-requirements
+ ;
+
+lib gmpxx
+ : # sources
+ : # requirements
+ <name>gmpxx <search>$(GMP_PATH)/lib
+ : # default-build
+ : # usage-requirements
+ ;
+
+alias gmp_libs : gmp gmpxx : <link>static ;
+
+alias "mpz_arithmetic"
+ :
+ [ run mpz_arithmetic_test.cpp gmp_libs
+ : # args
+ : # input files
+ : # requirements
+ <include>$(GMP_PATH)/include
+ : # target-name
+ ]
+ ;
Added: sandbox/SOC/2010/sweepline/libs/sweepline/test/mpz_arithmetic_test.cpp
==============================================================================
--- (empty file)
+++ sandbox/SOC/2010/sweepline/libs/sweepline/test/mpz_arithmetic_test.cpp 2011-03-17 21:35:19 EDT (Thu, 17 Mar 2011)
@@ -0,0 +1,67 @@
+// Boost sweepline library mpz_arithmetic_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 <stdlib.h>
+#include <time.h>
+
+#include "test_type_list.hpp"
+#include "boost/sweepline/detail/mpz_arithmetic.hpp"
+using namespace boost::sweepline::detail;
+
+#define BOOST_TEST_MODULE mpz_arithmetic_test
+#include <boost/test/test_case_template.hpp>
+
+BOOST_AUTO_TEST_CASE(mpz_swap_test) {
+ mpz_wrapper a(1.0);
+ mpz_wrapper b(2.0);
+ swap(a, b);
+ BOOST_CHECK(a == mpz_wrapper(2.0));
+ BOOST_CHECK(b == mpz_wrapper(1.0));
+}
+
+template <int N>
+struct test_sqr_evaluator {
+ static bool run() {
+ bool ret_val = true;
+ static mpz_wrapper A[N], B[N];
+ double a[N], b[N];
+ for (int i = 0; i < N; ++i) {
+ a[i] = rand() % 10 + 5;
+ int rand_number = rand() % 10 + 5;
+ b[i] = rand_number * rand_number;
+ }
+ int mask = (1 << N);
+ for (int i = 0; i < mask; i++) {
+ double expected_val = 0.0;
+ for (int j = 0; j < N; j++) {
+ if (i & (1 << j)) {
+ A[j] = a[j];
+ B[j] = b[j];
+ expected_val += a[j] * sqrt(b[j]);
+ } else {
+ A[j] = -a[j];
+ B[j] = b[j];
+ expected_val -= a[j] * sqrt(b[j]);
+ }
+ }
+ double received_val = sqr_expr_evaluator<N>::eval(A, B);
+ ret_val &= (fabs(expected_val - received_val) <= 1E-9);
+ }
+ return ret_val;
+ }
+};
+
+BOOST_AUTO_TEST_CASE(mpz_sqr_evaluator_test) {
+ srand(static_cast<unsigned int>(time(0)));
+ for (int i = 0; i < 100; i++) {
+ BOOST_CHECK(test_sqr_evaluator<2>::run());
+ BOOST_CHECK(test_sqr_evaluator<3>::run());
+ BOOST_CHECK(test_sqr_evaluator<4>::run());
+ }
+}
Modified: sandbox/SOC/2010/sweepline/libs/sweepline/test/sweepline_test.cpp
==============================================================================
--- sandbox/SOC/2010/sweepline/libs/sweepline/test/sweepline_test.cpp (original)
+++ sandbox/SOC/2010/sweepline/libs/sweepline/test/sweepline_test.cpp 2011-03-17 21:35:19 EDT (Thu, 17 Mar 2011)
@@ -600,8 +600,8 @@
BOOST_AUTO_TEST_CASE_TEMPLATE(segment_grid_test, T, test_types) {
voronoi_output<T> test_output_small, test_output_large;
std::vector< std::pair< point_2d<T>, point_2d<T> > > segm_vec_small, segm_vec_large;
- int grid_size[4] = {10, 33, 100, 333};
- int max_value[4] = {100, 330, 1000, 3330};
+ int grid_size[] = {10, 33, 100, 333};
+ int max_value[] = {100, 330, 1000, 3330};
int array_length = sizeof(grid_size) / sizeof(int);
for (int k = 0; k < array_length; k++) {
int cur_sz = grid_size[k];
@@ -637,85 +637,85 @@
}
#endif
-#ifdef NDEBUG
-BOOST_AUTO_TEST_CASE_TEMPLATE(segment_random_test1, T, test_types) {
- srand(static_cast<unsigned int>(time(NULL)));
- voronoi_output<T> test_output;
- std::vector< point_2d<T> > point_vec;
- std::vector< std::pair< point_2d<T>, point_2d<T> > > segm_vec;
- int num_runs = 100000;
- int num_segments = 3;
- point_vec.push_back(make_point_2d<T>(-100, -100));
- point_vec.push_back(make_point_2d<T>(-100, 100));
- point_vec.push_back(make_point_2d<T>(100, -100));
- point_vec.push_back(make_point_2d<T>(100, 100));
- for (int i = 0; i < num_runs; i++) {
- for (int j = 0; j < num_segments; j++) {
- T x1 = 0, y1 = 0, x2 = 0, y2 = 0;
- while (x1 == x2 && y1 == y2) {
- x1 = (rand() % 100) - 50;
- y1 = (rand() % 100) - 50;
- x2 = (rand() % 100) - 50;
- y2 = (rand() % 100) - 50;
- }
- point_2d<T> point1(x1, y1);
- point_2d<T> point2(x2, y2);
- segm_vec.push_back(std::make_pair(point1, point2));
- }
- remove_intersections(segm_vec);
- build_voronoi(point_vec, segm_vec, test_output);
- VERIFY_VORONOI_OUTPUT(test_output, NO_HALF_EDGE_INTERSECTIONS);
- segm_vec.clear();
- }
-}
-#endif
+//#ifdef NDEBUG
+//BOOST_AUTO_TEST_CASE_TEMPLATE(segment_random_test1, T, test_types) {
+// srand(static_cast<unsigned int>(time(NULL)));
+// voronoi_output<T> test_output;
+// std::vector< point_2d<T> > point_vec;
+// std::vector< std::pair< point_2d<T>, point_2d<T> > > segm_vec;
+// int num_runs = 10000;
+// int num_segments = 5;
+// point_vec.push_back(make_point_2d<T>(-100, -100));
+// point_vec.push_back(make_point_2d<T>(-100, 100));
+// point_vec.push_back(make_point_2d<T>(100, -100));
+// point_vec.push_back(make_point_2d<T>(100, 100));
+// for (int i = 0; i < num_runs; i++) {
+// for (int j = 0; j < num_segments; j++) {
+// T x1 = 0, y1 = 0, x2 = 0, y2 = 0;
+// while (x1 == x2 && y1 == y2) {
+// x1 = (rand() % 100) - 50;
+// y1 = (rand() % 100) - 50;
+// x2 = (rand() % 100) - 50;
+// y2 = (rand() % 100) - 50;
+// }
+// point_2d<T> point1(x1, y1);
+// point_2d<T> point2(x2, y2);
+// segm_vec.push_back(std::make_pair(point1, point2));
+// }
+// remove_intersections(segm_vec);
+// build_voronoi(point_vec, segm_vec, test_output);
+// VERIFY_VORONOI_OUTPUT(test_output, NO_HALF_EDGE_INTERSECTIONS);
+// segm_vec.clear();
+// }
+//}
+//#endif
-#ifdef NDEBUG
-BOOST_AUTO_TEST_CASE_TEMPLATE(segment_random_test2, T, test_types) {
- srand(static_cast<unsigned int>(time(NULL)));
- voronoi_output<T> test_output_small, test_output_large;
- std::vector< std::pair< point_2d<T>, point_2d<T> > > segm_vec;
- int num_segments[] = {10, 100, 1000, 10000};
- int num_runs[] = {1000, 100, 10, 1};
- int mod_koef1[] = {100, 1000, 10000, 100000};
- int mod_koef2[] = {100, 200, 300, 400};
- int max_value[] = {100, 600, 5150, 50200};
- int array_length = sizeof(num_segments) / sizeof(int);
- for (int k = 3; k < array_length; k++) {
- int koef = std::numeric_limits<int>::max() / max_value[k];
- for (int i = 0; i < num_runs[k]; i++) {
- for (int j = 0; j < num_segments[k]; j++) {
- T x1 = (rand() % (mod_koef1[k] / 100)) - mod_koef1[k] / 2;
- T y1 = (rand() % (mod_koef1[k] / 100)) - mod_koef1[k] / 2;
- T dx = 0, dy = 0;
- while (dx == 0 && dy == 0) {
- dx = (rand() % mod_koef2[k]) - mod_koef2[k] / 2;
- dy = (rand() % mod_koef2[k]) - mod_koef2[k] / 2;
- }
- T x2 = x1 + dx;
- T y2 = y1 + dy;
- point_2d<T> point1_small(x1, y1);
- point_2d<T> point2_small(x2, y2);
- segm_vec.push_back(std::make_pair(point1_small, point2_small));
- }
- remove_intersections(segm_vec);
- build_voronoi(segm_vec, test_output_small);
- for (size_t j = 0; j < segm_vec.size(); j++) {
- segm_vec[j].first.x(segm_vec[j].first.x() * koef);
- segm_vec[j].first.y(segm_vec[j].first.y() * koef);
- segm_vec[j].second.x(segm_vec[j].second.x() * koef);
- segm_vec[j].second.y(segm_vec[j].second.y() * koef);
- }
- build_voronoi(segm_vec, test_output_large);
- VERIFY_VORONOI_OUTPUT(test_output_small, NO_HALF_EDGE_INTERSECTIONS);
- BOOST_CHECK_EQUAL(test_output_small.num_cell_records(),
- test_output_large.num_cell_records());
- BOOST_CHECK_EQUAL(test_output_small.num_vertex_records(),
- test_output_large.num_vertex_records());
- BOOST_CHECK_EQUAL(test_output_small.num_edge_records(),
- test_output_large.num_edge_records());
- segm_vec.clear();
- }
- }
-}
-#endif
+//#ifdef NDEBUG
+//BOOST_AUTO_TEST_CASE_TEMPLATE(segment_random_test2, T, test_types) {
+// srand(static_cast<unsigned int>(time(NULL)));
+// voronoi_output<T> test_output_small, test_output_large;
+// std::vector< std::pair< point_2d<T>, point_2d<T> > > segm_vec;
+// int num_segments[] = {10, 100, 1000, 10000};
+// int num_runs[] = {1000, 100, 10, 1};
+// int mod_koef1[] = {100, 1000, 10000, 100000};
+// int mod_koef2[] = {100, 200, 300, 400};
+// int max_value[] = {100, 600, 5150, 50200};
+// int array_length = sizeof(num_segments) / sizeof(int);
+// for (int k = 0; k < 4; k++) {
+// int koef = std::numeric_limits<int>::max() / max_value[k];
+// for (int i = 0; i < num_runs[k]; i++) {
+// for (int j = 0; j < num_segments[k]; j++) {
+// T x1 = (rand() % (mod_koef1[k] / 100)) - mod_koef1[k] / 2;
+// T y1 = (rand() % (mod_koef1[k] / 100)) - mod_koef1[k] / 2;
+// T dx = 0, dy = 0;
+// while (dx == 0 && dy == 0) {
+// dx = (rand() % mod_koef2[k]) - mod_koef2[k] / 2;
+// dy = (rand() % mod_koef2[k]) - mod_koef2[k] / 2;
+// }
+// T x2 = x1 + dx;
+// T y2 = y1 + dy;
+// point_2d<T> point1_small(x1, y1);
+// point_2d<T> point2_small(x2, y2);
+// segm_vec.push_back(std::make_pair(point1_small, point2_small));
+// }
+// remove_intersections(segm_vec);
+// build_voronoi(segm_vec, test_output_small);
+// for (size_t j = 0; j < segm_vec.size(); j++) {
+// segm_vec[j].first.x(segm_vec[j].first.x() * koef);
+// segm_vec[j].first.y(segm_vec[j].first.y() * koef);
+// segm_vec[j].second.x(segm_vec[j].second.x() * koef);
+// segm_vec[j].second.y(segm_vec[j].second.y() * koef);
+// }
+// build_voronoi(segm_vec, test_output_large);
+// VERIFY_VORONOI_OUTPUT(test_output_large, NO_HALF_EDGE_INTERSECTIONS);
+// BOOST_CHECK_EQUAL(test_output_small.num_cell_records(),
+// test_output_large.num_cell_records());
+// BOOST_CHECK_EQUAL(test_output_small.num_vertex_records(),
+// test_output_large.num_vertex_records());
+// BOOST_CHECK_EQUAL(test_output_small.num_edge_records(),
+// test_output_large.num_edge_records());
+// segm_vec.clear();
+// }
+// }
+//}
+//#endif