Subject: [Boost-commit] svn:boost r75024 - in sandbox/gtl: boost/polygon/detail libs/polygon/test
From: sydorchuk.andriy_at_[hidden]
Date: 2011-10-17 18:17:06
Author: asydorchuk
Date: 2011-10-17 18:17:05 EDT (Mon, 17 Oct 2011)
New Revision: 75024
URL: http://svn.boost.org/trac/boost/changeset/75024
Log:
Integrated circle formation predicate tests.
Refactored circle formation predicates.
Unified the ULPS uaage.
Text files modified:
sandbox/gtl/boost/polygon/detail/voronoi_calc_kernel.hpp | 593 +++++++++++++++++++++++----------------
sandbox/gtl/libs/polygon/test/Jamfile.v2 | 2
sandbox/gtl/libs/polygon/test/voronoi_calc_kernel_test.cpp | 109 +++++++
3 files changed, 452 insertions(+), 252 deletions(-)
Modified: sandbox/gtl/boost/polygon/detail/voronoi_calc_kernel.hpp
==============================================================================
--- sandbox/gtl/boost/polygon/detail/voronoi_calc_kernel.hpp (original)
+++ sandbox/gtl/boost/polygon/detail/voronoi_calc_kernel.hpp 2011-10-17 18:17:05 EDT (Mon, 17 Oct 2011)
@@ -66,6 +66,11 @@
public:
typedef double fpt_type;
typedef unsigned long long ulong_long_type;
+
+ static const unsigned int ULPS = 64;
+ static const unsigned int ULPSx2 = (ULPS << 1);
+ static const fpt_type fULPS;
+ static const fpt_type fULPSx2;
enum kOrientation {
RIGHT = -1,
@@ -160,9 +165,6 @@
typedef Site site_type;
typedef Circle circle_type;
- static const unsigned int ULPS = 64;
- static const unsigned int ULPSx2 = (ULPS << 1);
-
bool operator()(const site_type &lhs, const site_type &rhs) const {
if (lhs.x0() != rhs.x0()) return lhs.x0() < rhs.x0();
if (!lhs.is_segment()) {
@@ -286,8 +288,7 @@
// result without any further computations.
return static_cast<fpt_type>(left_point.y()) +
static_cast<fpt_type>(right_point.y()) <
- static_cast<fpt_type>(2) *
- static_cast<fpt_type>(new_point.y());
+ 2.0 * static_cast<fpt_type>(new_point.y());
}
fpt_type dist1 = find_distance_to_point_arc(left_site, new_point);
@@ -374,8 +375,8 @@
static_cast<fpt_type>(segment0.y());
fpt_type k = std::sqrt(a1 * a1 + b1 * b1);
// Avoid substraction while computing k.
- if (b1 >= static_cast<fpt_type>(0)) {
- k = static_cast<fpt_type>(1) / (b1 + k);
+ if (b1 >= static_cast<fpt_type>(0.0)) {
+ k = static_cast<fpt_type>(1.0) / (b1 + k);
} else {
k = (k - b1) / (a1 * a1);
}
@@ -583,27 +584,36 @@
typedef mpt_wrapper<mpf_class, 2> mpf_type;
typedef robust_sqrt_expr<mpt_type, mpf_type> robust_sqrt_expr_type;
- bool ppp(const site_type &site1,
+ void ppp(const site_type &site1,
const site_type &site2,
const site_type &site3,
circle_type &circle,
bool recompute_c_x = true,
bool recompute_c_y = true,
bool recompute_lower_x = true) {
- typedef mpt_wrapper<mpz_class, 8> mpt_type;
static mpt_type mpz_dif_x[3], mpz_dif_y[3], mpz_sum_x[2], mpz_sum_y[2],
mpz_numerator[3], mpz_c_x, mpz_c_y, mpz_sqr_r, denom,
cA[2], cB[2];
- 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();
- 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_dif_x[0] = static_cast<fpt_type>(site1.x()) -
+ static_cast<fpt_type>(site2.x());
+ mpz_dif_x[1] = static_cast<fpt_type>(site2.x()) -
+ static_cast<fpt_type>(site3.x());
+ mpz_dif_x[2] = static_cast<fpt_type>(site1.x()) -
+ static_cast<fpt_type>(site3.x());
+ mpz_dif_y[0] = static_cast<fpt_type>(site1.y()) -
+ static_cast<fpt_type>(site2.y());
+ mpz_dif_y[1] = static_cast<fpt_type>(site2.y()) -
+ static_cast<fpt_type>(site3.y());
+ mpz_dif_y[2] = static_cast<fpt_type>(site1.y()) -
+ static_cast<fpt_type>(site3.y());
+ mpz_sum_x[0] = static_cast<fpt_type>(site1.x()) +
+ static_cast<fpt_type>(site2.x());
+ mpz_sum_x[1] = static_cast<fpt_type>(site2.x()) +
+ static_cast<fpt_type>(site3.x());
+ mpz_sum_y[0] = static_cast<fpt_type>(site1.y()) +
+ static_cast<fpt_type>(site2.y());
+ mpz_sum_y[1] = static_cast<fpt_type>(site2.y()) +
+ static_cast<fpt_type>(site3.y());
denom = (mpz_dif_x[0] * mpz_dif_y[1] - mpz_dif_x[1] * mpz_dif_y[0]) * 2.0;
mpz_numerator[1] = mpz_dif_x[0] * mpz_sum_x[0] + mpz_dif_y[0] * mpz_sum_y[0];
@@ -618,17 +628,17 @@
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];
- double r = std::sqrt(mpz_sqr_r.get_d());
+ fpt_type r = std::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).
// To guarantee epsilon relative error.
- if (circle.x() >= 0) {
+ if (circle.x() >= static_cast<fpt_type>(0.0)) {
circle.lower_x(circle.x() + r / fabs(denom.get_d()));
} else {
mpz_numerator[0] = mpz_c_x * mpz_c_x - mpz_sqr_r;
- double lower_x = mpz_numerator[0].get_d() /
- (denom.get_d() * (mpz_c_x.get_d() + r));
+ fpt_type lower_x = mpz_numerator[0].get_d() /
+ (denom.get_d() * (mpz_c_x.get_d() + r));
circle.lower_x(lower_x);
}
}
@@ -638,12 +648,10 @@
mpz_c_y = mpz_numerator[2] * mpz_dif_x[0] - mpz_numerator[1] * mpz_dif_x[1];
circle.y(mpz_c_y.get_d() / denom.get_d());
}
-
- return true;
}
// Recompute parameters of the circle event using high-precision library.
- bool pps(const site_type &site1,
+ void pps(const site_type &site1,
const site_type &site2,
const site_type &site3,
int segment_index,
@@ -651,26 +659,33 @@
bool recompute_c_x = true,
bool recompute_c_y = true,
bool recompute_lower_x = true) {
- typedef mpt_wrapper<mpz_class, 8> mpt_type;
- typedef mpt_wrapper<mpf_class, 2> mpf_type;
static mpt_type 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();
+ line_a = static_cast<fpt_type>(site3.point1(true).y()) -
+ static_cast<fpt_type>(site3.point0(true).y());
+ line_b = static_cast<fpt_type>(site3.point0(true).x()) -
+ static_cast<fpt_type>(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();
+ vec_x = static_cast<fpt_type>(site2.y()) -
+ static_cast<fpt_type>(site1.y());
+ vec_y = static_cast<fpt_type>(site1.x()) -
+ static_cast<fpt_type>(site2.x());
+ sum_x = static_cast<fpt_type>(site1.x()) +
+ static_cast<fpt_type>(site2.x());
+ sum_y = static_cast<fpt_type>(site1.y()) +
+ static_cast<fpt_type>(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();
+ dif[0] = static_cast<fpt_type>(site3.point1().y()) -
+ static_cast<fpt_type>(site1.y());
+ dif[1] = static_cast<fpt_type>(site1.x()) -
+ static_cast<fpt_type>(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();
+ dif[0] = static_cast<fpt_type>(site3.point1().y()) -
+ static_cast<fpt_type>(site2.y());
+ dif[1] = static_cast<fpt_type>(site2.x()) -
+ static_cast<fpt_type>(site3.point1().x());
B = line_a * dif[1] - line_b * dif[0];
sum_AB = A + B;
@@ -692,11 +707,11 @@
c_event.lower_x(0.25 * sqrt_expr_.eval2(cA, cB).get_d() /
(denom.get_d() * std::sqrt(segm_len.get_d())));
}
- return true;
+ return;
}
det = (teta * teta + denom * denom) * A * B * 4;
- double denom_sqr = denom.get_d() * denom.get_d();
+ fpt_type denom_sqr = denom.get_d() * denom.get_d();
if (recompute_c_x || recompute_lower_x) {
cA[0] = sum_x * denom * denom + teta * sum_AB * vec_x;
@@ -729,12 +744,10 @@
c_event.lower_x(0.5 * sqrt_expr_.eval4(cA, cB).get_d() /
(denom_sqr * std::sqrt(segm_len.get_d())));
}
-
- return true;
}
// Recompute parameters of the circle event using high-precision library.
- bool pss(const site_type &site1,
+ void pss(const site_type &site1,
const site_type &site2,
const site_type &site3,
int point_index,
@@ -742,59 +755,72 @@
bool recompute_c_x = true,
bool recompute_c_y = true,
bool recompute_lower_x = true) {
- typedef mpt_wrapper<mpz_class, 8> mpt_type;
- typedef mpt_wrapper<mpf_class, 2> mpf_type;
static mpt_type a[2], b[2], c[2], cA[4], cB[4], orientation, dx, dy, ix, iy;
-
const point_type &segm_start1 = site2.point1(true);
const point_type &segm_end1 = site2.point0(true);
const point_type &segm_start2 = site3.point0(true);
const point_type &segm_end2 = site3.point1(true);
- a[0] = segm_end1.x() - segm_start1.x();
- b[0] = segm_end1.y() - segm_start1.y();
- a[1] = segm_end2.x() - segm_start2.x();
- b[1] = segm_end2.y() - segm_start2.y();
+ a[0] = static_cast<fpt_type>(segm_end1.x()) -
+ static_cast<fpt_type>(segm_start1.x());
+ b[0] = static_cast<fpt_type>(segm_end1.y()) -
+ static_cast<fpt_type>(segm_start1.y());
+ a[1] = static_cast<fpt_type>(segm_end2.x()) -
+ static_cast<fpt_type>(segm_start2.x());
+ b[1] = static_cast<fpt_type>(segm_end2.y()) -
+ static_cast<fpt_type>(segm_start2.y());
orientation = a[1] * b[0] - a[0] * b[1];
if (orientation == 0) {
- double denom = (a[0] * a[0] + b[0] * b[0]).get_d() * 2.0;
- c[0] = b[0] * (segm_start2.x() - segm_start1.x()) -
- a[0] * (segm_start2.y() - segm_start1.y());
- dx = a[0] * (site1.y() - segm_start1.y()) -
- b[0] * (site1.x() - segm_start1.x());
- dy = b[0] * (site1.x() - segm_start2.x()) -
- a[0] * (site1.y() - segm_start2.y());
+ fpt_type denom = (a[0] * a[0] + b[0] * b[0]).get_d() * 2.0;
+ c[0] = b[0] * (static_cast<fpt_type>(segm_start2.x()) -
+ static_cast<fpt_type>(segm_start1.x())) -
+ a[0] * (static_cast<fpt_type>(segm_start2.y()) -
+ static_cast<fpt_type>(segm_start1.y()));
+ dx = a[0] * (static_cast<fpt_type>(site1.y()) -
+ static_cast<fpt_type>(segm_start1.y())) -
+ b[0] * (static_cast<fpt_type>(site1.x()) -
+ static_cast<fpt_type>(segm_start1.x()));
+ dy = b[0] * (static_cast<fpt_type>(site1.x()) -
+ static_cast<fpt_type>(segm_start2.x())) -
+ a[0] * (static_cast<fpt_type>(site1.y()) -
+ static_cast<fpt_type>(segm_start2.y()));
cB[0] = dx * dy;
cB[1] = 1;
if (recompute_c_y) {
cA[0] = b[0] * ((point_index == 2) ? 2 : -2);
- cA[1] = a[0] * a[0] * (segm_start1.y() + segm_start2.y()) -
- a[0] * b[0] * (segm_start1.x() + segm_start2.x() - 2.0 * site1.x()) +
- b[0] * b[0] * (2.0 * site1.y());
- double c_y = sqrt_expr_.eval2(cA, cB).get_d();
+ cA[1] = a[0] * a[0] * (static_cast<fpt_type>(segm_start1.y()) +
+ static_cast<fpt_type>(segm_start2.y())) -
+ a[0] * b[0] * (static_cast<fpt_type>(segm_start1.x()) +
+ static_cast<fpt_type>(segm_start2.x()) -
+ 2.0 * static_cast<fpt_type>(site1.x())) +
+ b[0] * b[0] * (2.0 * static_cast<fpt_type>(site1.y()));
+ fpt_type c_y = sqrt_expr_.eval2(cA, cB).get_d();
c_event.y(c_y / denom);
}
if (recompute_c_x || recompute_lower_x) {
cA[0] = a[0] * ((point_index == 2) ? 2 : -2);
- cA[1] = b[0] * b[0] * (segm_start1.x() + segm_start2.x()) -
- a[0] * b[0] * (segm_start1.y() + segm_start2.y() - 2.0 * site1.y()) +
- a[0] * a[0] * (2.0 * site1.x());
+ cA[1] = b[0] * b[0] * (static_cast<fpt_type>(segm_start1.x()) +
+ static_cast<fpt_type>(segm_start2.x())) -
+ a[0] * b[0] * (static_cast<fpt_type>(segm_start1.y()) +
+ static_cast<fpt_type>(segm_start2.y()) -
+ 2.0 * static_cast<fpt_type>(site1.y())) +
+ a[0] * a[0] * (2.0 * static_cast<fpt_type>(site1.x()));
if (recompute_c_x) {
- double c_x = sqrt_expr_.eval2(cA, cB).get_d();
+ fpt_type c_x = sqrt_expr_.eval2(cA, cB).get_d();
c_event.x(c_x / denom);
}
if (recompute_lower_x) {
cA[2] = (c[0] < 0) ? -c[0] : c[0];
cB[2] = a[0] * a[0] + b[0] * b[0];
- double lower_x = sqrt_expr_.eval3(cA, cB).get_d();
+ fpt_type lower_x = sqrt_expr_.eval3(cA, cB).get_d();
c_event.lower_x(lower_x / denom);
}
}
- return true;
+ return;
}
c[0] = b[0] * segm_end1.x() - a[0] * segm_end1.y();
c[1] = a[1] * segm_end2.y() - b[1] * segm_end2.x();
@@ -803,14 +829,14 @@
dx = ix - orientation * site1.x();
dy = iy - orientation * site1.y();
if ((dx == 0) && (dy == 0)) {
- double denom = orientation.get_d();
- double c_x = ix.get_d() / denom;
- double c_y = iy.get_d() / denom;
+ fpt_type denom = orientation.get_d();
+ fpt_type c_x = ix.get_d() / denom;
+ fpt_type c_y = iy.get_d() / denom;
c_event = circle_type(c_x, c_y, c_x);
- return true;
+ return;
}
- double sign = ((point_index == 2) ? 1 : -1) * ((orientation < 0) ? 1 : -1);
+ fpt_type sign = ((point_index == 2) ? 1 : -1) * ((orientation < 0) ? 1 : -1);
cA[0] = a[1] * -dx + b[1] * -dy;
cA[1] = a[0] * -dx + b[0] * -dy;
cA[2] = sign;
@@ -819,14 +845,14 @@
cB[1] = a[1] * a[1] + b[1] * b[1];
cB[2] = a[0] * a[1] + b[0] * b[1];
cB[3] = (a[0] * dy - b[0] * dx) * (a[1] * dy - b[1] * dx) * -2;
- double temp = sqrt_expr_evaluator_pss<mpt_type, mpf_type>(cA, cB).get_d();
- double denom = temp * orientation.get_d();
+ fpt_type temp = sqrt_expr_evaluator_pss<mpt_type, mpf_type>(cA, cB).get_d();
+ fpt_type denom = temp * orientation.get_d();
if (recompute_c_y) {
cA[0] = b[1] * (dx * dx + dy * dy) - iy * (dx * a[1] + dy * b[1]);
cA[1] = b[0] * (dx * dx + dy * dy) - iy * (dx * a[0] + dy * b[0]);
cA[2] = iy * sign;
- double cy = sqrt_expr_evaluator_pss<mpt_type, mpf_type>(cA, cB).get_d();
+ fpt_type cy = sqrt_expr_evaluator_pss<mpt_type, mpf_type>(cA, cB).get_d();
c_event.y(cy / denom);
}
@@ -836,39 +862,40 @@
cA[2] = ix * sign;
if (recompute_c_x) {
- double cx = sqrt_expr_evaluator_pss<mpt_type, mpf_type>(cA, cB).get_d();
+ fpt_type cx = sqrt_expr_evaluator_pss<mpt_type, mpf_type>(cA, cB).get_d();
c_event.x(cx / denom);
}
if (recompute_lower_x) {
cA[3] = orientation * (dx * dx + dy * dy) * (temp < 0 ? -1 : 1);
- double lower_x = sqrt_expr_evaluator_pss<mpt_type, mpf_type>(cA, cB).get_d();
+ fpt_type lower_x = sqrt_expr_evaluator_pss<mpt_type, mpf_type>(cA, cB).get_d();
c_event.lower_x(lower_x / denom);
}
}
-
- return true;
}
// Recompute parameters of the circle event using high-precision library.
- bool sss(const site_type &site1,
+ void sss(const site_type &site1,
const site_type &site2,
const site_type &site3,
circle_type &c_event,
bool recompute_c_x = true,
bool recompute_c_y = true,
bool recompute_lower_x = true) {
- typedef mpt_wrapper<mpz_class, 8> mpt_type;
- typedef mpt_wrapper<mpf_class, 2> mpf_type;
static mpt_type 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);
+ a[0] = static_cast<fpt_type>(site1.x1(true)) -
+ static_cast<fpt_type>(site1.x0(true));
+ a[1] = static_cast<fpt_type>(site2.x1(true)) -
+ static_cast<fpt_type>(site2.x0(true));
+ a[2] = static_cast<fpt_type>(site3.x1(true)) -
+ static_cast<fpt_type>(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);
+ b[0] = static_cast<fpt_type>(site1.y1(true)) -
+ static_cast<fpt_type>(site1.y0(true));
+ b[1] = static_cast<fpt_type>(site2.y1(true)) -
+ static_cast<fpt_type>(site2.y0(true));
+ b[2] = static_cast<fpt_type>(site3.y1(true)) -
+ static_cast<fpt_type>(site3.y0(true));
c[0] = mpt_cross(site1.x0(true), site1.y0(true), site1.x1(true), site1.y1(true));
c[1] = mpt_cross(site2.x0(true), site2.y0(true), site2.x1(true), site2.y1(true));
@@ -883,7 +910,7 @@
int k = (i+2) % 3;
cross[i] = a[j] * b[k] - a[k] * b[j];
}
- double denom = sqrt_expr_.eval3(cross, sqr_len).get_d();
+ fpt_type denom = sqrt_expr_.eval3(cross, sqr_len).get_d();
if (recompute_c_y) {
for (int i = 0; i < 3; ++i) {
@@ -891,7 +918,7 @@
int k = (i+2) % 3;
cross[i] = b[j] * c[k] - b[k] * c[j];
}
- double c_y = sqrt_expr_.eval3(cross, sqr_len).get_d();
+ fpt_type c_y = sqrt_expr_.eval3(cross, sqr_len).get_d();
c_event.y(c_y / denom);
}
@@ -907,24 +934,22 @@
}
if (recompute_c_x) {
- double c_x = sqrt_expr_.eval3(cross, sqr_len).get_d();
+ fpt_type c_x = sqrt_expr_.eval3(cross, sqr_len).get_d();
c_event.x(c_x / denom);
}
if (recompute_lower_x) {
sqr_len[3] = 1;
- double lower_x = sqrt_expr_.eval4(cross, sqr_len).get_d();
+ fpt_type lower_x = sqrt_expr_.eval4(cross, sqr_len).get_d();
c_event.lower_x(lower_x / denom);
}
}
-
- return true;
}
private:
template <typename T>
mpt_wrapper<mpz_class, 8> &mpt_cross(T a0, T b0, T a1, T b1) {
- static mpt_wrapper<mpz_class, 8> temp[2];
+ static mpt_type temp[2];
temp[0] = a0;
temp[1] = b0;
temp[0] = temp[0] * b1;
@@ -994,6 +1019,8 @@
template <typename Site, typename Circle>
class lazy_circle_formation_functor {
public:
+ typedef robust_fpt<fpt_type> robust_fpt_type;
+ typedef robust_dif<robust_fpt_type> robust_dif_type;
typedef typename Site::point_type point_type;
typedef Site site_type;
typedef Circle circle_type;
@@ -1002,108 +1029,128 @@
// Find parameters of the inscribed circle that is tangent to three
// point sites.
- bool ppp(const site_type &site1,
+ void ppp(const site_type &site1,
const site_type &site2,
const site_type &site3,
circle_type &c_event) {
- double dif_x1 = site1.x() - site2.x();
- double dif_x2 = site2.x() - site3.x();
- double dif_y1 = site1.y() - site2.y();
- double dif_y2 = site2.y() - site3.y();
- double orientation = robust_cross_product(dif_x1, dif_y1, dif_x2, dif_y2);
- robust_fpt<double> inv_orientation(0.5 / orientation, 3.0);
- double sum_x1 = site1.x() + site2.x();
- double sum_x2 = site2.x() + site3.x();
- double sum_y1 = site1.y() + site2.y();
- double sum_y2 = site2.y() + site3.y();
- double dif_x3 = site1.x() - site3.x();
- double dif_y3 = site1.y() - site3.y();
- robust_dif< robust_fpt<double> > c_x, c_y;
- c_x += robust_fpt<double>(dif_x1 * sum_x1 * dif_y2, 2.0);
- c_x += robust_fpt<double>(dif_y1 * sum_y1 * dif_y2, 2.0);
- c_x -= robust_fpt<double>(dif_x2 * sum_x2 * dif_y1, 2.0);
- c_x -= robust_fpt<double>(dif_y2 * sum_y2 * dif_y1, 2.0);
- c_y += robust_fpt<double>(dif_x2 * sum_x2 * dif_x1, 2.0);
- c_y += robust_fpt<double>(dif_y2 * sum_y2 * dif_x1, 2.0);
- c_y -= robust_fpt<double>(dif_x1 * sum_x1 * dif_x2, 2.0);
- c_y -= robust_fpt<double>(dif_y1 * sum_y1 * dif_x2, 2.0);
- robust_dif< robust_fpt<double> > lower_x(c_x);
- lower_x -= robust_fpt<double>(std::sqrt(sqr_distance(dif_x1, dif_y1) *
- sqr_distance(dif_x2, dif_y2) *
- sqr_distance(dif_x3, dif_y3)), 5.0);
+ fpt_type dif_x1 = static_cast<fpt_type>(site1.x()) -
+ static_cast<fpt_type>(site2.x());
+ fpt_type dif_x2 = static_cast<fpt_type>(site2.x()) -
+ static_cast<fpt_type>(site3.x());
+ fpt_type dif_y1 = static_cast<fpt_type>(site1.y()) -
+ static_cast<fpt_type>(site2.y());
+ 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);
+ 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()) +
+ static_cast<fpt_type>(site3.x());
+ fpt_type sum_y1 = static_cast<fpt_type>(site1.y()) +
+ static_cast<fpt_type>(site2.y());
+ fpt_type sum_y2 = static_cast<fpt_type>(site2.y()) +
+ static_cast<fpt_type>(site3.y());
+ fpt_type dif_x3 = static_cast<fpt_type>(site1.x()) -
+ static_cast<fpt_type>(site3.x());
+ fpt_type dif_y3 = static_cast<fpt_type>(site1.y()) -
+ static_cast<fpt_type>(site3.y());
+ robust_dif_type c_x, c_y;
+ c_x += robust_fpt_type(dif_x1 * sum_x1 * dif_y2, 2.0);
+ c_x += robust_fpt_type(dif_y1 * sum_y1 * dif_y2, 2.0);
+ c_x -= robust_fpt_type(dif_x2 * sum_x2 * dif_y1, 2.0);
+ c_x -= robust_fpt_type(dif_y2 * sum_y2 * dif_y1, 2.0);
+ c_y += robust_fpt_type(dif_x2 * sum_x2 * dif_x1, 2.0);
+ c_y += robust_fpt_type(dif_y2 * sum_y2 * dif_x1, 2.0);
+ c_y -= robust_fpt_type(dif_x1 * sum_x1 * dif_x2, 2.0);
+ c_y -= robust_fpt_type(dif_y1 * sum_y1 * dif_x2, 2.0);
+ robust_dif_type lower_x(c_x);
+ lower_x -= robust_fpt_type(std::sqrt(sqr_distance(dif_x1, dif_y1) *
+ sqr_distance(dif_x2, dif_y2) *
+ sqr_distance(dif_x3, dif_y3)), 5.0);
c_event = circle_type(c_x.dif().fpv() * inv_orientation.fpv(),
- c_y.dif().fpv() * inv_orientation.fpv(),
- lower_x.dif().fpv() * inv_orientation.fpv());
- bool recompute_c_x = c_x.dif().ulp() >= 128;
- bool recompute_c_y = c_y.dif().ulp() >= 128;
- bool recompute_lower_x = lower_x.dif().ulp() >= 128;
+ c_y.dif().fpv() * inv_orientation.fpv(),
+ lower_x.dif().fpv() * inv_orientation.fpv());
+ bool recompute_c_x = c_x.dif().ulp() > fULPS;
+ bool recompute_c_y = c_y.dif().ulp() > fULPS;
+ bool recompute_lower_x = lower_x.dif().ulp() > fULPS;
if (recompute_c_x || recompute_c_y || recompute_lower_x) {
- return exact_circle_formation_functor_.ppp(
+ exact_circle_formation_functor_.ppp(
site1, site2, site3, c_event, recompute_c_x, recompute_c_y, recompute_lower_x);
}
- return true;
}
// Find parameters of the inscribed circle that is tangent to two
// point sites and on segment site.
- bool pps(const site_type &site1,
+ void pps(const site_type &site1,
const site_type &site2,
const site_type &site3,
int segment_index,
circle_type &c_event) {
- 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();
- robust_fpt<double> teta(robust_cross_product(line_a, line_b, -vec_y, vec_x), 2.0);
- robust_fpt<double> A(robust_cross_product(line_a, line_b,
- site3.point1().y() - site1.y(),
- site1.x() - site3.point1().x()), 2.0);
- robust_fpt<double> B(robust_cross_product(line_a, line_b,
- site3.point1().y() - site2.y(),
- site2.x() - site3.point1().x()), 2.0);
- robust_fpt<double> denom(robust_cross_product(vec_x, vec_y, line_a, line_b), 2.0);
- robust_fpt<double> inv_segm_len(1.0 / std::sqrt(sqr_distance(line_a, line_b)), 3.0);
- robust_dif< robust_fpt<double> > t;
+ fpt_type line_a = static_cast<fpt_type>(site3.point1(true).y()) -
+ static_cast<fpt_type>(site3.point0(true).y());
+ fpt_type line_b = static_cast<fpt_type>(site3.point0(true).x()) -
+ static_cast<fpt_type>(site3.point1(true).x());
+ fpt_type vec_x = static_cast<fpt_type>(site2.y()) -
+ 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 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);
+ 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);
+ 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) {
- t += teta / (robust_fpt<double>(8.0) * A);
- t -= A / (robust_fpt<double>(2.0) * teta);
+ t += teta / (robust_fpt_type(8.0, false) * A);
+ t -= A / (robust_fpt_type(2.0, false) * teta);
} else {
- robust_fpt<double> det = ((teta * teta + denom * denom) * A * B).sqrt();
+ robust_fpt_type det = ((teta * teta + denom * denom) * A * B).sqrt();
if (segment_index == 2) {
t -= det / (denom * denom);
} else {
t += det / (denom * denom);
}
- t += teta * (A + B) / (robust_fpt<double>(2.0) * denom * denom);
+ t += teta * (A + B) / (robust_fpt_type(2.0, false) * denom * denom);
}
- robust_dif< robust_fpt<double> > c_x, c_y;
- c_x += robust_fpt<double>(0.5 * (site1.x() + site2.x()));
- c_x += robust_fpt<double>(vec_x) * t;
- c_y += robust_fpt<double>(0.5 * (site1.y() + site2.y()));
- c_y += robust_fpt<double>(vec_y) * t;
- robust_dif< robust_fpt<double> > r, lower_x(c_x);
- r -= robust_fpt<double>(line_a) * robust_fpt<double>(site3.x0());
- r -= robust_fpt<double>(line_b) * robust_fpt<double>(site3.y0());
- r += robust_fpt<double>(line_a) * c_x;
- r += robust_fpt<double>(line_b) * c_y;
+ robust_dif_type c_x, c_y;
+ c_x += robust_fpt_type(0.5 * (
+ static_cast<fpt_type>(site1.x()) +
+ static_cast<fpt_type>(site2.x())), false);
+ c_x += robust_fpt_type(vec_x, false) * t;
+ c_y += robust_fpt_type(0.5 * (
+ static_cast<fpt_type>(site1.y()) +
+ static_cast<fpt_type>(site2.y())), false);
+ c_y += robust_fpt_type(vec_y, false) * t;
+ robust_dif_type r, lower_x(c_x);
+ r -= robust_fpt_type(line_a, false) * robust_fpt_type(site3.x0(), false);
+ r -= robust_fpt_type(line_b, false) * robust_fpt_type(site3.y0(), false);
+ r += robust_fpt_type(line_a, false) * c_x;
+ r += robust_fpt_type(line_b, false) * c_y;
r.abs();
lower_x += r * inv_segm_len;
c_event = circle_type(c_x.dif().fpv(), c_y.dif().fpv(), lower_x.dif().fpv());
- bool recompute_c_x = c_x.dif().ulp() >= 128;
- bool recompute_c_y = c_y.dif().ulp() >= 128;
- bool recompute_lower_x = lower_x.dif().ulp() >= 128;
+ bool recompute_c_x = c_x.dif().ulp() > fULPS;
+ bool recompute_c_y = c_y.dif().ulp() > fULPS;
+ bool recompute_lower_x = lower_x.dif().ulp() > fULPS;
if (recompute_c_x || recompute_c_y || recompute_lower_x) {
- return exact_circle_formation_functor_.pps(
+ exact_circle_formation_functor_.pps(
site1, site2, site3, segment_index, c_event,
recompute_c_x, recompute_c_y, recompute_lower_x);
}
- return true;
}
// Find parameters of the inscribed circle that is tangent to one
// point site and two segment sites.
- bool pss(const site_type &site1,
+ void pss(const site_type &site1,
const site_type &site2,
const site_type &site3,
int point_index,
@@ -1112,69 +1159,106 @@
const point_type &segm_end1 = site2.point0(true);
const point_type &segm_start2 = site3.point0(true);
const point_type &segm_end2 = site3.point1(true);
- double a1 = segm_end1.x() - segm_start1.x();
- double b1 = segm_end1.y() - segm_start1.y();
- double a2 = segm_end2.x() - segm_start2.x();
- double b2 = segm_end2.y() - segm_start2.y();
+ fpt_type a1 = static_cast<fpt_type>(segm_end1.x()) -
+ static_cast<fpt_type>(segm_start1.x());
+ fpt_type b1 = static_cast<fpt_type>(segm_end1.y()) -
+ static_cast<fpt_type>(segm_start1.y());
+ fpt_type a2 = static_cast<fpt_type>(segm_end2.x()) -
+ static_cast<fpt_type>(segm_start2.x());
+ 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<double> orientation(robust_cross_product(b1, a1, b2, a2), 2.0);
+ robust_fpt_type orientation(robust_cross_product(b1, a1, b2, a2), 2.0);
if (get_orientation(orientation) == COLLINEAR) {
- robust_fpt<double> a(a1 * a1 + b1 * b1, 2.0);
- robust_fpt<double> c(robust_cross_product(b1, a1, segm_start2.y() - segm_start1.y(),
- segm_start2.x() - segm_start1.x()), 2.0);
- robust_fpt<double> det(robust_cross_product(a1, b1, site1.x() - segm_start1.x(),
- site1.y() - segm_start1.y()) *
- robust_cross_product(b1, a1, site1.y() - segm_start2.y(),
- site1.x() - segm_start2.x()), 5.0);
- robust_dif< robust_fpt<double> > t;
- t -= robust_fpt<double>(a1) * robust_fpt<double>((segm_start1.x() + segm_start2.x()) * 0.5 - site1.x());
- t -= robust_fpt<double>(b1) * robust_fpt<double>((segm_start1.y() + segm_start2.y()) * 0.5 - site1.y());
+ 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);
+ robust_fpt_type det(robust_cross_product(a1, b1,
+ static_cast<fpt_type>(site1.x()) -
+ static_cast<fpt_type>(segm_start1.x()),
+ static_cast<fpt_type>(site1.y()) -
+ static_cast<fpt_type>(segm_start1.y())) *
+ robust_cross_product(b1, a1,
+ 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);
+ robust_dif_type t;
+ t -= robust_fpt_type(a1, false) * robust_fpt_type(
+ (static_cast<fpt_type>(segm_start1.x()) +
+ static_cast<fpt_type>(segm_start2.x())) * 0.5 -
+ static_cast<fpt_type>(site1.x()), false);
+ t -= robust_fpt_type(b1, false) * robust_fpt_type((
+ static_cast<fpt_type>(segm_start1.y()) +
+ static_cast<fpt_type>(segm_start2.y())) * 0.5 -
+ static_cast<fpt_type>(site1.y()), false);
if (point_index == 2) {
t += det.sqrt();
} else {
t -= det.sqrt();
}
t /= a;
- robust_dif< robust_fpt<double> > c_x, c_y;
- c_x += robust_fpt<double>(0.5 * (segm_start1.x() + segm_start2.x()));
- c_x += robust_fpt<double>(a1) * t;
- c_y += robust_fpt<double>(0.5 * (segm_start1.y() + segm_start2.y()));
- c_y += robust_fpt<double>(b1) * t;
- robust_dif< robust_fpt<double> > lower_x(c_x);
- lower_x += robust_fpt<double>(0.5) * c.fabs() / a.sqrt();
- recompute_c_x = c_x.dif().ulp() > 128;
- recompute_c_y = c_y.dif().ulp() > 128;
- recompute_lower_x = lower_x.dif().ulp() > 128;
+ robust_dif_type c_x, c_y;
+ c_x += robust_fpt_type(0.5 * (
+ static_cast<fpt_type>(segm_start1.x()) +
+ static_cast<fpt_type>(segm_start2.x())), false);
+ c_x += robust_fpt_type(a1, false) * t;
+ c_y += robust_fpt_type(0.5 * (
+ static_cast<fpt_type>(segm_start1.y()) +
+ static_cast<fpt_type>(segm_start2.y())), false);
+ c_y += robust_fpt_type(b1, false) * t;
+ robust_dif_type lower_x(c_x);
+ lower_x += robust_fpt_type(0.5, false) * c.fabs() / a.sqrt();
+ recompute_c_x = c_x.dif().ulp() > fULPS;
+ recompute_c_y = c_y.dif().ulp() > fULPS;
+ recompute_lower_x = lower_x.dif().ulp() > fULPS;
c_event = circle_type(c_x.dif().fpv(), c_y.dif().fpv(), lower_x.dif().fpv());
} else {
- robust_fpt<double> sqr_sum1(std::sqrt(a1 * a1 + b1 * b1), 2.0);
- robust_fpt<double> sqr_sum2(std::sqrt(a2 * a2 + b2 * b2), 2.0);
- robust_fpt<double> a(robust_cross_product(a1, b1, -b2, a2), 2.0);
+ 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);
if (a >= 0) {
a += sqr_sum1 * sqr_sum2;
} else {
a = (orientation * orientation) / (sqr_sum1 * sqr_sum2 - a);
}
- robust_fpt<double> or1(robust_cross_product(b1, a1, segm_end1.y() - site1.y(),
- segm_end1.x() - site1.x()), 2.0);
- robust_fpt<double> or2(robust_cross_product(a2, b2, segm_end2.x() - site1.x(),
- segm_end2.y() - site1.y()), 2.0);
- robust_fpt<double> det = robust_fpt<double>(2.0) * a * or1 * or2;
- robust_fpt<double> c1(robust_cross_product(b1, a1, segm_end1.y(), segm_end1.x()), 2.0);
- robust_fpt<double> c2(robust_cross_product(a2, b2, segm_end2.x(), segm_end2.y()), 2.0);
- robust_fpt<double> inv_orientation = robust_fpt<double>(1.0) / orientation;
- robust_dif< robust_fpt<double> > t, b, ix, iy;
- ix += robust_fpt<double>(a2) * c1 * inv_orientation;
- ix += robust_fpt<double>(a1) * c2 * inv_orientation;
- iy += robust_fpt<double>(b1) * c2 * inv_orientation;
- iy += robust_fpt<double>(b2) * c1 * inv_orientation;
-
- b += ix * (robust_fpt<double>(a1) * sqr_sum2);
- b += ix * (robust_fpt<double>(a2) * sqr_sum1);
- b += iy * (robust_fpt<double>(b1) * sqr_sum2);
- b += iy * (robust_fpt<double>(b2) * sqr_sum1);
- b -= sqr_sum1 * robust_fpt<double>(robust_cross_product(a2, b2, -site1.y(), site1.x()), 2.0);
- b -= sqr_sum2 * robust_fpt<double>(robust_cross_product(a1, b1, -site1.y(), site1.x()), 2.0);
+ robust_fpt_type or1(robust_cross_product(b1, a1,
+ 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);
+ 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);
+ 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);
+ 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);
+ 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;
+ ix += robust_fpt_type(a1, false) * c2 * inv_orientation;
+ iy += robust_fpt_type(b1, false) * c2 * inv_orientation;
+ iy += robust_fpt_type(b2, false) * c1 * inv_orientation;
+
+ b += ix * (robust_fpt_type(a1, false) * sqr_sum2);
+ b += ix * (robust_fpt_type(a2, false) * sqr_sum1);
+ 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,
+ static_cast<fpt_type>(-site1.y()),
+ static_cast<fpt_type>(site1.x())), 2.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);
t -= b;
if (point_index == 2) {
t += det.sqrt();
@@ -1182,52 +1266,57 @@
t -= det.sqrt();
}
t /= (a * a);
- robust_dif< robust_fpt<double> > c_x(ix), c_y(iy);
- c_x += t * (robust_fpt<double>(a1) * sqr_sum2);
- c_x += t * (robust_fpt<double>(a2) * sqr_sum1);
- c_y += t * (robust_fpt<double>(b1) * sqr_sum2);
- c_y += t * (robust_fpt<double>(b2) * sqr_sum1);
+ robust_dif_type c_x(ix), c_y(iy);
+ c_x += t * (robust_fpt_type(a1, false) * sqr_sum2);
+ c_x += t * (robust_fpt_type(a2, false) * sqr_sum1);
+ c_y += t * (robust_fpt_type(b1, false) * sqr_sum2);
+ c_y += t * (robust_fpt_type(b2, false) * sqr_sum1);
t.abs();
- robust_dif< robust_fpt<double> > lower_x(c_x);
+ robust_dif_type lower_x(c_x);
lower_x += t * orientation.fabs();
- recompute_c_x = c_x.dif().ulp() > 128;
- recompute_c_y = c_y.dif().ulp() > 128;
- recompute_lower_x = lower_x.dif().ulp() > 128;
+ recompute_c_x = c_x.dif().ulp() > fULPS;
+ recompute_c_y = c_y.dif().ulp() > fULPS;
+ recompute_lower_x = lower_x.dif().ulp() > fULPS;
c_event = circle_type(c_x.dif().fpv(), c_y.dif().fpv(), lower_x.dif().fpv());
}
if (recompute_c_x || recompute_c_y || recompute_lower_x) {
- return exact_circle_formation_functor_.pss(
+ exact_circle_formation_functor_.pss(
site1, site2, site3, point_index, c_event,
recompute_c_x, recompute_c_y, recompute_lower_x);
}
- return true;
}
// Find parameters of the inscribed circle that is tangent to three
// segment sites.
- bool sss(const site_type &site1,
+ void sss(const site_type &site1,
const site_type &site2,
const site_type &site3,
circle_type &c_event) {
- robust_fpt<double> a1(site1.x1(true) - site1.x0(true), 0.0);
- robust_fpt<double> b1(site1.y1(true) - site1.y0(true), 0.0);
- robust_fpt<double> c1(robust_cross_product(site1.x0(true), site1.y0(true), site1.x1(true), site1.y1(true)), 2.0);
+ robust_fpt_type a1(static_cast<fpt_type>(site1.x1(true)) -
+ 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<double> a2(site2.x1(true) - site2.x0(true), 0.0);
- robust_fpt<double> b2(site2.y1(true) - site2.y0(true), 0.0);
- robust_fpt<double> c2(robust_cross_product(site2.x0(true), site2.y0(true), site2.x1(true), site2.y1(true)), 2.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<double> a3(site3.x1(true) - site3.x0(true), 0.0);
- robust_fpt<double> b3(site3.y1(true) - site3.y0(true), 0.0);
- robust_fpt<double> c3(robust_cross_product(site3.x0(true), site3.y0(true), site3.x1(true), site3.y1(true)), 2.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<double> len1 = (a1 * a1 + b1 * b1).sqrt();
- robust_fpt<double> len2 = (a2 * a2 + b2 * b2).sqrt();
- robust_fpt<double> len3 = (a3 * a3 + b3 * b3).sqrt();
- robust_fpt<double> cross_12(robust_cross_product(a1.fpv(), b1.fpv(), a2.fpv(), b2.fpv()), 2.0);
- robust_fpt<double> cross_23(robust_cross_product(a2.fpv(), b2.fpv(), a3.fpv(), b3.fpv()), 2.0);
- robust_fpt<double> cross_31(robust_cross_product(a3.fpv(), b3.fpv(), a1.fpv(), b1.fpv()), 2.0);
- robust_dif< robust_fpt<double> > denom, c_x, c_y, r;
+ 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_dif_type denom, c_x, c_y, r;
// denom = cross_12 * len3 + cross_23 * len1 + cross_31 * len2.
denom += cross_12 * len3;
@@ -1251,20 +1340,19 @@
c_y -= b3 * c2 * len1;
c_y += b3 * c1 * len2;
c_y -= b1 * c3 * len2;
- robust_dif< robust_fpt<double> > lower_x(c_x + r);
- bool recompute_c_x = c_x.dif().ulp() > 128;
- bool recompute_c_y = c_y.dif().ulp() > 128;
- bool recompute_lower_x = lower_x.dif().ulp() > 128;
- bool recompute_denom = denom.dif().ulp() > 128;
+ robust_dif_type lower_x(c_x + r);
+ bool recompute_c_x = c_x.dif().ulp() > fULPS;
+ bool recompute_c_y = c_y.dif().ulp() > fULPS;
+ bool recompute_lower_x = lower_x.dif().ulp() > fULPS;
+ bool recompute_denom = denom.dif().ulp() > fULPS;
c_event = circle_type(c_x.dif().fpv() / denom.dif().fpv(),
- c_y.dif().fpv() / denom.dif().fpv(),
- lower_x.dif().fpv() / denom.dif().fpv());
+ c_y.dif().fpv() / denom.dif().fpv(),
+ lower_x.dif().fpv() / denom.dif().fpv());
if (recompute_c_x || recompute_c_y || recompute_lower_x || recompute_denom) {
- return exact_circle_formation_functor_.sss(
+ exact_circle_formation_functor_.sss(
site1, site2, site3, c_event,
recompute_c_x, recompute_c_y, recompute_lower_x);
}
- return true;
}
private:
@@ -1369,6 +1457,11 @@
}
};
+const voronoi_calc_kernel<int>::fpt_type fpt_type voronoi_calc_kernel<int>::fULPS =
+ voronoi_calc_kernel<int>::ULPS;
+const voronoi_calc_kernel<int>::fpt_type fpt_type voronoi_calc_kernel<int>::fULPSx2 =
+ voronoi_calc_kernel<int>::ULPSx2;
+
} // detail
} // polygon
} // boost
Modified: sandbox/gtl/libs/polygon/test/Jamfile.v2
==============================================================================
--- sandbox/gtl/libs/polygon/test/Jamfile.v2 (original)
+++ sandbox/gtl/libs/polygon/test/Jamfile.v2 2011-10-17 18:17:05 EDT (Mon, 17 Oct 2011)
@@ -21,7 +21,7 @@
[ run gtl_boost_unit_test.cpp ]
;
-alias "benchmark"
+alias "voronoi-benchmark"
:
[ run voronoi_benchmark_test.cpp ]
;
Modified: sandbox/gtl/libs/polygon/test/voronoi_calc_kernel_test.cpp
==============================================================================
--- sandbox/gtl/libs/polygon/test/voronoi_calc_kernel_test.cpp (original)
+++ sandbox/gtl/libs/polygon/test/voronoi_calc_kernel_test.cpp 2011-10-17 18:17:05 EDT (Mon, 17 Oct 2011)
@@ -31,7 +31,11 @@
distance_predicate_type distance_predicate;
node_comparison_type node_comparison;
-//typedef VCK::circle_formation_predicate lazy_predicate
+typedef VCK::circle_existence_predicate<site_type> CEP_type;
+typedef VCK::gmp_circle_formation_functor<site_type, circle_type> GMP_CFF_type;
+typedef VCK::lazy_circle_formation_functor<site_type, circle_type> lazy_CFF_type;
+VCK::circle_formation_predicate<site_type, circle_type, CEP_type, GMP_CFF_type> gmp_predicate;
+VCK::circle_formation_predicate<site_type, circle_type, CEP_type, lazy_CFF_type> lazy_predicate;
#define CHECK_ORIENTATION(P1, P2, P3, R1, R2) \
BOOST_CHECK_EQUAL(VCK::get_orientation(P1, P2, P3) == R1, true); \
@@ -55,6 +59,22 @@
BOOST_CHECK_EQUAL(node_comparison(nodes[i], node), !res[i]); \
}
+#define CHECK_CIRCLE(circle, c_x, c_y, l_x) \
+ BOOST_CHECK_EQUAL(almost_equal(c1.x(), c_x, 10), true); \
+ BOOST_CHECK_EQUAL(almost_equal(c1.y(), c_y, 10), true); \
+ BOOST_CHECK_EQUAL(almost_equal(c1.lower_x(), l_x, 10), true)
+
+#define CHECK_CIRCLE_EXISTENCE(s1, s2, s3, RES) \
+ { circle_type c1; \
+ BOOST_CHECK_EQUAL(lazy_predicate(s1, s2, s3, c1), RES); }
+
+#define CHECK_CIRCLE_FORMATION_PREDICATE(s1, s2, s3, c_x, c_y, l_x) \
+ { circle_type c1, c2; \
+ BOOST_CHECK_EQUAL(gmp_predicate(s1, s2, s3, c1), true); \
+ BOOST_CHECK_EQUAL(lazy_predicate(s1, s2, s3, c2), true); \
+ CHECK_CIRCLE(c1, c_x, c_y, l_x); \
+ CHECK_CIRCLE(c2, c_x, c_y, l_x); }
+
BOOST_AUTO_TEST_CASE(orientation_test) {
int min_int = std::numeric_limits<int>::min();
int max_int = std::numeric_limits<int>::max();
@@ -343,3 +363,90 @@
bool res[] = {false, true, true, true};
CHECK_NODE_COMPARISON(node, nodes, res, 4);
}
+
+BOOST_AUTO_TEST_CASE(circle_formation_predicate_test1) {
+ site_type site1(0, 0, 0);
+ site_type site2(-8, 0, 1);
+ site_type site3(0, 6, 2);
+ CHECK_CIRCLE_FORMATION_PREDICATE(site1, site2, site3, -4.0, 3.0, 1.0);
+}
+
+BOOST_AUTO_TEST_CASE(circle_formation_predicate_test2) {
+ int min_int = std::numeric_limits<int>::min();
+ int max_int = std::numeric_limits<int>::max();
+ site_type site1(min_int, min_int, 0);
+ site_type site2(min_int, max_int, 1);
+ site_type site3(max_int-1, max_int-1, 2);
+ site_type site4(max_int, max_int, 3);
+ CHECK_CIRCLE_EXISTENCE(site1, site2, site4, true);
+ CHECK_CIRCLE_EXISTENCE(site1, site3, site4, false);
+}
+
+BOOST_AUTO_TEST_CASE(circle_formation_predicate_test3) {
+ site_type site1(-4, 0, 0);
+ site_type site2(0, 4, 1);
+ site_type site3(site1.point0(), site2.point0(), 2);
+ CHECK_CIRCLE_EXISTENCE(site1, site3, site2, false);
+ site_type site4(-2, 0, 3);
+ site_type site5(0, 2, 4);
+ CHECK_CIRCLE_EXISTENCE(site3, site4, site5, false);
+ CHECK_CIRCLE_EXISTENCE(site4, site5, site3, false);
+}
+
+BOOST_AUTO_TEST_CASE(circle_formation_predicate_test5) {
+ site_type site1(-4, 0, -4, 20, 0);
+ site_type site2(-2, 10, 1);
+ site_type site3(4, 10, 0);
+ CHECK_CIRCLE_FORMATION_PREDICATE(site1, site2, site3, 1.0, 6.0, 6.0);
+ CHECK_CIRCLE_FORMATION_PREDICATE(site3, site2, site1, 1.0, 14.0, 6.0);
+}
+
+BOOST_AUTO_TEST_CASE(circle_formation_predicate_test6) {
+ site_type site1(1, 0, 7, 0, 0);
+ site1.inverse();
+ site_type site2(-2, 4, 10, 4, 1);
+ site_type site3(6, 2, 2);
+ site_type site4(1, 0, 2);
+ CHECK_CIRCLE_FORMATION_PREDICATE(site3, site1, site2, 4.0, 2.0, 6.0);
+ CHECK_CIRCLE_FORMATION_PREDICATE(site4, site2, site1, 1.0, 2.0, 3.0);
+}
+
+BOOST_AUTO_TEST_CASE(circle_formation_predicate_test7) {
+ site_type site1(-1, 2, 8, -10, 0);
+ site1.inverse();
+ site_type site2(-1, 0, 8, 12, 1);
+ site_type site3(1, 1, 2);
+ CHECK_CIRCLE_FORMATION_PREDICATE(site3, site2, site1, 6.0, 1.0, 11.0);
+}
+
+BOOST_AUTO_TEST_CASE(circle_formation_predicate_test8) {
+ site_type site1(1, 0, 6, 0, 0);
+ site1.inverse();
+ site_type site2(-6, 4, 0, 12, 1);
+ site_type site3(1, 0, 2);
+ CHECK_CIRCLE_FORMATION_PREDICATE(site3, site2, site1, 1.0, 5.0, 6.0);
+}
+
+BOOST_AUTO_TEST_CASE(circle_formation_predicate_test9) {
+ site_type site1(1, 0, 5, 0, 0);
+ site1.inverse();
+ site_type site2(8, 6, 0, 12, 1);
+ site_type site3(1, 0, 2);
+ CHECK_CIRCLE_FORMATION_PREDICATE(site3, site2, site1, 1.0, 5.0, 6.0);
+}
+
+BOOST_AUTO_TEST_CASE(circle_formation_predicate_test10) {
+ site_type site1(4, 0, 0, 0, 0);
+ site_type site2(0, 0, 0, 4, 1);
+ site_type site3(0, 4, 4, 4, 2);
+ site1.inverse();
+ CHECK_CIRCLE_FORMATION_PREDICATE(site1, site2, site3, 2.0, 2.0, 4.0);
+}
+
+BOOST_AUTO_TEST_CASE(circle_formation_predicate_test11) {
+ site_type site1(41, 30, 1, 0, 0);
+ site_type site2(-39, 30, 1, 60, 1);
+ site_type site3(1, 60, 41, 30, 2);
+ site1.inverse();
+ CHECK_CIRCLE_FORMATION_PREDICATE(site1, site2, site3, 1.0, 30.0, 25.0);
+}