|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r76148 - sandbox/gtl/boost/polygon/detail
From: sydorchuk.andriy_at_[hidden]
Date: 2011-12-24 18:06:48
Author: asydorchuk
Date: 2011-12-24 18:06:47 EST (Sat, 24 Dec 2011)
New Revision: 76148
URL: http://svn.boost.org/trac/boost/changeset/76148
Log:
It was math proven that relative error for the robust cross product is atmost 1 EPS.
Text files modified:
sandbox/gtl/boost/polygon/detail/voronoi_calc_utils.hpp | 59 +++++++++++++++++++--------------------
sandbox/gtl/boost/polygon/detail/voronoi_structures.hpp | 8 ++--
2 files changed, 33 insertions(+), 34 deletions(-)
Modified: sandbox/gtl/boost/polygon/detail/voronoi_calc_utils.hpp
==============================================================================
--- sandbox/gtl/boost/polygon/detail/voronoi_calc_utils.hpp (original)
+++ sandbox/gtl/boost/polygon/detail/voronoi_calc_utils.hpp 2011-12-24 18:06:47 EST (Sat, 24 Dec 2011)
@@ -19,7 +19,7 @@
template <typename T>
class voronoi_calc_utils;
-// Predicate utils. Operates with the coordinate types that could
+// Predicate utilities. Operates with the coordinate types that could
// be converted to the 32-bit signed integer without precision loss.
template <>
class voronoi_calc_utils<int32> {
@@ -55,7 +55,8 @@
}
// Compute robust cross_product: a1 * b2 - b1 * a2.
- // The result is correct with epsilon relative error equal to 1EPS.
+ // It was mathematically proven that the result is correct
+ // with epsilon relative error equal to 1EPS.
template <typename T>
static fpt_type robust_cross_product(T a1_, T b1_, T a2_, T b2_) {
uint64 a1, b1, a2, b2;
@@ -78,15 +79,13 @@
if (!expr_l_plus) {
if (expr_r_plus)
- return -static_cast<fpt_type>(expr_l) -
- static_cast<fpt_type>(expr_r);
+ return -static_cast<fpt_type>(expr_l - expr_r);
else return (expr_l > expr_r) ?
-static_cast<fpt_type>(expr_l - expr_r) :
static_cast<fpt_type>(expr_r - expr_l);
} else {
if (!expr_r_plus)
- return static_cast<fpt_type>(expr_l) +
- static_cast<fpt_type>(expr_r);
+ return static_cast<fpt_type>(expr_l + expr_r);
else
return (expr_l < expr_r) ?
-static_cast<fpt_type>(expr_r - expr_l) :
@@ -282,7 +281,7 @@
fpt_type dist1 = find_distance_to_point_arc(left_site, new_site.point0());
fpt_type dist2 = find_distance_to_segment_arc(right_site, new_site.point0());
- // The undefined ulp range is equal to 3EPS + 8EPS <= 11ULP.
+ // The undefined ulp range is equal to 3EPS + 7EPS <= 10ULP.
return reverse_order ^ (dist1 < dist2);
}
@@ -300,7 +299,7 @@
fpt_type dist1 = find_distance_to_segment_arc(left_site, new_site.point0());
fpt_type dist2 = find_distance_to_segment_arc(right_site, new_site.point0());
- // The undefined ulp range is equal to 8EPS + 8EPS <= 16ULP.
+ // The undefined ulp range is equal to 7EPS + 7EPS <= 14ULP.
return dist1 < dist2;
}
@@ -335,7 +334,7 @@
} else {
k = (k - b1) / (a1 * a1);
}
- // The relative error is atmost 8EPS.
+ // The relative error is atmost 7EPS.
return robust_cross_product(a1, b1, a3, b3) * k;
}
}
@@ -984,7 +983,7 @@
fpt_type dif_y2 = static_cast<fpt_type>(site2.y()) -
static_cast<fpt_type>(site3.y());
fpt_type orientation = robust_cross_product(dif_x1, dif_y1, dif_x2, dif_y2);
- robust_fpt_type inv_orientation(0.5 / orientation, 3.0);
+ robust_fpt_type inv_orientation(0.5 / orientation, 2.0);
fpt_type sum_x1 = static_cast<fpt_type>(site1.x()) +
static_cast<fpt_type>(site2.x());
fpt_type sum_x2 = static_cast<fpt_type>(site2.x()) +
@@ -1035,18 +1034,18 @@
static_cast<fpt_type>(site1.y());
fpt_type vec_y = static_cast<fpt_type>(site1.x()) -
static_cast<fpt_type>(site2.x());
- robust_fpt_type teta(robust_cross_product(line_a, line_b, -vec_y, vec_x), 2.0);
+ robust_fpt_type teta(robust_cross_product(line_a, line_b, -vec_y, vec_x), 1.0);
robust_fpt_type A(robust_cross_product(line_a, line_b,
static_cast<fpt_type>(site3.point1().y()) -
static_cast<fpt_type>(site1.y()),
static_cast<fpt_type>(site1.x()) -
- static_cast<fpt_type>(site3.point1().x())), 2.0);
+ static_cast<fpt_type>(site3.point1().x())), 1.0);
robust_fpt_type B(robust_cross_product(line_a, line_b,
static_cast<fpt_type>(site3.point1().y()) -
static_cast<fpt_type>(site2.y()),
static_cast<fpt_type>(site2.x()) -
- static_cast<fpt_type>(site3.point1().x())), 2.0);
- robust_fpt_type denom(robust_cross_product(vec_x, vec_y, line_a, line_b), 2.0);
+ static_cast<fpt_type>(site3.point1().x())), 1.0);
+ robust_fpt_type denom(robust_cross_product(vec_x, vec_y, line_a, line_b), 1.0);
robust_fpt_type inv_segm_len(1.0 / std::sqrt(sqr_distance(line_a, line_b)), 3.0);
robust_dif_type t;
if (get_orientation(denom) == COLLINEAR) {
@@ -1106,14 +1105,14 @@
fpt_type b2 = static_cast<fpt_type>(segm_end2.y()) -
static_cast<fpt_type>(segm_start2.y());
bool recompute_c_x, recompute_c_y, recompute_lower_x;
- robust_fpt_type orientation(robust_cross_product(b1, a1, b2, a2), 2.0);
+ robust_fpt_type orientation(robust_cross_product(b1, a1, b2, a2), 1.0);
if (get_orientation(orientation) == COLLINEAR) {
robust_fpt_type a(a1 * a1 + b1 * b1, 2.0);
robust_fpt_type c(robust_cross_product(b1, a1,
static_cast<fpt_type>(segm_start2.y()) -
static_cast<fpt_type>(segm_start1.y()),
static_cast<fpt_type>(segm_start2.x()) -
- static_cast<fpt_type>(segm_start1.x())), 2.0);
+ static_cast<fpt_type>(segm_start1.x())), 1.0);
robust_fpt_type det(robust_cross_product(a1, b1,
static_cast<fpt_type>(site1.x()) -
static_cast<fpt_type>(segm_start1.x()),
@@ -1123,7 +1122,7 @@
static_cast<fpt_type>(site1.y()) -
static_cast<fpt_type>(segm_start2.y()),
static_cast<fpt_type>(site1.x()) -
- static_cast<fpt_type>(segm_start2.x())), 5.0);
+ static_cast<fpt_type>(segm_start2.x())), 3.0);
robust_dif_type t;
t -= robust_fpt_type(a1, false) * robust_fpt_type(
(static_cast<fpt_type>(segm_start1.x()) +
@@ -1157,7 +1156,7 @@
} else {
robust_fpt_type sqr_sum1(std::sqrt(a1 * a1 + b1 * b1), 2.0);
robust_fpt_type sqr_sum2(std::sqrt(a2 * a2 + b2 * b2), 2.0);
- robust_fpt_type a(robust_cross_product(a1, b1, -b2, a2), 2.0);
+ robust_fpt_type a(robust_cross_product(a1, b1, -b2, a2), 1.0);
if (a >= 0) {
a += sqr_sum1 * sqr_sum2;
} else {
@@ -1167,19 +1166,19 @@
static_cast<fpt_type>(segm_end1.y()) -
static_cast<fpt_type>(site1.y()),
static_cast<fpt_type>(segm_end1.x()) -
- static_cast<fpt_type>(site1.x())), 2.0);
+ static_cast<fpt_type>(site1.x())), 1.0);
robust_fpt_type or2(robust_cross_product(a2, b2,
static_cast<fpt_type>(segm_end2.x()) -
static_cast<fpt_type>(site1.x()),
static_cast<fpt_type>(segm_end2.y()) -
- static_cast<fpt_type>(site1.y())), 2.0);
+ static_cast<fpt_type>(site1.y())), 1.0);
robust_fpt_type det = robust_fpt_type(2.0, false) * a * or1 * or2;
robust_fpt_type c1(robust_cross_product(b1, a1,
static_cast<fpt_type>(segm_end1.y()),
- static_cast<fpt_type>(segm_end1.x())), 2.0);
+ static_cast<fpt_type>(segm_end1.x())), 1.0);
robust_fpt_type c2(robust_cross_product(a2, b2,
static_cast<fpt_type>(segm_end2.x()),
- static_cast<fpt_type>(segm_end2.y())), 2.0);
+ static_cast<fpt_type>(segm_end2.y())), 1.0);
robust_fpt_type inv_orientation = robust_fpt_type(1.0, false) / orientation;
robust_dif_type t, b, ix, iy;
ix += robust_fpt_type(a2, false) * c1 * inv_orientation;
@@ -1193,10 +1192,10 @@
b += iy * (robust_fpt_type(b2, false) * sqr_sum1);
b -= sqr_sum1 * robust_fpt_type(robust_cross_product(a2, b2,
static_cast<fpt_type>(-site1.y()),
- static_cast<fpt_type>(site1.x())), 2.0);
+ static_cast<fpt_type>(site1.x())), 1.0);
b -= sqr_sum2 * robust_fpt_type(robust_cross_product(a1, b1,
static_cast<fpt_type>(-site1.y()),
- static_cast<fpt_type>(site1.x())), 2.0);
+ static_cast<fpt_type>(site1.x())), 1.0);
t -= b;
if (point_index == 2) {
t += det.sqrt();
@@ -1232,26 +1231,26 @@
static_cast<fpt_type>(site1.x0(true)), 0.0);
robust_fpt_type b1(static_cast<fpt_type>(site1.y1(true)) -
static_cast<fpt_type>(site1.y0(true)), 0.0);
- robust_fpt_type c1(robust_cross_product(site1.x0(true), site1.y0(true), site1.x1(true), site1.y1(true)), 2.0);
+ robust_fpt_type c1(robust_cross_product(site1.x0(true), site1.y0(true), site1.x1(true), site1.y1(true)), 1.0);
robust_fpt_type a2(static_cast<fpt_type>(site2.x1(true)) -
static_cast<fpt_type>(site2.x0(true)), 0.0);
robust_fpt_type b2(static_cast<fpt_type>(site2.y1(true)) -
static_cast<fpt_type>(site2.y0(true)), 0.0);
- robust_fpt_type c2(robust_cross_product(site2.x0(true), site2.y0(true), site2.x1(true), site2.y1(true)), 2.0);
+ robust_fpt_type c2(robust_cross_product(site2.x0(true), site2.y0(true), site2.x1(true), site2.y1(true)), 1.0);
robust_fpt_type a3(static_cast<fpt_type>(site3.x1(true)) -
static_cast<fpt_type>(site3.x0(true)), 0.0);
robust_fpt_type b3(static_cast<fpt_type>(site3.y1(true)) -
static_cast<fpt_type>(site3.y0(true)), 0.0);
- robust_fpt_type c3(robust_cross_product(site3.x0(true), site3.y0(true), site3.x1(true), site3.y1(true)), 2.0);
+ robust_fpt_type c3(robust_cross_product(site3.x0(true), site3.y0(true), site3.x1(true), site3.y1(true)), 1.0);
robust_fpt_type len1 = (a1 * a1 + b1 * b1).sqrt();
robust_fpt_type len2 = (a2 * a2 + b2 * b2).sqrt();
robust_fpt_type len3 = (a3 * a3 + b3 * b3).sqrt();
- robust_fpt_type cross_12(robust_cross_product(a1.fpv(), b1.fpv(), a2.fpv(), b2.fpv()), 2.0);
- robust_fpt_type cross_23(robust_cross_product(a2.fpv(), b2.fpv(), a3.fpv(), b3.fpv()), 2.0);
- robust_fpt_type cross_31(robust_cross_product(a3.fpv(), b3.fpv(), a1.fpv(), b1.fpv()), 2.0);
+ robust_fpt_type cross_12(robust_cross_product(a1.fpv(), b1.fpv(), a2.fpv(), b2.fpv()), 1.0);
+ robust_fpt_type cross_23(robust_cross_product(a2.fpv(), b2.fpv(), a3.fpv(), b3.fpv()), 1.0);
+ robust_fpt_type cross_31(robust_cross_product(a3.fpv(), b3.fpv(), a1.fpv(), b1.fpv()), 1.0);
robust_dif_type denom, c_x, c_y, r;
// denom = cross_12 * len3 + cross_23 * len1 + cross_31 * len2.
Modified: sandbox/gtl/boost/polygon/detail/voronoi_structures.hpp
==============================================================================
--- sandbox/gtl/boost/polygon/detail/voronoi_structures.hpp (original)
+++ sandbox/gtl/boost/polygon/detail/voronoi_structures.hpp 2011-12-24 18:06:47 EST (Sat, 24 Dec 2011)
@@ -214,9 +214,9 @@
// Circle event is made by the two consequtive nodes in the beach line data
// structure. In case another node was inserted during algorithm execution
// between the given two nodes circle event becomes inactive.
- // Variables: center_x_ - numerator of the center's x-coordinate.
- // center_y_ - numerator of the center's y-coordinate.
- // lower_x_ - numerator of the bottom point's x-coordinate.
+ // Variables: center_x_ - center x-coordinate.
+ // center_y_ - center y-coordinate.
+ // lower_x_ - leftmost x-coordinate.
// is_active_ - states whether circle event is still active.
// NOTE: lower_y coordinate is always equal to center_y.
template <typename T>
@@ -339,7 +339,7 @@
// the site and from the sweepline. If the site is a point then arc is
// a parabola, otherwise it's a line segment. A segment site event will
// produce different bisectors based on its direction.
- // In the general case two sites will create two opposite bisectors. That's
+ // In general case two sites will create two opposite bisectors. That's
// why the order of the sites is important to define the unique bisector.
// The one site is considered to be newer than the other one if it was
// processed by the algorithm later (has greater index).
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