Subject: [Boost-commit] svn:boost r77613 - sandbox/gtl/boost/polygon/detail
From: sydorchuk.andriy_at_[hidden]
Date: 2012-03-28 16:41:19
Author: asydorchuk
Date: 2012-03-28 16:41:18 EDT (Wed, 28 Mar 2012)
New Revision: 77613
URL: http://svn.boost.org/trac/boost/changeset/77613
Log:
Fixing voronoi_predicates for advanced Voronoi tutorial.
Text files modified:
sandbox/gtl/boost/polygon/detail/voronoi_predicates.hpp | 145 +++++++++++++++++++++++++++------------
1 files changed, 100 insertions(+), 45 deletions(-)
Modified: sandbox/gtl/boost/polygon/detail/voronoi_predicates.hpp
==============================================================================
--- sandbox/gtl/boost/polygon/detail/voronoi_predicates.hpp (original)
+++ sandbox/gtl/boost/polygon/detail/voronoi_predicates.hpp 2012-03-28 16:41:18 EDT (Wed, 28 Mar 2012)
@@ -49,8 +49,10 @@
// Compute robust cross_product: a1 * b2 - b1 * a2.
// 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_) {
+ static fpt_type robust_cross_product(int_x2_type a1_,
+ int_x2_type b1_,
+ int_x2_type a2_,
+ int_x2_type b2_) {
static to_fpt_converter to_fpt;
uint_x2_type a1 = static_cast<uint_x2_type>(is_neg(a1_) ? -a1_ : a1_);
uint_x2_type b1 = static_cast<uint_x2_type>(is_neg(b1_) ? -b1_ : b1_);
@@ -90,8 +92,10 @@
return (is_neg(value)) ? RIGHT : LEFT;
}
- template <typename T>
- static Orientation eval(T dif_x1_, T dif_y1_, T dif_x2_, T dif_y2_) {
+ static Orientation eval(int_x2_type dif_x1_,
+ int_x2_type dif_y1_,
+ int_x2_type dif_x2_,
+ int_x2_type dif_y2_) {
return eval(robust_cross_product(dif_x1_, dif_y1_, dif_x2_, dif_y2_));
}
@@ -306,8 +310,6 @@
const point_type &segment1 = site.point1(true);
fpt_type a1 = to_fpt(segment1.x()) - to_fpt(segment0.x());
fpt_type b1 = to_fpt(segment1.y()) - to_fpt(segment0.y());
- fpt_type a3 = to_fpt(point.x()) - to_fpt(segment0.x());
- fpt_type b3 = to_fpt(point.y()) - to_fpt(segment0.y());
fpt_type k = get_sqrt(a1 * a1 + b1 * b1);
// Avoid subtraction while computing k.
if (!is_neg(b1)) {
@@ -316,7 +318,11 @@
k = (k - b1) / (a1 * a1);
}
// The relative error is at most 7EPS.
- return robust_cross_product(a1, b1, a3, b3) * k;
+ return k * robust_cross_product(
+ static_cast<int_x2_type>(segment1.x()) - static_cast<int_x2_type>(segment0.x()),
+ static_cast<int_x2_type>(segment1.y()) - static_cast<int_x2_type>(segment0.y()),
+ static_cast<int_x2_type>(point.x()) - static_cast<int_x2_type>(segment0.x()),
+ static_cast<int_x2_type>(point.y()) - static_cast<int_x2_type>(segment0.y()));
}
}
@@ -343,7 +349,11 @@
return LESS;
return UNDEFINED;
} else {
- typename ot::Orientation orientation = ot::eval(a, b, dif_x, dif_y);
+ typename ot::Orientation orientation = ot::eval(
+ static_cast<int_x2_type>(segment_end.x()) - static_cast<int_x2_type>(segment_start.x()),
+ static_cast<int_x2_type>(segment_end.y()) - static_cast<int_x2_type>(segment_start.y()),
+ static_cast<int_x2_type>(new_point.x()) - static_cast<int_x2_type>(site_point.x()),
+ static_cast<int_x2_type>(new_point.y()) - static_cast<int_x2_type>(site_point.y()));
if (orientation == ot::LEFT) {
if (!right_site.is_inverse())
return reverse_order ? LESS : UNDEFINED;
@@ -715,7 +725,7 @@
a[0] * (static_cast<int_x2_type>(segm_start2.y()) -
static_cast<int_x2_type>(segm_start1.y()));
big_int_type dx = a[0] * (static_cast<int_x2_type>(site1.y()) -
- static_cast<int_x2_type>(segm_start1.y())) -
+ static_cast<int_x2_type>(segm_start1.y())) -
b[0] * (static_cast<int_x2_type>(site1.x()) -
static_cast<int_x2_type>(segm_start1.x()));
big_int_type dy = b[0] * (static_cast<int_x2_type>(site1.x()) -
@@ -989,8 +999,11 @@
fpt_type dif_x2 = to_fpt(site2.x()) - to_fpt(site3.x());
fpt_type dif_y1 = to_fpt(site1.y()) - to_fpt(site2.y());
fpt_type dif_y2 = to_fpt(site2.y()) - to_fpt(site3.y());
- fpt_type orientation =
- robust_cross_product(dif_x1, dif_y1, dif_x2, dif_y2);
+ fpt_type orientation = robust_cross_product(
+ static_cast<int_x2_type>(site1.x()) - static_cast<int_x2_type>(site2.x()),
+ static_cast<int_x2_type>(site2.x()) - static_cast<int_x2_type>(site3.x()),
+ static_cast<int_x2_type>(site1.y()) - static_cast<int_x2_type>(site2.y()),
+ static_cast<int_x2_type>(site2.y()) - static_cast<int_x2_type>(site3.y()));
robust_fpt_type inv_orientation(to_fpt(0.5) / orientation, 2.0);
fpt_type sum_x1 = to_fpt(site1.x()) + to_fpt(site2.x());
fpt_type sum_x2 = to_fpt(site2.x()) + to_fpt(site3.x());
@@ -1037,18 +1050,26 @@
to_fpt(site3.point1(true).x());
fpt_type vec_x = to_fpt(site2.y()) - to_fpt(site1.y());
fpt_type vec_y = to_fpt(site1.x()) - to_fpt(site2.x());
- robust_fpt_type teta(
- robust_cross_product(line_a, line_b, -vec_y, vec_x), 1.0);
+ robust_fpt_type teta(robust_cross_product(
+ static_cast<int_x2_type>(site3.point1(true).y()) - static_cast<int_x2_type>(site3.point0(true).y()),
+ static_cast<int_x2_type>(site3.point0(true).x()) - static_cast<int_x2_type>(site3.point1(true).x()),
+ static_cast<int_x2_type>(site2.x()) - static_cast<int_x2_type>(site1.x()),
+ static_cast<int_x2_type>(site2.y()) - static_cast<int_x2_type>(site1.y())), 1.0);
robust_fpt_type A(robust_cross_product(
- line_a, line_b,
- to_fpt(site3.point1().y()) - to_fpt(site1.y()),
- to_fpt(site1.x()) - to_fpt(site3.point1().x())), 1.0);
+ static_cast<int_x2_type>(site3.point1(true).y()) - static_cast<int_x2_type>(site3.point0(true).y()),
+ static_cast<int_x2_type>(site3.point0(true).x()) - static_cast<int_x2_type>(site3.point1(true).x()),
+ static_cast<int_x2_type>(site3.point1().y()) - static_cast<int_x2_type>(site1.y()),
+ static_cast<int_x2_type>(site1.x()) - static_cast<int_x2_type>(site3.point1().x())), 1.0);
robust_fpt_type B(robust_cross_product(
- line_a, line_b,
- to_fpt(site3.point1().y()) - to_fpt(site2.y()),
- to_fpt(site2.x()) - to_fpt(site3.point1().x())), 1.0);
- robust_fpt_type denom(
- robust_cross_product(vec_x, vec_y, line_a, line_b), 1.0);
+ static_cast<int_x2_type>(site3.point1(true).y()) - static_cast<int_x2_type>(site3.point0(true).y()),
+ static_cast<int_x2_type>(site3.point0(true).x()) - static_cast<int_x2_type>(site3.point1(true).x()),
+ static_cast<int_x2_type>(site3.point1().y()) - static_cast<int_x2_type>(site2.y()),
+ static_cast<int_x2_type>(site2.x()) - static_cast<int_x2_type>(site3.point1().x())), 1.0);
+ robust_fpt_type denom(robust_cross_product(
+ static_cast<int_x2_type>(site2.y()) - static_cast<int_x2_type>(site1.y()),
+ static_cast<int_x2_type>(site1.x()) - static_cast<int_x2_type>(site2.x()),
+ static_cast<int_x2_type>(site3.point1(true).y()) - static_cast<int_x2_type>(site3.point0(true).y()),
+ static_cast<int_x2_type>(site3.point0(true).x()) - static_cast<int_x2_type>(site3.point1(true).x())), 1.0);
robust_fpt_type inv_segm_len(to_fpt(1.0) /
get_sqrt(line_a * line_a + line_b * line_b), 3.0);
robust_dif_type t;
@@ -1105,22 +1126,29 @@
fpt_type a2 = to_fpt(segm_end2.x()) - to_fpt(segm_start2.x());
fpt_type b2 = to_fpt(segm_end2.y()) - to_fpt(segm_start2.y());
bool recompute_c_x, recompute_c_y, recompute_lower_x;
- robust_fpt_type orientation(robust_cross_product(b1, a1, b2, a2), 1.0);
+ robust_fpt_type orientation(robust_cross_product(
+ static_cast<int_x2_type>(segm_end1.y()) - static_cast<int_x2_type>(segm_start1.y()),
+ static_cast<int_x2_type>(segm_end1.x()) - static_cast<int_x2_type>(segm_start1.x()),
+ static_cast<int_x2_type>(segm_end2.y()) - static_cast<int_x2_type>(segm_start2.y()),
+ static_cast<int_x2_type>(segm_end2.x()) - static_cast<int_x2_type>(segm_start2.x())), 1.0);
if (ot::eval(orientation) == ot::COLLINEAR) {
robust_fpt_type a(a1 * a1 + b1 * b1, 2.0);
robust_fpt_type c(robust_cross_product(
- b1, a1,
- to_fpt(segm_start2.y()) - to_fpt(segm_start1.y()),
- to_fpt(segm_start2.x()) - to_fpt(segm_start1.x())), 1.0);
+ static_cast<int_x2_type>(segm_end1.y()) - static_cast<int_x2_type>(segm_start1.y()),
+ static_cast<int_x2_type>(segm_end1.x()) - static_cast<int_x2_type>(segm_start1.x()),
+ static_cast<int_x2_type>(segm_start2.y()) - static_cast<int_x2_type>(segm_start1.y()),
+ static_cast<int_x2_type>(segm_start2.x()) - static_cast<int_x2_type>(segm_start1.x())), 1.0);
robust_fpt_type det(
robust_cross_product(
- a1, b1,
- to_fpt(site1.x()) - to_fpt(segm_start1.x()),
- to_fpt(site1.y()) - to_fpt(segm_start1.y())) *
+ static_cast<int_x2_type>(segm_end1.x()) - static_cast<int_x2_type>(segm_start1.x()),
+ static_cast<int_x2_type>(segm_end1.y()) - static_cast<int_x2_type>(segm_start1.y()),
+ static_cast<int_x2_type>(site1.x()) - static_cast<int_x2_type>(segm_start1.x()),
+ static_cast<int_x2_type>(site1.y()) - static_cast<int_x2_type>(segm_start1.y())) *
robust_cross_product(
- b1, a1,
- to_fpt(site1.y()) - to_fpt(segm_start2.y()),
- to_fpt(site1.x()) - to_fpt(segm_start2.x())),
+ static_cast<int_x2_type>(segm_end1.y()) - static_cast<int_x2_type>(segm_start1.y()),
+ static_cast<int_x2_type>(segm_end1.x()) - static_cast<int_x2_type>(segm_start1.x()),
+ static_cast<int_x2_type>(site1.y()) - static_cast<int_x2_type>(segm_start2.y()),
+ static_cast<int_x2_type>(site1.x()) - static_cast<int_x2_type>(segm_start2.x())),
3.0);
robust_dif_type t;
t -= robust_fpt_type(a1, false) * robust_fpt_type((
@@ -1156,25 +1184,37 @@
} else {
robust_fpt_type sqr_sum1(get_sqrt(a1 * a1 + b1 * b1), 2.0);
robust_fpt_type sqr_sum2(get_sqrt(a2 * a2 + b2 * b2), 2.0);
- robust_fpt_type a(robust_cross_product(a1, b1, -b2, a2), 1.0);
+ robust_fpt_type a(robust_cross_product(
+ static_cast<int_x2_type>(segm_end1.x()) - static_cast<int_x2_type>(segm_start1.x()),
+ static_cast<int_x2_type>(segm_end1.y()) - static_cast<int_x2_type>(segm_start1.y()),
+ static_cast<int_x2_type>(segm_start2.y()) - static_cast<int_x2_type>(segm_end2.y()),
+ static_cast<int_x2_type>(segm_end2.x()) - static_cast<int_x2_type>(segm_start2.x())), 1.0);
if (!is_neg(a)) {
a += sqr_sum1 * sqr_sum2;
} else {
a = (orientation * orientation) / (sqr_sum1 * sqr_sum2 - a);
}
robust_fpt_type or1(robust_cross_product(
- b1, a1,
- to_fpt(segm_end1.y()) - to_fpt(site1.y()),
- to_fpt(segm_end1.x()) - to_fpt(site1.x())), 1.0);
+ static_cast<int_x2_type>(segm_end1.y()) - static_cast<int_x2_type>(segm_start1.y()),
+ static_cast<int_x2_type>(segm_end1.x()) - static_cast<int_x2_type>(segm_start1.x()),
+ static_cast<int_x2_type>(segm_end1.y()) - static_cast<int_x2_type>(site1.y()),
+ static_cast<int_x2_type>(segm_end1.x()) - static_cast<int_x2_type>(site1.x())), 1.0);
robust_fpt_type or2(robust_cross_product(
- a2, b2,
- to_fpt(segm_end2.x()) - to_fpt(site1.x()),
- to_fpt(segm_end2.y()) - to_fpt(site1.y())), 1.0);
+ static_cast<int_x2_type>(segm_end2.x()) - static_cast<int_x2_type>(segm_start2.x()),
+ static_cast<int_x2_type>(segm_end2.y()) - static_cast<int_x2_type>(segm_start2.y()),
+ static_cast<int_x2_type>(segm_end2.x()) - static_cast<int_x2_type>(site1.x()),
+ static_cast<int_x2_type>(segm_end2.y()) - static_cast<int_x2_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, to_fpt(segm_end1.y()), to_fpt(segm_end1.x())), 1.0);
+ static_cast<int_x2_type>(segm_end1.y()) - static_cast<int_x2_type>(segm_start1.y()),
+ static_cast<int_x2_type>(segm_end1.x()) - static_cast<int_x2_type>(segm_start1.x()),
+ static_cast<int_x2_type>(segm_end1.y()),
+ static_cast<int_x2_type>(segm_end1.x())), 1.0);
robust_fpt_type c2(robust_cross_product(
- a2, b2, to_fpt(segm_end2.x()), to_fpt(segm_end2.y())), 1.0);
+ static_cast<int_x2_type>(segm_end2.x()) - static_cast<int_x2_type>(segm_start2.x()),
+ static_cast<int_x2_type>(segm_end2.y()) - static_cast<int_x2_type>(segm_start2.y()),
+ static_cast<int_x2_type>(segm_end2.x()),
+ static_cast<int_x2_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;
@@ -1188,9 +1228,15 @@
b += iy * (robust_fpt_type(b1, false) * sqr_sum2);
b += iy * (robust_fpt_type(b2, false) * sqr_sum1);
b -= sqr_sum1 * robust_fpt_type(robust_cross_product(
- a2, b2, to_fpt(-site1.y()), to_fpt(site1.x())), 1.0);
+ static_cast<int_x2_type>(segm_end2.x()) - static_cast<int_x2_type>(segm_start2.x()),
+ static_cast<int_x2_type>(segm_end2.y()) - static_cast<int_x2_type>(segm_start2.y()),
+ static_cast<int_x2_type>(-site1.y()),
+ static_cast<int_x2_type>(site1.x())), 1.0);
b -= sqr_sum2 * robust_fpt_type(robust_cross_product(
- a1, b1, to_fpt(-site1.y()), to_fpt(site1.x())), 1.0);
+ static_cast<int_x2_type>(segm_end1.x()) - static_cast<int_x2_type>(segm_start1.x()),
+ static_cast<int_x2_type>(segm_end1.y()) - static_cast<int_x2_type>(segm_start1.y()),
+ static_cast<int_x2_type>(-site1.y()),
+ static_cast<int_x2_type>(site1.x())), 1.0);
t -= b;
if (point_index == 2) {
t += det.sqrt();
@@ -1251,11 +1297,20 @@
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()), 1.0);
+ static_cast<int_x2_type>(site1.x1(true)) - static_cast<int_x2_type>(site1.x0(true)),
+ static_cast<int_x2_type>(site1.y1(true)) - static_cast<int_x2_type>(site1.y0(true)),
+ static_cast<int_x2_type>(site2.x1(true)) - static_cast<int_x2_type>(site2.x0(true)),
+ static_cast<int_x2_type>(site2.y1(true)) - static_cast<int_x2_type>(site2.y0(true))), 1.0);
robust_fpt_type cross_23(robust_cross_product(
- a2.fpv(), b2.fpv(), a3.fpv(), b3.fpv()), 1.0);
+ static_cast<int_x2_type>(site2.x1(true)) - static_cast<int_x2_type>(site2.x0(true)),
+ static_cast<int_x2_type>(site2.y1(true)) - static_cast<int_x2_type>(site2.y0(true)),
+ static_cast<int_x2_type>(site3.x1(true)) - static_cast<int_x2_type>(site3.x0(true)),
+ static_cast<int_x2_type>(site3.y1(true)) - static_cast<int_x2_type>(site3.y0(true))), 1.0);
robust_fpt_type cross_31(robust_cross_product(
- a3.fpv(), b3.fpv(), a1.fpv(), b1.fpv()), 1.0);
+ static_cast<int_x2_type>(site3.x1(true)) - static_cast<int_x2_type>(site3.x0(true)),
+ static_cast<int_x2_type>(site3.y1(true)) - static_cast<int_x2_type>(site3.y0(true)),
+ static_cast<int_x2_type>(site1.x1(true)) - static_cast<int_x2_type>(site1.x0(true)),
+ static_cast<int_x2_type>(site1.y1(true)) - static_cast<int_x2_type>(site1.y0(true))), 1.0);
robust_dif_type denom, c_x, c_y, r;
// denom = cross_12 * len3 + cross_23 * len1 + cross_31 * len2.