Boost logo

Boost-Commit :

Subject: [Boost-commit] svn:boost r70236 - sandbox/SOC/2010/sweepline/boost/sweepline/detail
From: sydorchuk.andriy_at_[hidden]
Date: 2011-03-20 19:14:25


Author: asydorchuk
Date: 2011-03-20 19:14:24 EDT (Sun, 20 Mar 2011)
New Revision: 70236
URL: http://svn.boost.org/trac/boost/changeset/70236

Log:
Added (point, point, segment) predicate implementation using robust sqrt evaluator class.
Text files modified:
   sandbox/SOC/2010/sweepline/boost/sweepline/detail/mpz_arithmetic.hpp | 12 ++
   sandbox/SOC/2010/sweepline/boost/sweepline/detail/voronoi_formation.hpp | 223 ++++++++++++++++++++++++++++-----------
   2 files changed, 169 insertions(+), 66 deletions(-)

Modified: sandbox/SOC/2010/sweepline/boost/sweepline/detail/mpz_arithmetic.hpp
==============================================================================
--- sandbox/SOC/2010/sweepline/boost/sweepline/detail/mpz_arithmetic.hpp (original)
+++ sandbox/SOC/2010/sweepline/boost/sweepline/detail/mpz_arithmetic.hpp 2011-03-20 19:14:24 EDT (Sun, 20 Mar 2011)
@@ -103,6 +103,12 @@
                         return temp_[next_cur()];
                 }
 
+ mpz_wrapper& operator*(double that) const {
+ temp_[cur_].m_ = that;
+ temp_[cur_].m_ *= this->m_;
+ return temp_[next_cur()];
+ }
+
                 mpz_wrapper& operator+=(const mpz_wrapper& that) {
                         this->m_ += that.m_;
                         return *this;
@@ -181,6 +187,8 @@
                         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 (ret_val == 0.0)
+ return 0.0;
                         if (!A[1].is_negative())
                                 return ret_val;
                         if (!A[0].is_positive())
@@ -200,6 +208,8 @@
                         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 (res == 0.0)
+ return res;
                         if (!A[2].is_negative())
                                 return res;
                         if (!A[0].is_positive())
@@ -234,6 +244,8 @@
                                                  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 (res == 0.0)
+ return res;
                         if (!A[3].is_negative())
                                 return res;
                         if (!A[0].is_positive())

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-20 19:14:24 EDT (Sun, 20 Mar 2011)
@@ -1295,7 +1295,8 @@
     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;
@@ -1316,8 +1317,7 @@
                 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]);
+ double orientation = (mpz_dif_x[0] * mpz_dif_y[1] - mpz_dif_x[1] * mpz_dif_y[0]).get_d();
 
                 // If use this function only to recompute parameters of the circle
                 // event, this check won't be required.
@@ -1346,11 +1346,11 @@
 
                 // 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;
+ double c_x = mpz_c_x.get_d() * inv_orientation;
+ double c_y = mpz_c_y.get_d() * inv_orientation;
 
                 // r(r) = 1.5 * EPS < 2 * EPS.
- double r = sqrt(static_cast<double>(mpz_sqr_r));
+ double r = sqrt(mpz_sqr_r.get_d());
 
                 // If c_x >= 0 then lower_x = c_x + r,
                 // else lower_x = (c_x * c_x - r * r) / (c_x - r).
@@ -1361,13 +1361,157 @@
                 } 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);
+ double lower_x = mpz_numerator[0].get_d() * fabs(inv_orientation);
+ lower_x /= (!mpz_c_x.is_negative()) ? (-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;
         }
-*/
+
+ template <typename T>
+ static bool create_circle_event_pps_gmpxx(const site_event<T> &site1,
+ const site_event<T> &site2,
+ const site_event<T> &site3,
+ int segment_index,
+ circle_event<T> &c_event) {
+ // TODO(asydorchuk): Checks are not required if we recompute event parameters.
+ if (segment_index != 2) {
+ kOrientation orient1 = orientation_test(site1.point0(),
+ site2.point0(), site3.point0(true));
+ kOrientation orient2 = orientation_test(site1.point0(),
+ site2.point0(), site3.point1(true));
+ if (segment_index == 1 && site1.x0() >= site2.x0()) {
+ if (orient1 != RIGHT_ORIENTATION)
+ return false;
+ } else if (segment_index == 3 && site2.x0() >= site1.x0()) {
+ if (orient2 != RIGHT_ORIENTATION)
+ return false;
+ } else if (orient1 != RIGHT_ORIENTATION && orient2 != RIGHT_ORIENTATION) {
+ return false;
+ }
+ } else {
+ if (site3.point0(true) == site1.point0() &&
+ site3.point1(true) == site2.point0())
+ return false;
+ }
+
+ static mpz_wrapper line_a, line_b, segm_len, vec_x, vec_y, sum_x, sum_y, teta,
+ denom, A, B, sum_AB, dif[2], numer, cA[4], cB[4], det;
+
+ line_a = site3.point1(true).y() - site3.point0(true).y();
+ line_b = site3.point0(true).x() - site3.point1(true).x();
+ segm_len = line_a * line_a + line_b * line_b;
+ vec_x = site2.y() - site1.y();
+ vec_y = site1.x() - site2.x();
+ sum_x = site1.x() + site2.x();
+ sum_y = site1.y() + site2.y();
+ teta = line_a * vec_x + line_b * vec_y;
+ denom = vec_x * line_b - vec_y * line_a;
+
+ dif[0] = site3.point1().y() - site1.y();
+ dif[1] = site1.x() - site3.point1().x();
+ A = line_a * dif[1] - line_b * dif[0];
+ dif[0] = site3.point1().y() - site2.y();
+ dif[1] = site2.x() - site3.point1().x();
+ B = line_a * dif[1] - line_b * dif[0];
+ sum_AB = A + B;
+
+ if (!denom.is_negative() && !denom.is_positive()) {
+ numer = teta * teta - sum_AB * sum_AB;
+ denom = teta * sum_AB;
+ cA[0] = denom * sum_x * 2.0 + numer * vec_x;
+ cB[0] = segm_len;
+ cA[1] = denom * sum_AB * 2.0 + numer * teta;
+ cB[1] = 1;
+ cA[2] = denom * sum_y * 2.0 + numer * vec_y;
+ double c_x = 0.25 * cA[0].get_d() / denom.get_d();
+ double c_y = 0.25 * cA[2].get_d() / denom.get_d();
+ double lower_x = 0.25 * sqr_expr_evaluator<2>::eval(cA, cB) /
+ (denom.get_d() * sqrt(segm_len.get_d()));
+ c_event = circle_event<double>(c_x, c_y, lower_x);
+ return true;
+ }
+
+ det = (teta * teta + denom * denom) * A * B * 4.0;
+ cA[0] = sum_x * denom * denom + teta * sum_AB * vec_x;
+ cB[0] = 1.0;
+ cA[1] = (segment_index == 2) ? -vec_x : vec_x;
+ cB[1] = det;
+ double c_x = 0.5 * sqr_expr_evaluator<2>::eval(cA, cB) / (denom.get_d() * denom.get_d());
+
+ cA[2] = sum_y * denom * denom + teta * sum_AB * vec_y;
+ cB[2] = 1.0;
+ cA[3] = (segment_index == 2) ? -vec_y : vec_y;
+ cB[3] = det;
+ double c_y = 0.5 * sqr_expr_evaluator<2>::eval(&cA[2], &cB[2]) / (denom.get_d() * denom.get_d());
+
+ cB[0] *= segm_len;
+ cB[1] *= segm_len;
+ cA[2] = sum_AB * (denom * denom + teta * teta);
+ cB[2] = 1.0;
+ cA[3] = (segment_index == 2) ? -teta : teta;
+ cB[3] = det;
+ double lower_x = 0.5 * sqr_expr_evaluator<4>::eval(cA, cB) /
+ (denom.get_d() * denom.get_d() * sqrt(segm_len.get_d()));
+
+ c_event = circle_event<double>(c_x, c_y, lower_x);
+ 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
     // point sites.
     template <typename T>
@@ -1406,7 +1550,7 @@
         c_event = circle_event<double>(c_x, c_y, lower_x);
         return true;
     }
-
+
     // Find parameters of the inscribed circle that is tangent to two
     // point sites and on segment site.
     template <typename T>
@@ -1435,8 +1579,8 @@
                 return false;
         }
 
- double line_a = site3.point1().y() - site3.point0().y();
- double line_b = site3.point0().x() - site3.point1().x();
+ double line_a = site3.point1(true).y() - site3.point0(true).y();
+ double line_b = site3.point0(true).x() - site3.point1(true).x();
         double vec_x = site2.y() - site1.y();
         double vec_y = site1.x() - site2.x();
         double teta = robust_cross_product(line_a, line_b, -vec_y, vec_x);
@@ -1471,7 +1615,7 @@
         r -= line_b * site3.y0();
         r += line_a * c_x;
         r += line_b * c_y;
- r.abs();
+ r.abs();
         lower_x += r * inv_segm_len;
         c_event = circle_event<double>(c_x, c_y, lower_x);
         return true;
@@ -1581,60 +1725,7 @@
         }
         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>


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