|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r75406 - in sandbox/gtl: boost/polygon boost/polygon/detail libs/polygon/test
From: sydorchuk.andriy_at_[hidden]
Date: 2011-11-08 09:16:59
Author: asydorchuk
Date: 2011-11-08 09:16:57 EST (Tue, 08 Nov 2011)
New Revision: 75406
URL: http://svn.boost.org/trac/boost/changeset/75406
Log:
Renaming voronoi_calc_kernel->voronoi_calc_utils.
Renaming voronoi_fpt_kernel->voronoi_robust_fpt.
Minor comment updates, code refactoring.
Added:
sandbox/gtl/boost/polygon/detail/voronoi_calc_utils.hpp
- copied, changed from r75182, /sandbox/gtl/boost/polygon/detail/voronoi_calc_kernel.hpp
sandbox/gtl/boost/polygon/detail/voronoi_robust_fpt.hpp
- copied, changed from r75144, /sandbox/gtl/boost/polygon/detail/voronoi_fpt_kernel.hpp
sandbox/gtl/libs/polygon/test/voronoi_calc_utils_test.cpp
- copied, changed from r75182, /sandbox/gtl/libs/polygon/test/voronoi_calc_kernel_test.cpp
sandbox/gtl/libs/polygon/test/voronoi_robust_fpt_test.cpp
- copied, changed from r75083, /sandbox/gtl/libs/polygon/test/voronoi_fpt_kernel_test.cpp
Removed:
sandbox/gtl/boost/polygon/detail/voronoi_calc_kernel.hpp
sandbox/gtl/boost/polygon/detail/voronoi_fpt_kernel.hpp
sandbox/gtl/libs/polygon/test/voronoi_calc_kernel_test.cpp
sandbox/gtl/libs/polygon/test/voronoi_fpt_kernel_test.cpp
Text files modified:
sandbox/gtl/boost/polygon/detail/voronoi_calc_utils.hpp | 22 ++++++++++----------
sandbox/gtl/boost/polygon/detail/voronoi_robust_fpt.hpp | 19 ++++++++++-------
sandbox/gtl/boost/polygon/detail/voronoi_structures.hpp | 12 ++++++++--
sandbox/gtl/boost/polygon/voronoi_builder.hpp | 19 ++++++++---------
sandbox/gtl/boost/polygon/voronoi_diagram.hpp | 2
sandbox/gtl/libs/polygon/test/Jamfile.v2 | 4 +-
sandbox/gtl/libs/polygon/test/voronoi_calc_utils_test.cpp | 42 ++++++++++++++++++++--------------------
sandbox/gtl/libs/polygon/test/voronoi_robust_fpt_test.cpp | 8 +++---
8 files changed, 68 insertions(+), 60 deletions(-)
Deleted: sandbox/gtl/boost/polygon/detail/voronoi_calc_kernel.hpp
==============================================================================
--- sandbox/gtl/boost/polygon/detail/voronoi_calc_kernel.hpp 2011-11-08 09:16:57 EST (Tue, 08 Nov 2011)
+++ (empty file)
@@ -1,1406 +0,0 @@
-// Boost.Polygon library detail/voronoi_calc_kernel.hpp header file
-
-// Copyright Andrii Sydorchuk 2010-2011.
-// 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_POLYGON_VORONOI_CALC_KERNEL
-#define BOOST_POLYGON_VORONOI_CALC_KERNEL
-
-#pragma warning(disable:4800)
-#include <gmpxx.h>
-
-#include "mpt_wrapper.hpp"
-#include "voronoi_fpt_kernel.hpp"
-
-namespace boost {
-namespace polygon {
-namespace detail {
-
-template <typename T>
-class voronoi_calc_kernel;
-
-// Predicates kernel. Operates with the types that could
-// be converted to the int32 without precision loss.
-template <>
-class voronoi_calc_kernel<int> {
-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,
- COLLINEAR = 0,
- LEFT = 1,
- };
-
- // Value is a determinant of two vectors (e.g. x1 * y2 - x2 * y1).
- // Return orientation based on the sign of the determinant.
- template <typename T>
- static kOrientation get_orientation(T value) {
- if (value == static_cast<T>(0.0)) return COLLINEAR;
- return (value < static_cast<T>(0.0)) ? RIGHT : LEFT;
- }
-
- // Compute robust cross_product: a1 * b2 - b1 * a2.
- // The result is correct with epsilon relative error equal to 2EPS.
- template <typename T>
- static fpt_type robust_cross_product(T a1_, T b1_, T a2_, T b2_) {
- 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);
-
- ulong_long_type expr_l = a1 * b2;
- bool expr_l_plus = (a1_plus == b2_plus) ? true : false;
- ulong_long_type expr_r = b1 * a2;
- bool expr_r_plus = (a2_plus == b1_plus) ? true : false;
-
- if (expr_l == 0) expr_l_plus = true;
- if (expr_r == 0) expr_r_plus = true;
-
- if ((expr_l_plus == expr_r_plus) && (expr_l == expr_r))
- return static_cast<fpt_type>(0.0);
-
- if (!expr_l_plus) {
- if (expr_r_plus)
- return -static_cast<fpt_type>(expr_l) -
- static_cast<fpt_type>(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);
- else
- return (expr_l < expr_r) ?
- -static_cast<fpt_type>(expr_r - expr_l) :
- static_cast<fpt_type>(expr_l - expr_r);
- }
- }
-
- template <typename T>
- static kOrientation get_orientation(T dif_x1_, T dif_y1_, T dif_x2_, T dif_y2_) {
- return get_orientation(robust_cross_product(dif_x1_, dif_y1_, dif_x2_, dif_y2_));
- }
-
- template <typename Point>
- static kOrientation get_orientation(const Point &point1,
- const Point &point2,
- const Point &point3) {
- fpt_type dx1 = static_cast<fpt_type>(point1.x()) -
- static_cast<fpt_type>(point2.x());
- fpt_type dx2 = static_cast<fpt_type>(point2.x()) -
- static_cast<fpt_type>(point3.x());
- fpt_type dy1 = static_cast<fpt_type>(point1.y()) -
- static_cast<fpt_type>(point2.y());
- fpt_type dy2 = static_cast<fpt_type>(point2.y()) -
- static_cast<fpt_type>(point3.y());
- return get_orientation(robust_cross_product(dx1, dy1, dx2, dy2));
- }
-
- template <typename Point>
- static bool is_vertical(const Point &point1, const Point &point2) {
- return point1.x() == point2.x();
- }
-
- template <typename Site>
- static bool is_vertical(const Site &site) {
- return is_vertical(site.point0(), site.point1());
- }
-
- template <typename Point>
- class point_comparison_predicate {
- public:
- typedef Point point_type;
-
- bool operator()(const point_type &lhs, const point_type &rhs) const {
- if (lhs.x() == rhs.x()) {
- return lhs.y() < rhs.y();
- }
- return lhs.x() < rhs.x();
- }
- };
-
- template <typename Site, typename Circle>
- class event_comparison_predicate {
- public:
- typedef Site site_type;
- typedef Circle circle_type;
-
- 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()) {
- if (!rhs.is_segment()) {
- return lhs.y0() < rhs.y0();
- }
- if (is_vertical(rhs)) {
- return lhs.y0() <= rhs.y0();
- }
- return true;
- } else {
- if (is_vertical(rhs)) {
- if(is_vertical(lhs)) {
- return lhs.y0() < rhs.y0();
- }
- return false;
- }
- if (is_vertical(lhs)) {
- return true;
- }
- if (lhs.y0() != rhs.y0()) {
- return lhs.y0() < rhs.y0();
- }
- return get_orientation(lhs.point1(), lhs.point0(), rhs.point1()) == LEFT;
- }
- }
-
- bool operator()(const site_type &lhs, const circle_type &rhs) const {
- if (almost_equal(static_cast<fpt_type>(lhs.x()),
- static_cast<fpt_type>(rhs.lower_x()), ULPS)) {
- if (almost_equal(static_cast<fpt_type>(lhs.y()),
- static_cast<fpt_type>(rhs.lower_y()), ULPS)) {
- return false;
- }
- return static_cast<fpt_type>(lhs.y()) <
- static_cast<fpt_type>(rhs.lower_y());
- }
- return static_cast<fpt_type>(lhs.x()) <
- static_cast<fpt_type>(rhs.lower_x());
- }
-
- bool operator()(const circle_type &lhs, const site_type &rhs) const {
- if (almost_equal(static_cast<fpt_type>(lhs.lower_x()),
- static_cast<fpt_type>(rhs.x()), ULPS)) {
- if (almost_equal(static_cast<fpt_type>(lhs.lower_y()),
- static_cast<fpt_type>(rhs.y()), ULPS)) {
- return false;
- }
- return static_cast<fpt_type>(lhs.lower_y()) <
- static_cast<fpt_type>(rhs.y());
- }
- return static_cast<fpt_type>(lhs.lower_x()) <
- static_cast<fpt_type>(rhs.x());
- }
-
- bool operator()(const circle_type &lhs, const circle_type &rhs) const {
- if (almost_equal(static_cast<fpt_type>(lhs.lower_x()),
- static_cast<fpt_type>(rhs.lower_x()), ULPSx2)) {
- if (almost_equal(static_cast<fpt_type>(lhs.lower_y()),
- static_cast<fpt_type>(rhs.lower_y()), ULPSx2)) {
- return false;
- }
- return lhs.lower_y() < rhs.lower_y();
- }
- return lhs.lower_x() < rhs.lower_x();
- }
- };
-
- template <typename Site>
- class distance_predicate {
- public:
- typedef Site site_type;
-
- // Returns true if a horizontal line going through a new site intersects
- // right arc at first, else returns false. If horizontal line goes
- // through intersection point of the given two arcs returns false also.
- bool operator()(const site_type& left_site,
- const site_type& right_site,
- const site_type& new_site) const {
- if (!left_site.is_segment()) {
- if (!right_site.is_segment()) {
- return pp(left_site, right_site, new_site);
- } else {
- return ps(left_site, right_site, new_site, false);
- }
- } else {
- if (!right_site.is_segment()) {
- return ps(right_site, left_site, new_site, true);
- } else {
- return ss(left_site, right_site, new_site);
- }
- }
- }
- private:
- typedef typename Site::point_type point_type;
-
- // Robust predicate, avoids using high-precision libraries.
- // Returns true if a horizontal line going through the new point site
- // intersects right arc at first, else returns false. If horizontal line
- // goes through intersection point of the given two arcs returns false.
- bool pp(const site_type &left_site,
- const site_type &right_site,
- const site_type &new_site) const {
- const point_type &left_point = left_site.point0();
- const point_type &right_point = right_site.point0();
- const point_type &new_point = new_site.point0();
- if (left_point.x() > right_point.x()) {
- if (new_point.y() <= left_point.y())
- return false;
- } else if (left_point.x() < right_point.x()) {
- if (new_point.y() >= right_point.y())
- return true;
- } else {
- return static_cast<fpt_type>(left_point.y()) +
- static_cast<fpt_type>(right_point.y()) <
- 2.0 * static_cast<fpt_type>(new_point.y());
- }
-
- fpt_type dist1 = find_distance_to_point_arc(left_site, new_point);
- fpt_type dist2 = find_distance_to_point_arc(right_site, new_point);
-
- // The undefined ulp range is equal to 3EPS + 3EPS <= 6ULP.
- return dist1 < dist2;
- }
-
- bool ps(const site_type &left_site, const site_type &right_site,
- const site_type &new_site, bool reverse_order) const {
- kPredicateResult fast_res = fast_ps(
- left_site, right_site, new_site, reverse_order);
- if (fast_res != UNDEFINED) {
- return (fast_res == LESS);
- }
-
- 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.
- return reverse_order ^ (dist1 < dist2);
- }
-
- bool ss(const site_type &left_site,
- const site_type &right_site,
- const site_type &new_site) const {
- // Handle temporary segment sites.
- if (left_site.point0() == right_site.point0() &&
- left_site.point1() == right_site.point1()) {
- return get_orientation(left_site.point0(),
- left_site.point1(),
- new_site.point0()) == LEFT;
- }
-
- 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.
- return dist1 < dist2;
- }
-
- fpt_type find_distance_to_point_arc(const site_type &site,
- const point_type &point) const {
- fpt_type dx = site.x() - point.x();
- fpt_type dy = site.y() - point.y();
- // The relative error is atmost 3EPS.
- return (dx * dx + dy * dy) / (static_cast<fpt_type>(2.0) * dx);
- }
-
- fpt_type find_distance_to_segment_arc(const site_type &site,
- const point_type &point) const {
- if (is_vertical(site)) {
- return (static_cast<fpt_type>(site.x()) - static_cast<fpt_type>(point.x())) *
- static_cast<fpt_type>(0.5);
- } else {
- const point_type &segment0 = site.point0(true);
- const point_type &segment1 = site.point1(true);
- fpt_type a1 = static_cast<fpt_type>(segment1.x()) -
- static_cast<fpt_type>(segment0.x());
- fpt_type b1 = static_cast<fpt_type>(segment1.y()) -
- static_cast<fpt_type>(segment0.y());
- fpt_type a3 = static_cast<fpt_type>(point.x()) -
- static_cast<fpt_type>(segment0.x());
- fpt_type b3 = static_cast<fpt_type>(point.y()) -
- 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.0)) {
- k = static_cast<fpt_type>(1.0) / (b1 + k);
- } else {
- k = (k - b1) / (a1 * a1);
- }
- // The relative error is atmost 8EPS.
- return robust_cross_product(a1, b1, a3, b3) * k;
- }
- }
-
- kPredicateResult fast_ps(const site_type &left_site, const site_type &right_site,
- const site_type &new_site, bool reverse_order) const {
- const point_type &site_point = left_site.point0();
- const point_type &segment_start = right_site.point0(true);
- const point_type &segment_end = right_site.point1(true);
- const point_type &new_point = new_site.point0();
- if (get_orientation(segment_start, segment_end, new_point) != RIGHT) {
- return (!right_site.is_inverse()) ? LESS : MORE;
- }
-
- fpt_type dif_x = static_cast<fpt_type>(new_point.x()) -
- static_cast<fpt_type>(site_point.x());
- fpt_type dif_y = static_cast<fpt_type>(new_point.y()) -
- static_cast<fpt_type>(site_point.y());
- fpt_type a = static_cast<fpt_type>(segment_end.x()) -
- static_cast<fpt_type>(segment_start.x());
- fpt_type b = static_cast<fpt_type>(segment_end.y()) -
- static_cast<fpt_type>(segment_start.y());
-
- if (is_vertical(right_site)) {
- if (new_point.y() < site_point.y() && !reverse_order)
- return MORE;
- else if (new_point.y() > site_point.y() && reverse_order)
- return LESS;
- return UNDEFINED;
- } else {
- kOrientation orientation = get_orientation(a, b, dif_x, dif_y);
- if (orientation == LEFT) {
- if (!right_site.is_inverse())
- return reverse_order ? LESS : UNDEFINED;
- return reverse_order ? UNDEFINED : MORE;
- }
- }
-
- fpt_type fast_left_expr = a * (dif_y + dif_x) * (dif_y - dif_x);
- fpt_type fast_right_expr = (static_cast<fpt_type>(2) * b) * dif_x * dif_y;
- if (!(almost_equal(fast_left_expr, fast_right_expr, 4))) {
- if ((fast_left_expr > fast_right_expr) ^ reverse_order)
- return reverse_order ? LESS : MORE;
- return UNDEFINED;
- }
- return UNDEFINED;
- }
- };
-
- template <typename Node>
- class node_comparison_predicate {
- public:
- typedef Node node_type;
- typedef typename Node::site_type site_type;
- typedef typename site_type::coordinate_type coordinate_type;
- typedef distance_predicate<site_type> distance_predicate_type;
-
- // Compares nodes in the balanced binary search tree. Nodes are
- // compared based on the y coordinates of the arcs intersection points.
- // Nodes with less y coordinate of the intersection point go first.
- // Comparison is only called during the new site events processing.
- // That's why one of the nodes will always lie on the sweepline and may
- // be represented as a straight horizontal line.
- bool operator() (const node_type &node1,
- const node_type &node2) const {
- // Get x coordinate of the righmost site from both nodes.
- const site_type &site1 = get_comparison_site(node1);
- const site_type &site2 = get_comparison_site(node2);
-
- if (site1.x() < site2.x()) {
- // The second node contains a new site.
- return predicate_(node1.left_site(), node1.right_site(), site2);
- } else if (site1.x() > site2.x()) {
- // The first node contains a new site.
- return !predicate_(node2.left_site(), node2.right_site(), site1);
- } else {
- // This checks were evaluated experimentally.
- if (site1.index() == site2.index()) {
- // Both nodes are new (inserted during the same site event processing).
- return get_comparison_y(node1) < get_comparison_y(node2);
- } else if (site1.index() < site2.index()) {
- std::pair<coordinate_type, int> y1 = get_comparison_y(node1, false);
- std::pair<coordinate_type, int> y2 = get_comparison_y(node2, true);
- if (y1.first != y2.first) return y1.first < y2.first;
- return (!site1.is_segment()) ? (y1.second < 0) : false;
- } else {
- std::pair<coordinate_type, int> y1 = get_comparison_y(node1, true);
- std::pair<coordinate_type, int> y2 = get_comparison_y(node2, false);
- if (y1.first != y2.first) return y1.first < y2.first;
- return (!site2.is_segment()) ? (y2.second > 0) : true;
- }
- }
- }
- private:
- // Get the newer site.
- const site_type &get_comparison_site(const node_type &node) const {
- if (node.left_site().index() > node.right_site().index()) {
- return node.left_site();
- }
- return node.right_site();
- }
-
- // Get comparison pair: y coordinate and direction of the newer site.
- std::pair<coordinate_type, int> get_comparison_y(
- const node_type &node, bool is_new_node = true) const {
- if (node.left_site().index() == node.right_site().index()) {
- return std::make_pair(node.left_site().y(), 0);
- }
- if (node.left_site().index() > node.right_site().index()) {
- if (!is_new_node &&
- node.left_site().is_segment() &&
- is_vertical(node.left_site())) {
- return std::make_pair(node.left_site().y1(), 1);
- }
- return std::make_pair(node.left_site().y(), 1);
- }
- return std::make_pair(node.right_site().y(), -1);
- }
- private:
- distance_predicate_type predicate_;
- };
-
- template <typename Site>
- class circle_existence_predicate {
- public:
- typedef typename Site::point_type point_type;
- typedef Site site_type;
-
- bool ppp(const site_type &site1,
- const site_type &site2,
- const site_type &site3) const {
- return get_orientation(site1.point0(), site2.point0(), site3.point0()) == RIGHT;
- }
-
- bool pps(const site_type &site1,
- const site_type &site2,
- const site_type &site3,
- int segment_index) const {
- if (segment_index != 2) {
- kOrientation orient1 = get_orientation(site1.point0(),
- site2.point0(), site3.point0(true));
- kOrientation orient2 = get_orientation(site1.point0(),
- site2.point0(), site3.point1(true));
- if (segment_index == 1 && site1.x0() >= site2.x0()) {
- if (orient1 != RIGHT)
- return false;
- } else if (segment_index == 3 && site2.x0() >= site1.x0()) {
- if (orient2 != RIGHT)
- return false;
- } else if (orient1 != RIGHT && orient2 != RIGHT) {
- return false;
- }
- } else {
- if (site3.point0(true) == site1.point0() &&
- site3.point1(true) == site2.point0())
- return false;
- }
- return true;
- }
-
- bool pss(const site_type &site1,
- const site_type &site2,
- const site_type &site3,
- int point_index) const {
- if (site2.point0() == site3.point0() &&
- site2.point1() == site3.point1()) {
- return false;
- }
- if (point_index == 2) {
- if (!site2.is_inverse() && site3.is_inverse())
- return false;
- if (site2.is_inverse() == site3.is_inverse() &&
- get_orientation(site2.point0(true), site1.point0(), site3.point1(true)) != RIGHT)
- return false;
- }
- return true;
- }
-
- bool sss(const site_type &site1,
- const site_type &site2,
- const site_type &site3) const {
- if (site1.point0() == site2.point0() &&
- site1.point1() == site2.point1())
- return false;
- if (site2.point0() == site3.point0() &&
- site2.point1() == site3.point1())
- return false;
- return true;
- }
- };
-
- template <typename Site, typename Circle>
- class gmp_circle_formation_functor {
- public:
- typedef typename Site::point_type point_type;
- typedef Site site_type;
- typedef Circle circle_type;
-
- typedef mpt_wrapper<mpz_class, 8> mpt_type;
- typedef mpt_wrapper<mpf_class, 2> mpf_type;
- typedef robust_sqrt_expr<mpt_type, mpf_type> robust_sqrt_expr_type;
-
- 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) {
- 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] = 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];
- mpz_numerator[2] = mpz_dif_x[1] * mpz_sum_x[1] + mpz_dif_y[1] * mpz_sum_y[1];
-
- if (recompute_c_x || recompute_lower_x) {
- mpz_c_x = mpz_numerator[1] * mpz_dif_y[1] - mpz_numerator[2] * mpz_dif_y[0];
- circle.x(mpz_c_x.get_d() / denom.get_d());
-
- if (recompute_lower_x) {
- // 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];
- 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() >= 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;
- fpt_type lower_x = mpz_numerator[0].get_d() /
- (denom.get_d() * (mpz_c_x.get_d() + r));
- circle.lower_x(lower_x);
- }
- }
- }
-
- if (recompute_c_y) {
- 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());
- }
- }
-
- // Recompute parameters of the circle event using high-precision library.
- void pps(const site_type &site1,
- const site_type &site2,
- const site_type &site3,
- int segment_index,
- circle_type &c_event,
- bool recompute_c_x = true,
- bool recompute_c_y = true,
- bool recompute_lower_x = true) {
- 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 = 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 = 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] = 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] = 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;
-
- if (denom == 0) {
- numer = teta * teta - sum_AB * sum_AB;
- denom = teta * sum_AB;
- cA[0] = denom * sum_x * 2 + numer * vec_x;
- cB[0] = segm_len;
- cA[1] = denom * sum_AB * 2 + numer * teta;
- cB[1] = 1;
- cA[2] = denom * sum_y * 2 + numer * vec_y;
- if (recompute_c_x) {
- c_event.x(0.25 * cA[0].get_d() / denom.get_d());
- }
- if (recompute_c_y) {
- c_event.y(0.25 * cA[2].get_d() / denom.get_d());
- }
- if (recompute_lower_x) {
- c_event.lower_x(0.25 * sqrt_expr_.eval2(cA, cB).get_d() /
- (denom.get_d() * std::sqrt(segm_len.get_d())));
- }
- return;
- }
-
- det = (teta * teta + denom * denom) * A * B * 4;
- 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;
- cB[0] = 1;
- cA[1] = (segment_index == 2) ? -vec_x : vec_x;
- cB[1] = det;
- if (recompute_c_x) {
- c_event.x(0.5 * sqrt_expr_.eval2(cA, cB).get_d() / denom_sqr);
- }
- }
-
- if (recompute_c_y || recompute_lower_x) {
- cA[2] = sum_y * denom * denom + teta * sum_AB * vec_y;
- cB[2] = 1;
- cA[3] = (segment_index == 2) ? -vec_y : vec_y;
- cB[3] = det;
- if (recompute_c_y) {
- c_event.y(0.5 * sqrt_expr_.eval2(&cA[2], &cB[2]).get_d() /
- denom_sqr);
- }
- }
-
- if (recompute_lower_x) {
- cB[0] *= segm_len;
- cB[1] *= segm_len;
- cA[2] = sum_AB * (denom * denom + teta * teta);
- cB[2] = 1;
- cA[3] = (segment_index == 2) ? -teta : teta;
- cB[3] = det;
- c_event.lower_x(0.5 * sqrt_expr_.eval4(cA, cB).get_d() /
- (denom_sqr * std::sqrt(segm_len.get_d())));
- }
- }
-
- // Recompute parameters of the circle event using high-precision library.
- void pss(const site_type &site1,
- const site_type &site2,
- const site_type &site3,
- int point_index,
- circle_type &c_event,
- bool recompute_c_x = true,
- bool recompute_c_y = true,
- bool recompute_lower_x = true) {
- 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] = 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) {
- 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] * (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] * (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) {
- 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];
- fpt_type lower_x = sqrt_expr_.eval3(cA, cB).get_d();
- c_event.lower_x(lower_x / denom);
- }
- }
- 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();
- ix = a[0] * c[1] + a[1] * c[0];
- iy = b[0] * c[1] + b[1] * c[0];
- dx = ix - orientation * site1.x();
- dy = iy - orientation * site1.y();
- if ((dx == 0) && (dy == 0)) {
- 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;
- }
-
- 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;
- cA[3] = 0;
- cB[0] = a[0] * a[0] + b[0] * b[0];
- 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;
- 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;
- fpt_type cy = sqrt_expr_evaluator_pss<mpt_type, mpf_type>(cA, cB).get_d();
- c_event.y(cy / denom);
- }
-
- if (recompute_c_x || recompute_lower_x) {
- cA[0] = a[1] * (dx * dx + dy * dy) - ix * (dx * a[1] + dy * b[1]);
- cA[1] = a[0] * (dx * dx + dy * dy) - ix * (dx * a[0] + dy * b[0]);
- cA[2] = ix * sign;
-
- if (recompute_c_x) {
- 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);
- fpt_type lower_x = sqrt_expr_evaluator_pss<mpt_type, mpf_type>(cA, cB).get_d();
- c_event.lower_x(lower_x / denom);
- }
- }
- }
-
- // Recompute parameters of the circle event using high-precision library.
- 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) {
- static mpt_type a[3], b[3], c[3], cA[4], cB[4];
- // cA - corresponds to the cross product.
- // cB - corresponds to the squared length.
- 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] = 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));
- c[2] = mpt_cross(site3.x0(true), site3.y0(true), site3.x1(true), site3.y1(true));
-
- for (int i = 0; i < 3; ++i) {
- cB[i] = a[i] * a[i] + b[i] * b[i];
- }
-
- for (int i = 0; i < 3; ++i) {
- int j = (i+1) % 3;
- int k = (i+2) % 3;
- cA[i] = a[j] * b[k] - a[k] * b[j];
- }
- fpt_type denom = sqrt_expr_.eval3(cA, cB).get_d();
-
- if (recompute_c_y) {
- for (int i = 0; i < 3; ++i) {
- int j = (i+1) % 3;
- int k = (i+2) % 3;
- cA[i] = b[j] * c[k] - b[k] * c[j];
- }
- fpt_type c_y = sqrt_expr_.eval3(cA, cB).get_d();
- c_event.y(c_y / denom);
- }
-
- if (recompute_c_x || recompute_lower_x) {
- cA[3] = 0;
- for (int i = 0; i < 3; ++i) {
- int j = (i+1) % 3;
- int k = (i+2) % 3;
- cA[i] = a[j] * c[k] - a[k] * c[j];
- if (recompute_lower_x) {
- cA[3] += cA[i] * b[i];
- }
- }
-
- if (recompute_c_x) {
- fpt_type c_x = sqrt_expr_.eval3(cA, cB).get_d();
- c_event.x(c_x / denom);
- }
-
- if (recompute_lower_x) {
- cB[3] = 1;
- fpt_type lower_x = sqrt_expr_.eval4(cA, cB).get_d();
- c_event.lower_x(lower_x / denom);
- }
- }
- }
- private:
- template <typename T>
- mpt_wrapper<mpz_class, 8> &mpt_cross(T a0, T b0, T a1, T b1) {
- static mpt_type temp[2];
- temp[0] = a0;
- temp[1] = b0;
- temp[0] = temp[0] * b1;
- temp[1] = temp[1] * a1;
- temp[0] -= temp[1];
- return temp[0];
- }
-
- // Evaluates A[3] + A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) +
- // A[2] * sqrt(B[3] * (sqrt(B[0] * B[1]) + B[2]));
- template <typename mpt, typename mpf>
- mpf sqrt_expr_evaluator_pss(mpt *A, mpt *B) {
- static mpt cA[4], cB[4];
- static mpf lh, rh, numer;
- if (A[3] == 0) {
- lh = sqrt_expr_.eval2(A, B);
- cA[0] = 1;
- cB[0] = B[0] * B[1];
- cA[1] = B[2];
- cB[1] = 1;
- rh = A[2].get_d() * std::sqrt(B[3].get_d() *
- sqrt_expr_.eval2(cA, cB).get_d());
- if (((lh >= 0) && (rh >= 0)) || ((lh <= 0) && (rh <= 0))) {
- return lh + rh;
- }
- cA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1];
- cA[0] -= A[2] * A[2] * B[3] * B[2];
- cB[0] = 1;
- cA[1] = A[0] * A[1] * 2 - A[2] * A[2] * B[3];
- cB[1] = B[0] * B[1];
- numer = sqrt_expr_.eval2(cA, cB);
- return numer / (lh - rh);
- }
- cA[0] = A[3];
- cB[0] = 1;
- cA[1] = A[0];
- cB[1] = B[0];
- cA[2] = A[1];
- cB[2] = B[1];
- lh = sqrt_expr_.eval3(cA, cB);
- cA[0] = 1;
- cB[0] = B[0] * B[1];
- cA[1] = B[2];
- cB[1] = 1;
- rh = A[2].get_d() * std::sqrt(B[3].get_d() *
- sqrt_expr_.eval2(cA, cB).get_d());
- if (((lh >= 0) && (rh >= 0)) || ((lh <= 0) && (rh <= 0))) {
- return lh + rh;
- }
- cA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1];
- cA[0] += A[3] * A[3] - A[2] * A[2] * B[2] * B[3];
- cB[0] = 1;
- cA[1] = A[3] * A[0] * 2;
- cB[1] = B[0];
- cA[2] = A[3] * A[1] * 2;
- cB[2] = B[1];
- cA[3] = A[0] * A[1] * 2 - A[2] * A[2] * B[3];
- cB[3] = B[0] * B[1];
- numer = sqrt_expr_.eval4(cA, cB);
- return numer / (lh - rh);
- }
- private:
- robust_sqrt_expr_type sqrt_expr_;
- };
-
- 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;
- typedef gmp_circle_formation_functor<site_type, circle_type>
- exact_circle_formation_functor_type;
-
- void ppp(const site_type &site1,
- const site_type &site2,
- const site_type &site3,
- circle_type &c_event) {
- 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() > 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) {
- exact_circle_formation_functor_.ppp(
- site1, site2, site3, c_event, recompute_c_x, recompute_c_y, recompute_lower_x);
- }
- }
-
- void pps(const site_type &site1,
- const site_type &site2,
- const site_type &site3,
- int segment_index,
- circle_type &c_event) {
- 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_type(8.0, false) * A);
- t -= A / (robust_fpt_type(2.0, false) * teta);
- } else {
- 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_type(2.0, false) * denom * denom);
- }
- 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() > 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) {
- exact_circle_formation_functor_.pps(
- site1, site2, site3, segment_index, c_event,
- recompute_c_x, recompute_c_y, recompute_lower_x);
- }
- }
-
- void pss(const site_type &site1,
- const site_type &site2,
- const site_type &site3,
- int point_index,
- circle_type &c_event) {
- 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);
- 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_type orientation(robust_cross_product(b1, a1, b2, a2), 2.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);
- 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_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_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_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();
- } else {
- t -= det.sqrt();
- }
- t /= (a * a);
- 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_type lower_x(c_x);
- lower_x += t * orientation.fabs();
- 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) {
- exact_circle_formation_functor_.pss(
- site1, site2, site3, point_index, c_event,
- recompute_c_x, recompute_c_y, recompute_lower_x);
- }
- }
-
- void sss(const site_type &site1,
- const site_type &site2,
- const site_type &site3,
- circle_type &c_event) {
- 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_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 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 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;
- 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;
- c_x -= a3 * c2 * len1;
- c_x += a3 * c1 * len2;
- c_x -= a1 * c3 * len2;
- c_y += b1 * c2 * len3;
- c_y -= b2 * c1 * len3;
- c_y += b2 * c3 * len1;
- c_y -= b3 * c2 * len1;
- c_y += b3 * c1 * len2;
- c_y -= b1 * c3 * len2;
- 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());
- if (recompute_c_x || recompute_c_y || recompute_lower_x || recompute_denom) {
- exact_circle_formation_functor_.sss(
- site1, site2, site3, c_event,
- recompute_c_x, recompute_c_y, recompute_lower_x);
- }
- }
- private:
- template <typename T>
- T sqr_distance(T dif_x, T dif_y) {
- return dif_x * dif_x + dif_y * dif_y;
- }
- private:
- exact_circle_formation_functor_type exact_circle_formation_functor_;
- };
-
- template <typename Site,
- typename Circle,
- typename CEP = circle_existence_predicate<Site>,
- typename CFF = lazy_circle_formation_functor<Site, Circle> >
- class circle_formation_predicate {
- public:
- typedef Site site_type;
- typedef Circle circle_type;
- typedef CEP circle_existence_predicate_type;
- typedef CFF circle_formation_functor_type;
-
- // Create a circle event from the given three sites.
- // Returns true if the circle event exists, else false.
- // If exists circle event is saved into the c_event variable.
- bool operator()(const site_type &site1, const site_type &site2,
- const site_type &site3, circle_type &circle) {
- if (!site1.is_segment()) {
- if (!site2.is_segment()) {
- if (!site3.is_segment()) {
- // (point, point, point) sites.
- if (!circle_existence_predicate_.ppp(site1, site2, site3))
- return false;
- circle_formation_functor_.ppp(site1, site2, site3, circle);
- } else {
- // (point, point, segment) sites.
- if (!circle_existence_predicate_.pps(site1, site2, site3, 3))
- return false;
- circle_formation_functor_.pps(site1, site2, site3, 3, circle);
- }
- } else {
- if (!site3.is_segment()) {
- // (point, segment, point) sites.
- if (!circle_existence_predicate_.pps(site1, site3, site2, 2))
- return false;
- circle_formation_functor_.pps(site1, site3, site2, 2, circle);
- } else {
- // (point, segment, segment) sites.
- if (!circle_existence_predicate_.pss(site1, site2, site3, 1))
- return false;
- circle_formation_functor_.pss(site1, site2, site3, 1, circle);
- }
- }
- } else {
- if (!site2.is_segment()) {
- if (!site3.is_segment()) {
- // (segment, point, point) sites.
- if (!circle_existence_predicate_.pps(site2, site3, site1, 1))
- return false;
- circle_formation_functor_.pps(site2, site3, site1, 1, circle);
- } else {
- // (segment, point, segment) sites.
- if (!circle_existence_predicate_.pss(site2, site1, site3, 2))
- return false;
- circle_formation_functor_.pss(site2, site1, site3, 2, circle);
- }
- } else {
- if (!site3.is_segment()) {
- // (segment, segment, point) sites.
- if (!circle_existence_predicate_.pss(site3, site1, site2, 3))
- return false;
- circle_formation_functor_.pss(site3, site1, site2, 3, circle);
- } else {
- // (segment, segment, segment) sites.
- if (!circle_existence_predicate_.sss(site1, site2, site3))
- return false;
- circle_formation_functor_.sss(site1, site2, site3, circle);
- }
- }
- }
- return true;
- }
- private:
- circle_existence_predicate_type circle_existence_predicate_;
- circle_formation_functor_type circle_formation_functor_;
- };
-private:
- // Convert value to 64-bit unsigned integer.
- // Return true if the value is positive, else false.
- template <typename T>
- static bool convert_to_65_bit(T value, ulong_long_type &res) {
- if (value >= static_cast<T>(0)) {
- res = static_cast<ulong_long_type>(value);
- return true;
- } else {
- res = static_cast<ulong_long_type>(-value);
- return false;
- }
- }
-};
-
-const voronoi_calc_kernel<int>::fpt_type voronoi_calc_kernel<int>::fULPS =
- voronoi_calc_kernel<int>::ULPS;
-const voronoi_calc_kernel<int>::fpt_type voronoi_calc_kernel<int>::fULPSx2 =
- voronoi_calc_kernel<int>::ULPSx2;
-
-} // detail
-} // polygon
-} // boost
-
-#endif
Copied: sandbox/gtl/boost/polygon/detail/voronoi_calc_utils.hpp (from r75182, /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_utils.hpp 2011-11-08 09:16:57 EST (Tue, 08 Nov 2011)
@@ -1,4 +1,4 @@
-// Boost.Polygon library detail/voronoi_calc_kernel.hpp header file
+// Boost.Polygon library detail/voronoi_calc_utils.hpp header file
// Copyright Andrii Sydorchuk 2010-2011.
// Distributed under the Boost Software License, Version 1.0.
@@ -7,26 +7,26 @@
// See http://www.boost.org for updates, documentation, and revision history.
-#ifndef BOOST_POLYGON_VORONOI_CALC_KERNEL
-#define BOOST_POLYGON_VORONOI_CALC_KERNEL
+#ifndef BOOST_POLYGON_VORONOI_CALC_UTILS
+#define BOOST_POLYGON_VORONOI_CALC_UTILS
#pragma warning(disable:4800)
#include <gmpxx.h>
#include "mpt_wrapper.hpp"
-#include "voronoi_fpt_kernel.hpp"
+#include "voronoi_robust_fpt.hpp"
namespace boost {
namespace polygon {
namespace detail {
template <typename T>
-class voronoi_calc_kernel;
+class voronoi_calc_utils;
-// Predicates kernel. Operates with the types that could
+// Predicate utils. Operates with the types that could
// be converted to the int32 without precision loss.
template <>
-class voronoi_calc_kernel<int> {
+class voronoi_calc_utils<int> {
public:
typedef double fpt_type;
typedef unsigned long long ulong_long_type;
@@ -1394,10 +1394,10 @@
}
};
-const voronoi_calc_kernel<int>::fpt_type voronoi_calc_kernel<int>::fULPS =
- voronoi_calc_kernel<int>::ULPS;
-const voronoi_calc_kernel<int>::fpt_type voronoi_calc_kernel<int>::fULPSx2 =
- voronoi_calc_kernel<int>::ULPSx2;
+const voronoi_calc_utils<int>::fpt_type voronoi_calc_utils<int>::fULPS =
+ voronoi_calc_utils<int>::ULPS;
+const voronoi_calc_utils<int>::fpt_type voronoi_calc_utils<int>::fULPSx2 =
+ voronoi_calc_utils<int>::ULPSx2;
} // detail
} // polygon
Deleted: sandbox/gtl/boost/polygon/detail/voronoi_fpt_kernel.hpp
==============================================================================
--- sandbox/gtl/boost/polygon/detail/voronoi_fpt_kernel.hpp 2011-11-08 09:16:57 EST (Tue, 08 Nov 2011)
+++ (empty file)
@@ -1,618 +0,0 @@
-// Boost.Polygon library detail/voronoi_fpt_kernel.hpp header file
-
-// Copyright Andrii Sydorchuk 2010-2011.
-// 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_POLYGON_VORONOI_FPT_KERNEL
-#define BOOST_POLYGON_VORONOI_FPT_KERNEL
-
-// Geometry predicates with floating-point variables usually require
-// high-precision predicates to retrieve the correct result.
-// Epsilon robust predicates give the result within some epsilon relative
-// error, but are a lot faster than high-precision predicates.
-// To make algorithm robust and efficient epsilon robust predicates are
-// used at the first step. In case of the undefined result high-precision
-// arithmetic is used to produce required robustness. This approach
-// requires exact computation of epsilon intervals within which epsilon
-// robust predicates have undefined value.
-// There are two ways to measure an error of floating-point calculations:
-// relative error and ULPs (units in the last place).
-// Let EPS be machine epsilon, then next inequalities have place:
-// 1 EPS <= 1 ULP <= 2 EPS (1), 0.5 ULP <= 1 EPS <= 1 ULP (2).
-// ULPs are good for measuring rounding errors and comparing values.
-// Relative erros are good for computation of general relative
-// errors of formulas or expressions. So to calculate epsilon
-// intervals within which epsilon robust predicates have undefined result
-// next schema is used:
-// 1) Compute rounding errors of initial variables using ULPs;
-// 2) Transform ULPs to epsilons using upper bound of the (1);
-// 3) Compute relative error of the formula using epsilon arithmetic;
-// 4) Transform epsilon to ULPs using upper bound of the (2);
-// In case two values are inside undefined ULP range use high-precision
-// arithmetic to produce the correct result, else output the result.
-// Look at almost_equal function to see how two floating-point variables
-// are checked to fit in the ULP range.
-// If A has relative error of r(A) and B has relative error of r(B) then:
-// 1) r(A + B) <= max(r(A), r(B)), for A * B >= 0;
-// 2) r(A - B) <= B*r(A)+A*r(B)/(A-B), for A * B >= 0;
-// 2) r(A * B) <= r(A) + r(B);
-// 3) r(A / B) <= r(A) + r(B);
-// In addition rounding error should be added, that is always equal to
-// 0.5 ULP or atmost 1 epsilon. As you might see from the above formulas
-// substraction relative error may be extremely large, that's why
-// epsilon robust comparator class is used to store floating point values
-// and avoid substraction.
-// For further information about relative errors and ULPs try this link:
-// http://docs.sun.com/source/806-3568/ncg_goldberg.html
-
-namespace boost {
-namespace polygon {
-namespace detail {
- // Represents the result of the epsilon robust predicate.
- // If the result is undefined some further processing is usually required.
- enum kPredicateResult {
- LESS = -1,
- UNDEFINED = 0,
- MORE = 1,
- };
-
- // If two floating-point numbers in the same format are ordered (x < y),
- // then they are ordered the same way when their bits are reinterpreted as
- // sign-magnitude integers. Values are considered to be almost equal if
- // their integer reinterpretations differ in not more than maxUlps units.
- template <typename T>
- bool almost_equal(T a, T b, unsigned int ulps);
-
- template<>
- bool almost_equal<float>(float a, float b, unsigned int maxUlps) {
- unsigned int ll_a, ll_b;
-
- // Reinterpret double bits as 32-bit signed integer.
- memcpy(&ll_a, &a, sizeof(float));
- memcpy(&ll_b, &b, sizeof(float));
-
- if (ll_a < 0x80000000)
- ll_a = 0x80000000 - ll_a;
- if (ll_b < 0x80000000)
- ll_b = 0x80000000 - ll_b;
-
- if (ll_a > ll_b)
- return ll_a - ll_b <= maxUlps;
- return ll_b - ll_a <= maxUlps;
- }
-
- template<>
- bool almost_equal<double>(double a, double b, unsigned int maxUlps) {
- unsigned long long ll_a, ll_b;
-
- // Reinterpret double bits as 64-bit signed integer.
- memcpy(&ll_a, &a, sizeof(double));
- memcpy(&ll_b, &b, sizeof(double));
-
- // Positive 0.0 is integer zero. Negative 0.0 is 0x8000000000000000.
- // Map negative zero to an integer zero representation - making it
- // identical to positive zero - the smallest negative number is
- // represented by negative one, and downwards from there.
- if (ll_a < 0x8000000000000000ULL)
- ll_a = 0x8000000000000000ULL - ll_a;
- if (ll_b < 0x8000000000000000ULL)
- ll_b = 0x8000000000000000ULL - ll_b;
-
- // Compare 64-bit signed integer representations of input values.
- // Difference in 1 Ulp is equivalent to a relative error of between
- // 1/4,000,000,000,000,000 and 1/8,000,000,000,000,000.
- if (ll_a > ll_b)
- return ll_a - ll_b <= maxUlps;
- return ll_b - ll_a <= maxUlps;
- }
-
- template <typename T>
- double get_d(const T& value) {
- return value.get_d();
- }
-
- template <>
- double get_d(const float& value) {
- return value;
- }
-
- template <>
- double get_d(const double& value) {
- return value;
- }
-
- template <>
- double get_d(const long double& value) {
- return static_cast<double>(value);
- }
-
- template <typename FPT>
- class robust_fpt {
- public:
- typedef FPT floating_point_type;
- typedef FPT relative_error_type;
-
- // Rounding error is at most 1 EPS.
- static const relative_error_type ROUNDING_ERROR;
-
- robust_fpt() : fpv_(0.0), re_(0) {}
- explicit robust_fpt(int fpv) : fpv_(fpv), re_(0) {}
- explicit robust_fpt(floating_point_type fpv,
- bool rounded = true) : fpv_(fpv) {
- re_ = rounded ? ROUNDING_ERROR : 0;
- }
- robust_fpt(floating_point_type fpv, relative_error_type error) :
- fpv_(fpv), re_(error) {}
-
- floating_point_type fpv() const { return fpv_; }
- relative_error_type re() const { return re_; }
- relative_error_type ulp() const { return re_; }
-
- template <typename T>
- bool operator==(T that) const {
- floating_point_type value = static_cast<floating_point_type>(that);
- return almost_equal(this->fpv_,
- value,
- static_cast<unsigned int>(this->ulp()));
- }
-
- template <typename T>
- bool operator!=(T that) const {
- return !(*this == that);
- }
-
- template <typename T>
- bool operator<(T that) const {
- if (*this == that) return false;
- return this->fpv_ < static_cast<floating_point_type>(that);
- }
-
- template <typename T>
- bool operator<=(T that) const {
- if (*this == that) return true;
- return this->fpv_ < static_cast<floating_point_type>(that);
- }
-
- template <typename T>
- bool operator>(T that) const {
- if (*this == that) return false;
- return this->fpv_ > static_cast<floating_point_type>(that);
- }
-
- template <typename T>
- bool operator>=(T that) const {
- if (*this == that) return true;
- return this->fpv_ > static_cast<floating_point_type>(that);
- }
-
- bool operator==(const robust_fpt &that) const {
- unsigned int ulp = static_cast<unsigned int>(this->re_ + that.re_);
- return almost_equal(this->fpv_, that.fpv_, ulp);
- }
-
- bool operator!=(const robust_fpt &that) const {
- return !(*this == that);
- }
-
- bool operator<(const robust_fpt &that) const {
- if (*this == that)
- return false;
- return this->fpv_ < that.fpv_;
- }
-
- bool operator>(const robust_fpt &that) const {
- return that < *this;
- }
-
- bool operator<=(const robust_fpt &that) const {
- return !(that < *this);
- }
-
- bool operator>=(const robust_fpt &that) const {
- return !(*this < that);
- }
-
- robust_fpt operator-() const {
- return robust_fpt(-fpv_, re_);
- }
-
- robust_fpt& operator=(const robust_fpt &that) {
- this->fpv_ = that.fpv_;
- this->re_ = that.re_;
- return *this;
- }
-
- robust_fpt& operator+=(const robust_fpt &that) {
- floating_point_type fpv = this->fpv_ + that.fpv_;
- if ((this->fpv_ >= 0 && that.fpv_ >= 0) ||
- (this->fpv_ <= 0 && that.fpv_ <= 0))
- this->re_ = (std::max)(this->re_, that.re_) + ROUNDING_ERROR;
- else {
- floating_point_type temp =
- (this->fpv_ * this->re_ - that.fpv_ * that.re_) / fpv;
- this->re_ = std::fabs(temp) + ROUNDING_ERROR;
- }
- this->fpv_ = fpv;
- return *this;
- }
-
- robust_fpt& operator-=(const robust_fpt &that) {
- floating_point_type fpv = this->fpv_ - that.fpv_;
- if ((this->fpv_ >= 0 && that.fpv_ <= 0) ||
- (this->fpv_ <= 0 && that.fpv_ >= 0))
- this->re_ = (std::max)(this->re_, that.re_) + ROUNDING_ERROR;
- else {
- floating_point_type temp =
- (this->fpv_ * this->re_ + that.fpv_ * that.re_) / fpv;
- this->re_ = std::fabs(temp) + ROUNDING_ERROR;
- }
- this->fpv_ = fpv;
- return *this;
- }
-
- robust_fpt& operator*=(const robust_fpt &that) {
- this->re_ += that.re_ + ROUNDING_ERROR;
- this->fpv_ *= that.fpv_;
- return *this;
- }
-
- robust_fpt& operator/=(const robust_fpt &that) {
- this->re_ += that.re_ + ROUNDING_ERROR;
- this->fpv_ /= that.fpv_;
- return *this;
- }
-
- robust_fpt operator+(const robust_fpt &that) const {
- floating_point_type fpv = this->fpv_ + that.fpv_;
- relative_error_type re;
- if ((this->fpv_ >= 0 && that.fpv_ >= 0) ||
- (this->fpv_ <= 0 && that.fpv_ <= 0))
- re = (std::max)(this->re_, that.re_) + ROUNDING_ERROR;
- else {
- floating_point_type temp =
- (this->fpv_ * this->re_ - that.fpv_ * that.re_) / fpv;
- re = std::fabs(temp) + ROUNDING_ERROR;
- }
- return robust_fpt(fpv, re);
- }
-
- robust_fpt operator-(const robust_fpt &that) const {
- floating_point_type fpv = this->fpv_ - that.fpv_;
- relative_error_type re;
- if ((this->fpv_ >= 0 && that.fpv_ <= 0) ||
- (this->fpv_ <= 0 && that.fpv_ >= 0))
- re = (std::max)(this->re_, that.re_) + ROUNDING_ERROR;
- else {
- floating_point_type temp =
- (this->fpv_ * this->re_ + that.fpv_ * that.re_) / fpv;
- re = std::fabs(temp) + ROUNDING_ERROR;
- }
- return robust_fpt(fpv, re);
- }
-
- robust_fpt operator*(const robust_fpt &that) const {
- floating_point_type fpv = this->fpv_ * that.fpv_;
- relative_error_type re = this->re_ + that.re_ + ROUNDING_ERROR;
- return robust_fpt(fpv, re);
- }
-
- robust_fpt operator/(const robust_fpt &that) const {
- floating_point_type fpv = this->fpv_ / that.fpv_;
- relative_error_type re = this->re_ + that.re_ + ROUNDING_ERROR;
- return robust_fpt(fpv, re);
- }
-
- robust_fpt sqrt() const {
- return robust_fpt(std::sqrt(fpv_), re_ * 0.5 + ROUNDING_ERROR);
- }
-
- robust_fpt fabs() const {
- return (fpv_ >= 0) ? *this : -(*this);
- }
- private:
- floating_point_type fpv_;
- relative_error_type re_;
- };
-
- template <typename T>
- const typename robust_fpt<T>::relative_error_type
- robust_fpt<T>::ROUNDING_ERROR = 1;
-
- // robust_dif consists of two not negative values: value1 and value2.
- // The resulting expression is equal to the value1 - value2.
- // Substraction of a positive value is equivalent to the addition to value2
- // and substraction of a negative value is equivalent to the addition to
- // value1. The structure implicitly avoids difference computation.
- template <typename T>
- class robust_dif {
- public:
- robust_dif() :
- positive_sum_(0),
- negative_sum_(0) {}
-
- robust_dif(const T &value) :
- positive_sum_((value>0)?value:0),
- negative_sum_((value<0)?-value:0) {}
-
- robust_dif(const T &pos, const T &neg) :
- positive_sum_(pos),
- negative_sum_(neg) {}
-
- T dif() const {
- return positive_sum_ - negative_sum_;
- }
-
- T pos() const {
- return positive_sum_;
- }
-
- T neg() const {
- return negative_sum_;
- }
-
- // Equivalent to the unary minus.
- void swap() {
- (std::swap)(positive_sum_, negative_sum_);
- }
-
- bool abs() {
- if (positive_sum_ < negative_sum_) {
- swap();
- return true;
- }
- return false;
- }
-
- robust_dif<T> operator-() const {
- return robust_dif(negative_sum_, positive_sum_);
- }
-
- robust_dif<T> &operator+=(const T &val) {
- if (val >= 0)
- positive_sum_ += val;
- else
- negative_sum_ -= val;
- return *this;
- }
-
- robust_dif<T> &operator+=(const robust_dif<T> &that) {
- positive_sum_ += that.positive_sum_;
- negative_sum_ += that.negative_sum_;
- return *this;
- }
-
- robust_dif<T> &operator-=(const T &val) {
- if (val >= 0)
- negative_sum_ += val;
- else
- positive_sum_ -= val;
- return *this;
- }
-
- robust_dif<T> &operator-=(const robust_dif<T> &that) {
- positive_sum_ += that.negative_sum_;
- negative_sum_ += that.positive_sum_;
- return *this;
- }
-
- robust_dif<T> &operator*=(const T &val) {
- if (val >= 0) {
- positive_sum_ *= val;
- negative_sum_ *= val;
- } else {
- positive_sum_ *= -val;
- negative_sum_ *= -val;
- swap();
- }
- return *this;
- }
-
- robust_dif<T> &operator*=(const robust_dif<T> &that) {
- T positive_sum = this->positive_sum_ * that.positive_sum_ +
- this->negative_sum_ * that.negative_sum_;
- T negative_sum = this->positive_sum_ * that.negative_sum_ +
- this->negative_sum_ * that.positive_sum_;
- positive_sum_ = positive_sum;
- negative_sum_ = negative_sum;
- return *this;
- }
-
- robust_dif<T> &operator/=(const T &val) {
- if (val >= 0) {
- positive_sum_ /= val;
- negative_sum_ /= val;
- } else {
- positive_sum_ /= -val;
- negative_sum_ /= -val;
- swap();
- }
- return *this;
- }
- private:
- T positive_sum_;
- T negative_sum_;
- };
-
- template<typename T>
- robust_dif<T> operator+(const robust_dif<T>& lhs,
- const robust_dif<T>& rhs) {
- return robust_dif<T>(lhs.pos() + rhs.pos(),
- lhs.neg() + rhs.neg());
- }
-
- template<typename T>
- robust_dif<T> operator+(const robust_dif<T>& lhs, const T& rhs) {
- if (rhs >= 0) {
- return robust_dif<T>(lhs.pos() + rhs, lhs.neg());
- } else {
- return robust_dif<T>(lhs.pos(), lhs.neg() - rhs);
- }
- }
-
- template<typename T>
- robust_dif<T> operator+(const T& lhs, const robust_dif<T>& rhs) {
- if (lhs >= 0) {
- return robust_dif<T>(lhs + rhs.pos(), rhs.neg());
- } else {
- return robust_dif<T>(rhs.pos(), rhs.neg() - lhs);
- }
- }
-
- template<typename T>
- robust_dif<T> operator-(const robust_dif<T>& lhs,
- const robust_dif<T>& rhs) {
- return robust_dif<T>(lhs.pos() + rhs.neg(), lhs.neg() + rhs.pos());
- }
-
- template<typename T>
- robust_dif<T> operator-(const robust_dif<T>& lhs, const T& rhs) {
- if (rhs >= 0) {
- return robust_dif<T>(lhs.pos(), lhs.neg() + rhs);
- } else {
- return robust_dif<T>(lhs.pos() - rhs, lhs.neg());
- }
- }
-
- template<typename T>
- robust_dif<T> operator-(const T& lhs, const robust_dif<T>& rhs) {
- if (lhs >= 0) {
- return robust_dif<T>(lhs + rhs.neg(), rhs.pos());
- } else {
- return robust_dif<T>(rhs.neg(), rhs.pos() - lhs);
- }
- }
-
- template<typename T>
- robust_dif<T> operator*(const robust_dif<T>& lhs,
- const robust_dif<T>& rhs) {
- T res_pos = lhs.pos() * rhs.pos() + lhs.neg() * rhs.neg();
- T res_neg = lhs.pos() * rhs.neg() + lhs.neg() * rhs.pos();
- return robust_dif<T>(res_pos, res_neg);
- }
-
- template<typename T>
- robust_dif<T> operator*(const robust_dif<T>& lhs, const T& val) {
- if (val >= 0) {
- return robust_dif<T>(lhs.pos() * val, lhs.neg() * val);
- } else {
- return robust_dif<T>(-lhs.neg() * val, -lhs.pos() * val);
- }
- }
-
- template<typename T>
- robust_dif<T> operator*(const T& val, const robust_dif<T>& rhs) {
- if (val >= 0) {
- return robust_dif<T>(val * rhs.pos(), val * rhs.neg());
- } else {
- return robust_dif<T>(-val * rhs.neg(), -val * rhs.pos());
- }
- }
-
- template<typename T>
- robust_dif<T> operator/(const robust_dif<T>& lhs, const T& val) {
- if (val >= 0) {
- return robust_dif<T>(lhs.pos() / val, lhs.neg() / val);
- } else {
- return robust_dif<T>(-lhs.neg() / val, -lhs.pos() / val);
- }
- }
-
- // Used to compute expressions that operate with sqrts with predefined
- // relative error. Evaluates expressions of the next type:
- // sum(i = 1 .. n)(A[i] * sqrt(B[i])), 1 <= n <= 4.
- template <typename mpt, typename mpf>
- class robust_sqrt_expr {
- public:
- // Evaluates expression (re = 4 EPS):
- // A[0] * sqrt(B[0]).
- mpf& eval1(mpt *A, mpt *B) {
- a[0] = A[0];
- b[0] = B[0];
- b[0] = sqrt(b[0]);
- return b[0] *= a[0];
- }
-
- // Evaluates expression (re = 7 EPS):
- // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]).
- mpf& eval2(mpt *A, mpt *B) {
- a[1] = eval1(A, B);
- b[1] = eval1(A + 1, B + 1);
- if ((a[1] >= 0 && b[1] >= 0) || (a[1] <= 0 && b[1] <= 0))
- return b[1] += a[1];
- for (int i = 0; i < 2; ++i) {
- temp[i] = A[i] * A[i];
- temp[i] *= B[i];
- }
- a[1] -= b[1];
- b[1] = temp[0] - temp[1];
- return b[1] /= a[1];
- }
-
- // Evaluates expression (re = 16 EPS):
- // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) + A[2] * sqrt(B[2]).
- mpf& eval3(mpt *A, mpt *B) {
- a[2] = eval2(A, B);
- b[2] = eval1(A + 2, B + 2);
- if ((a[2] >= 0 && b[2] >= 0) || (a[2] <= 0 && b[2] <= 0))
- return b[2] += a[2];
- for (int i = 0; i < 3; ++i) {
- temp[i] = A[i] * A[i];
- temp[i] *= B[i];
- }
- cA[0] = temp[0] + temp[1];
- cA[0] -= temp[2];
- cB[0] = 1;
- cA[1] = A[0] * A[1];
- cA[1] += cA[1];
- cB[1] = B[0] * B[1];
- a[2] -= b[2];
- b[2] = eval2(cA, cB);
- return b[2] /= a[2];
- }
-
-
- // Evaluates expression (re = 25 EPS):
- // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) +
- // A[2] * sqrt(B[2]) + A[3] * sqrt(B[3]).
- mpf& eval4(mpt *A, mpt *B) {
- a[3] = eval2(A, B);
- b[3] = eval2(A + 2, B + 2);
- if ((a[3] >= 0 && b[3] >= 0) || (a[3] <= 0 && b[3] <= 0))
- return b[3] += a[3];
- for (int i = 0; i < 4; ++i) {
- temp[i] = A[i] * A[i];
- temp[i] *= B[i];
- }
- dA[0] = temp[0] + temp[1];
- dA[0] -= temp[2];
- dA[0] -= temp[3];
- dB[0] = 1;
- dA[1] = A[0] * A[1];
- dA[1] += dA[1];
- dB[1] = B[0] * B[1];
- dA[2] = A[2] * A[3];
- dA[2] = -dA[2];
- dA[2] += dA[2];
- dB[2] = B[2] * B[3];
- a[3] -= b[3];
- b[3] = eval3(dA, dB);
- return b[3] /= a[3];
- }
- private:
- mpf a[4];
- mpf b[4];
- mpt cA[2];
- mpt cB[2];
- mpt dA[3];
- mpt dB[3];
- mpt temp[4];
- };
-} // detail
-} // polygon
-} // boost
-
-#endif
Copied: sandbox/gtl/boost/polygon/detail/voronoi_robust_fpt.hpp (from r75144, /sandbox/gtl/boost/polygon/detail/voronoi_fpt_kernel.hpp)
==============================================================================
--- /sandbox/gtl/boost/polygon/detail/voronoi_fpt_kernel.hpp (original)
+++ sandbox/gtl/boost/polygon/detail/voronoi_robust_fpt.hpp 2011-11-08 09:16:57 EST (Tue, 08 Nov 2011)
@@ -1,4 +1,4 @@
-// Boost.Polygon library detail/voronoi_fpt_kernel.hpp header file
+// Boost.Polygon library detail/voronoi_robust_fpt.hpp header file
// Copyright Andrii Sydorchuk 2010-2011.
// Distributed under the Boost Software License, Version 1.0.
@@ -7,8 +7,8 @@
// See http://www.boost.org for updates, documentation, and revision history.
-#ifndef BOOST_POLYGON_VORONOI_FPT_KERNEL
-#define BOOST_POLYGON_VORONOI_FPT_KERNEL
+#ifndef BOOST_POLYGON_VORONOI_ROBUST_FPT
+#define BOOST_POLYGON_VORONOI_ROBUST_FPT
// Geometry predicates with floating-point variables usually require
// high-precision predicates to retrieve the correct result.
@@ -25,8 +25,8 @@
// 1 EPS <= 1 ULP <= 2 EPS (1), 0.5 ULP <= 1 EPS <= 1 ULP (2).
// ULPs are good for measuring rounding errors and comparing values.
// Relative erros are good for computation of general relative
-// errors of formulas or expressions. So to calculate epsilon
-// intervals within which epsilon robust predicates have undefined result
+// error of formulas or expressions. So to calculate epsilon
+// interval within which epsilon robust predicates have undefined result
// next schema is used:
// 1) Compute rounding errors of initial variables using ULPs;
// 2) Transform ULPs to epsilons using upper bound of the (1);
@@ -45,7 +45,7 @@
// 0.5 ULP or atmost 1 epsilon. As you might see from the above formulas
// substraction relative error may be extremely large, that's why
// epsilon robust comparator class is used to store floating point values
-// and avoid substraction.
+// and compute substraction as the final step of the evaluation.
// For further information about relative errors and ULPs try this link:
// http://docs.sun.com/source/806-3568/ncg_goldberg.html
@@ -139,8 +139,8 @@
// Rounding error is at most 1 EPS.
static const relative_error_type ROUNDING_ERROR;
- robust_fpt() : fpv_(0.0), re_(0) {}
- explicit robust_fpt(int fpv) : fpv_(fpv), re_(0) {}
+ robust_fpt() : fpv_(0.0), re_(0.0) {}
+ explicit robust_fpt(int fpv) : fpv_(fpv), re_(0.0) {}
explicit robust_fpt(floating_point_type fpv,
bool rounded = true) : fpv_(fpv) {
re_ = rounded ? ROUNDING_ERROR : 0;
@@ -313,6 +313,7 @@
robust_fpt fabs() const {
return (fpv_ >= 0) ? *this : -(*this);
}
+
private:
floating_point_type fpv_;
relative_error_type re_;
@@ -432,6 +433,7 @@
}
return *this;
}
+
private:
T positive_sum_;
T negative_sum_;
@@ -602,6 +604,7 @@
b[3] = eval3(dA, dB);
return b[3] /= a[3];
}
+
private:
mpf a[4];
mpf b[4];
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-11-08 09:16:57 EST (Tue, 08 Nov 2011)
@@ -54,6 +54,7 @@
y_ = y;
return *this;
}
+
private:
coordinate_type x_;
coordinate_type y_;
@@ -199,6 +200,7 @@
bool is_inverse() const {
return is_inverse_;
}
+
private:
point_type point0_;
point_type point1_;
@@ -215,7 +217,7 @@
// 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.
- // is_active_ - flag whether circle event is still active.
+ // is_active_ - states whether circle event is still active.
// NOTE: lower_y coordinate is always equal to center_y.
template <typename T>
class circle_event {
@@ -271,6 +273,7 @@
is_active_ = false;
return *this;
}
+
private:
coordinate_type center_x_;
coordinate_type center_y_;
@@ -309,15 +312,16 @@
c_.push(c_list_.begin());
return c_list_.front();
}
+
private:
typedef typename std::list<T>::iterator list_iterator_type;
struct comparison {
- Predicate cmp_;
bool operator() (const list_iterator_type &it1,
const list_iterator_type &it2) const {
return cmp_(*it1, *it2);
}
+ Predicate cmp_;
};
std::priority_queue< list_iterator_type,
@@ -338,7 +342,7 @@
// In the 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 (e.g. has greater index).
+ // processed by the algorithm later (has greater index).
template <typename Site>
class beach_line_node_key {
public:
@@ -382,6 +386,7 @@
right_site_ = site;
return *this;
}
+
private:
site_type left_site_;
site_type right_site_;
@@ -415,6 +420,7 @@
edge_ = new_edge;
return *this;
}
+
private:
Circle *circle_event_;
Edge *edge_;
Modified: sandbox/gtl/boost/polygon/voronoi_builder.hpp
==============================================================================
--- sandbox/gtl/boost/polygon/voronoi_builder.hpp (original)
+++ sandbox/gtl/boost/polygon/voronoi_builder.hpp 2011-11-08 09:16:57 EST (Tue, 08 Nov 2011)
@@ -20,8 +20,7 @@
#include "polygon.hpp"
-#include "detail/voronoi_calc_kernel.hpp"
-#include "detail/voronoi_fpt_kernel.hpp"
+#include "detail/voronoi_calc_utils.hpp"
#include "detail/voronoi_structures.hpp"
namespace boost {
@@ -155,26 +154,26 @@
site_events_.clear();
}
private:
- typedef detail::voronoi_calc_kernel<T> CKT;
- typedef typename CKT::fpt_type coordinate_type;
+ typedef detail::voronoi_calc_utils<T> VCU;
+ typedef typename VCU::fpt_type coordinate_type;
typedef detail::point_2d<coordinate_type> point_type;
typedef detail::site_event<coordinate_type> site_event_type;
typedef typename std::vector<site_event_type>::const_iterator
site_event_iterator_type;
typedef detail::circle_event<coordinate_type> circle_event_type;
- typedef typename CKT::template point_comparison_predicate<point_type>
+ typedef typename VCU::template point_comparison_predicate<point_type>
point_comparison_predicate;
- typedef typename CKT::
+ typedef typename VCU::
template event_comparison_predicate<site_event_type, circle_event_type>
event_comparison_predicate;
- typedef typename CKT::
+ typedef typename VCU::
template circle_formation_predicate<site_event_type, circle_event_type>
circle_formation_predicate_type;
typedef typename output_type::voronoi_edge_type edge_type;
typedef detail::beach_line_node_key<site_event_type> key_type;
typedef detail::beach_line_node_data<edge_type, circle_event_type> value_type;
- typedef typename CKT::
+ typedef typename VCU::
template node_comparison_predicate<key_type> node_comparer_type;
typedef std::map< key_type, value_type, node_comparer_type > beach_line_type;
typedef typename beach_line_type::iterator beach_line_iterator;
@@ -220,9 +219,9 @@
// Count the number of the sites to init the beach line.
while(site_event_iterator_ != site_events_.end() &&
- CKT::is_vertical(site_event_iterator_->point0(),
+ VCU::is_vertical(site_event_iterator_->point0(),
site_events_.begin()->point0()) &&
- CKT::is_vertical(*site_event_iterator_)) {
+ VCU::is_vertical(*site_event_iterator_)) {
++site_event_iterator_;
++skip;
}
Modified: sandbox/gtl/boost/polygon/voronoi_diagram.hpp
==============================================================================
--- sandbox/gtl/boost/polygon/voronoi_diagram.hpp (original)
+++ sandbox/gtl/boost/polygon/voronoi_diagram.hpp 2011-11-08 09:16:57 EST (Tue, 08 Nov 2011)
@@ -16,7 +16,7 @@
#include <vector>
#include "polygon.hpp"
-#include "detail/voronoi_fpt_kernel.hpp"
+#include "detail/voronoi_robust_fpt.hpp"
namespace boost {
namespace polygon {
Modified: sandbox/gtl/libs/polygon/test/Jamfile.v2
==============================================================================
--- sandbox/gtl/libs/polygon/test/Jamfile.v2 (original)
+++ sandbox/gtl/libs/polygon/test/Jamfile.v2 2011-11-08 09:16:57 EST (Tue, 08 Nov 2011)
@@ -29,9 +29,9 @@
alias "voronoi-unit"
:
[ run voronoi_builder_test.cpp ]
- [ run voronoi_calc_kernel_test.cpp ]
+ [ run voronoi_calc_utils_test.cpp ]
[ run voronoi_clipping_test.cpp ]
- [ run voronoi_fpt_kernel_test.cpp ]
+ [ run voronoi_robust_fpt_test.cpp ]
[ run voronoi_structures_test.cpp ]
;
Deleted: sandbox/gtl/libs/polygon/test/voronoi_calc_kernel_test.cpp
==============================================================================
--- sandbox/gtl/libs/polygon/test/voronoi_calc_kernel_test.cpp 2011-11-08 09:16:57 EST (Tue, 08 Nov 2011)
+++ (empty file)
@@ -1,456 +0,0 @@
-// Boost.Polygon library voronoi_calc_kernel_test.cpp file
-
-// Copyright Andrii Sydorchuk 2010-2011.
-// 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 <map>
-
-#define BOOST_TEST_MODULE voronoi_calc_kernel_test
-#include <boost/test/test_case_template.hpp>
-
-#include "boost/polygon/detail/voronoi_structures.hpp"
-#include "boost/polygon/detail/voronoi_calc_kernel.hpp"
-using namespace boost::polygon::detail;
-
-typedef voronoi_calc_kernel<int> VCK;
-typedef point_2d<int> point_type;
-typedef site_event<int> site_type;
-typedef circle_event<double> circle_type;
-VCK::event_comparison_predicate<site_type, circle_type> event_comparison;
-
-typedef beach_line_node_key<site_type> key_type;
-typedef VCK::distance_predicate<site_type> distance_predicate_type;
-typedef VCK::node_comparison_predicate<key_type> node_comparison_type;
-typedef std::map<key_type, int, node_comparison_type> beach_line_type;
-typedef beach_line_type::iterator bieach_line_iterator;
-distance_predicate_type distance_predicate;
-node_comparison_type node_comparison;
-
-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); \
- BOOST_CHECK_EQUAL(VCK::get_orientation(P1, P3, P2) == R2, true); \
- BOOST_CHECK_EQUAL(VCK::get_orientation(P2, P1, P3) == R2, true); \
- BOOST_CHECK_EQUAL(VCK::get_orientation(P2, P3, P1) == R1, true); \
- BOOST_CHECK_EQUAL(VCK::get_orientation(P3, P1, P2) == R1, true); \
- BOOST_CHECK_EQUAL(VCK::get_orientation(P3, P2, P1) == R2, true)
-
-#define CHECK_EVENT_COMPARISON(A, B, R1, R2) \
- BOOST_CHECK_EQUAL(event_comparison(A, B), R1); \
- BOOST_CHECK_EQUAL(event_comparison(B, A), R2)
-
-#define CHECK_DISTANCE_PREDICATE(S1, S2, S3, RES) \
- BOOST_CHECK_EQUAL(distance_predicate(S1, S2, S3), RES)
-
-#define CHECK_NODE_COMPARISON(node, nodes, res, sz) \
- for (int i = 0; i < sz; ++i) { \
- BOOST_CHECK_EQUAL(node_comparison(node, nodes[i]), res[i]); \
- 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();
- point_type point1(min_int, min_int);
- point_type point2(0, 0);
- point_type point3(max_int, max_int);
- point_type point4(min_int, max_int);
- point_type point5(max_int-1, max_int);
- CHECK_ORIENTATION(point1, point2, point3, VCK::COLLINEAR, VCK::COLLINEAR);
- CHECK_ORIENTATION(point1, point4, point3, VCK::RIGHT, VCK::LEFT);
- CHECK_ORIENTATION(point1, point5, point3, VCK::RIGHT, VCK::LEFT);
-}
-
-BOOST_AUTO_TEST_CASE(event_comparison_test1) {
- site_type site(1, 2);
- CHECK_EVENT_COMPARISON(site, site_type(0, 2), false, true);
- CHECK_EVENT_COMPARISON(site, site_type(1, 3), true, false);
- CHECK_EVENT_COMPARISON(site, site_type(1, 2), false, false);
-}
-
-BOOST_AUTO_TEST_CASE(event_comparison_test2) {
- site_type site(0, 0, 0, 2);
- CHECK_EVENT_COMPARISON(site, site_type(0, 2), true, false);
- CHECK_EVENT_COMPARISON(site, site_type(0, 0), false, true);
- CHECK_EVENT_COMPARISON(site, site_type(0, -2, 0, -1), false, true);
- CHECK_EVENT_COMPARISON(site, site_type(0, -2, 1, 1), true, false);
- CHECK_EVENT_COMPARISON(site, site_type(0, 0, 1, 1), true, false);
-}
-
-BOOST_AUTO_TEST_CASE(event_comparison_test3) {
- site_type site(0, 0, 10, 10);
- CHECK_EVENT_COMPARISON(site, site_type(0, 0), false, true);
- CHECK_EVENT_COMPARISON(site, site_type(0, -1), false, true);
- CHECK_EVENT_COMPARISON(site, site_type(0, 1), false, true);
- CHECK_EVENT_COMPARISON(site, site_type(0, 1, 0, 10), false, true);
- CHECK_EVENT_COMPARISON(site, site_type(0, -10, 0, -1), false, true);
- CHECK_EVENT_COMPARISON(site, site_type(0, 0, 10, 9), true, false);
- CHECK_EVENT_COMPARISON(site, site_type(0, 0, 9, 10), false, true);
-}
-
-BOOST_AUTO_TEST_CASE(event_comparison_test4) {
- circle_type circle(1, 2, 3);
- CHECK_EVENT_COMPARISON(circle, circle_type(1, 2, 3), false, false);
- CHECK_EVENT_COMPARISON(circle, circle_type(1, 3, 3), true, false);
- CHECK_EVENT_COMPARISON(circle, circle_type(1, 2, 4), true, false);
- CHECK_EVENT_COMPARISON(circle, circle_type(0, 2, 2), false, true);
- CHECK_EVENT_COMPARISON(circle, circle_type(-1, 2, 3), false, false);
-}
-
-BOOST_AUTO_TEST_CASE(event_comparison_test5) {
- circle_type circle(1, 2, 3);
- CHECK_EVENT_COMPARISON(circle, site_type(0, 100), false, true);
- CHECK_EVENT_COMPARISON(circle, site_type(3, 0), false, true);
- CHECK_EVENT_COMPARISON(circle, site_type(3, 2), false, false);
- CHECK_EVENT_COMPARISON(circle, site_type(3, 3), true, false);
- CHECK_EVENT_COMPARISON(circle, site_type(4, 2), true, false);
-}
-
-BOOST_AUTO_TEST_CASE(distance_predicate_test1) {
- site_type site1(-5, 0);
- site_type site2(-8, 9);
- site_type site3(-2, 1);
- CHECK_DISTANCE_PREDICATE(site1, site2, site_type(0, 5), false);
- CHECK_DISTANCE_PREDICATE(site3, site1, site_type(0, 5), false);
- CHECK_DISTANCE_PREDICATE(site1, site2, site_type(0, 4), false);
- CHECK_DISTANCE_PREDICATE(site3, site1, site_type(0, 4), false);
- CHECK_DISTANCE_PREDICATE(site1, site2, site_type(0, 6), true);
- CHECK_DISTANCE_PREDICATE(site3, site1, site_type(0, 6), true);
-}
-
-BOOST_AUTO_TEST_CASE(distance_predicate_test2) {
- site_type site1(-4, 0, -4, 20);
- site_type site2(-2, 10);
- CHECK_DISTANCE_PREDICATE(site2, site1, site_type(0, 11), false);
- CHECK_DISTANCE_PREDICATE(site2, site1, site_type(0, 9), false);
- CHECK_DISTANCE_PREDICATE(site1, site2, site_type(0, 11), true);
- CHECK_DISTANCE_PREDICATE(site1, site2, site_type(0, 9), true);
-}
-
-BOOST_AUTO_TEST_CASE(disntace_predicate_test3) {
- site_type site1(-5, 5, 2, -2);
- site1.inverse();
- site_type site2(-2, 4);
- CHECK_DISTANCE_PREDICATE(site1, site2, site_type(0, -1), false);
- CHECK_DISTANCE_PREDICATE(site2, site1, site_type(0, -1), false);
- CHECK_DISTANCE_PREDICATE(site1, site2, site_type(0, 1), false);
- CHECK_DISTANCE_PREDICATE(site2, site1, site_type(0, 1), false);
- CHECK_DISTANCE_PREDICATE(site1, site2, site_type(0, 4), true);
- CHECK_DISTANCE_PREDICATE(site2, site1, site_type(0, 4), false);
- CHECK_DISTANCE_PREDICATE(site1, site2, site_type(0, 5), true);
- CHECK_DISTANCE_PREDICATE(site2, site1, site_type(0, 5), false);
-}
-
-BOOST_AUTO_TEST_CASE(distance_predicate_test4) {
- site_type site1(-5, 5, 2, -2);
- site_type site2(-2, -4);
- site_type site3(-4, 1);
- CHECK_DISTANCE_PREDICATE(site1, site2, site_type(0, 1), true);
- CHECK_DISTANCE_PREDICATE(site2, site1, site_type(0, 1), true);
- CHECK_DISTANCE_PREDICATE(site1, site3, site_type(0, 1), true);
- CHECK_DISTANCE_PREDICATE(site3, site1, site_type(0, 1), true);
- CHECK_DISTANCE_PREDICATE(site1, site2, site_type(0, -2), true);
- CHECK_DISTANCE_PREDICATE(site2, site1, site_type(0, -2), false);
- CHECK_DISTANCE_PREDICATE(site1, site3, site_type(0, -2), true);
- CHECK_DISTANCE_PREDICATE(site3, site1, site_type(0, -2), false);
- CHECK_DISTANCE_PREDICATE(site1, site2, site_type(0, -8), true);
- CHECK_DISTANCE_PREDICATE(site2, site1, site_type(0, -8), false);
- CHECK_DISTANCE_PREDICATE(site1, site3, site_type(0, -8), true);
- CHECK_DISTANCE_PREDICATE(site3, site1, site_type(0, -8), false);
- CHECK_DISTANCE_PREDICATE(site1, site2, site_type(0, -9), true);
- CHECK_DISTANCE_PREDICATE(site2, site1, site_type(0, -9), false);
- CHECK_DISTANCE_PREDICATE(site1, site3, site_type(0, -9), true);
- CHECK_DISTANCE_PREDICATE(site3, site1, site_type(0, -9), false);
-}
-
-BOOST_AUTO_TEST_CASE(disntace_predicate_test5) {
- site_type site1(-5, 5, 2, -2);
- site_type site2 = site1;
- site2.inverse();
- site_type site3(-2, 4);
- site_type site4(-2, -4);
- site_type site5(-4, 1);
- CHECK_DISTANCE_PREDICATE(site3, site2, site_type(0, 1), false);
- CHECK_DISTANCE_PREDICATE(site3, site2, site_type(0, 4), false);
- CHECK_DISTANCE_PREDICATE(site3, site2, site_type(0, 5), false);
- CHECK_DISTANCE_PREDICATE(site3, site2, site_type(0, 7), true);
- CHECK_DISTANCE_PREDICATE(site4, site1, site_type(0, -2), false);
- CHECK_DISTANCE_PREDICATE(site5, site1, site_type(0, -2), false);
- CHECK_DISTANCE_PREDICATE(site4, site1, site_type(0, -8), false);
- CHECK_DISTANCE_PREDICATE(site5, site1, site_type(0, -8), false);
- CHECK_DISTANCE_PREDICATE(site4, site1, site_type(0, -9), false);
- CHECK_DISTANCE_PREDICATE(site5, site1, site_type(0, -9), false);
- CHECK_DISTANCE_PREDICATE(site4, site1, site_type(0, -18), false);
- CHECK_DISTANCE_PREDICATE(site5, site1, site_type(0, -18), false);
- CHECK_DISTANCE_PREDICATE(site4, site1, site_type(0, -1), true);
- CHECK_DISTANCE_PREDICATE(site5, site1, site_type(0, -1), true);
-}
-
-BOOST_AUTO_TEST_CASE(distance_predicate_test6) {
- site_type site1(-5, 0, 2, 7);
- site_type site2 = site1;
- site2.inverse();
- CHECK_DISTANCE_PREDICATE(site1, site2, site_type(2, 7), false);
- CHECK_DISTANCE_PREDICATE(site1, site2, site_type(1, 5), false);
- CHECK_DISTANCE_PREDICATE(site1, site2, site_type(-1, 5), true);
-}
-
-BOOST_AUTO_TEST_CASE(distance_predicate_test7) {
- site_type site1(-5, 5, 2, -2);
- site1.inverse();
- site_type site2(-5, 5, 0, 6);
- site_type site3(-2, 4, 0, 4);
- site_type site4(0, 2);
- site_type site5(0, 5);
- site_type site6(0, 6);
- site_type site7(0, 8);
- CHECK_DISTANCE_PREDICATE(site1, site2, site4, false);
- CHECK_DISTANCE_PREDICATE(site1, site2, site5, true);
- CHECK_DISTANCE_PREDICATE(site1, site2, site6, true);
- CHECK_DISTANCE_PREDICATE(site1, site2, site7, true);
- CHECK_DISTANCE_PREDICATE(site1, site3, site4, false);
- CHECK_DISTANCE_PREDICATE(site1, site3, site5, true);
- CHECK_DISTANCE_PREDICATE(site1, site3, site6, true);
- CHECK_DISTANCE_PREDICATE(site1, site3, site7, true);
- site3.inverse();
- CHECK_DISTANCE_PREDICATE(site3, site1, site4, false);
- CHECK_DISTANCE_PREDICATE(site3, site1, site5, false);
- CHECK_DISTANCE_PREDICATE(site3, site1, site6, false);
- CHECK_DISTANCE_PREDICATE(site3, site1, site7, true);
-}
-
-BOOST_AUTO_TEST_CASE(distatnce_predicate_test8) {
- site_type site1(-5, 3, -2, 2);
- site1.inverse();
- site_type site2(-5, 5, -2, 2);
- CHECK_DISTANCE_PREDICATE(site1, site2, site_type(-4, 2), false);
-}
-
-BOOST_AUTO_TEST_CASE(node_comparison_test1) {
- beach_line_type beach_line;
- site_type site1(0, 0);
- site1.index(0);
- site_type site2(0, 2);
- site2.index(1);
- site_type site3(1, 0);
- site3.index(2);
- beach_line[key_type(site1, site2)] = 2;
- beach_line[key_type(site1, site3)] = 0;
- beach_line[key_type(site3, site1)] = 1;
- int cur_index = 0;
- for (bieach_line_iterator it = beach_line.begin();
- it != beach_line.end(); ++it, ++cur_index) {
- BOOST_CHECK_EQUAL(it->second, cur_index);
- }
-}
-
-BOOST_AUTO_TEST_CASE(node_comparison_test2) {
- beach_line_type beach_line;
- site_type site1(0, 1);
- site1.index(0);
- site_type site2(2, 0);
- site2.index(1);
- site_type site3(2, 4);
- site3.index(2);
- beach_line[key_type(site1, site2)] = 0;
- beach_line[key_type(site2, site1)] = 1;
- beach_line[key_type(site1, site3)] = 2;
- beach_line[key_type(site3, site1)] = 3;
- int cur_index = 0;
- for (bieach_line_iterator it = beach_line.begin();
- it != beach_line.end(); ++it, ++cur_index) {
- BOOST_CHECK_EQUAL(it->second, cur_index);
- }
-}
-
-BOOST_AUTO_TEST_CASE(node_comparison_test3) {
- key_type node(site_type(1, 0).index(1), site_type(0, 2).index(0));
- key_type nodes[] = {
- key_type(site_type(2, -10).index(2)),
- key_type(site_type(2, -1).index(2)),
- key_type(site_type(2, 0).index(2)),
- key_type(site_type(2, 1).index(2)),
- key_type(site_type(2, 2).index(2)),
- key_type(site_type(2, 3).index(2)),
- };
- bool res[] = {false, false, false, false, true, true};
- CHECK_NODE_COMPARISON(node, nodes, res, 6);
-}
-
-BOOST_AUTO_TEST_CASE(node_comparison_test4) {
- key_type node(site_type(0, 1).index(0), site_type(1, 0).index(1));
- key_type nodes[] = {
- key_type(site_type(2, -3).index(2)),
- key_type(site_type(2, -2).index(2)),
- key_type(site_type(2, -1).index(2)),
- key_type(site_type(2, 0).index(2)),
- key_type(site_type(2, 1).index(2)),
- key_type(site_type(2, 3).index(2)),
- };
- bool res[] = {false, true, true, true, true, true};
- CHECK_NODE_COMPARISON(node, nodes, res, 6);
-}
-
-BOOST_AUTO_TEST_CASE(node_comparison_test5) {
- key_type node(site_type(0, 0).index(0), site_type(1, 2).index(1));
- key_type nodes[] = {
- key_type(site_type(2, -10).index(2)),
- key_type(site_type(2, 0).index(2)),
- key_type(site_type(2, 1).index(2)),
- key_type(site_type(2, 2).index(2)),
- key_type(site_type(2, 5).index(2)),
- key_type(site_type(2, 20).index(2)),
- };
- bool res[] = {false, false, true, true, true, true};
- CHECK_NODE_COMPARISON(node, nodes, res, 6);
-}
-
-BOOST_AUTO_TEST_CASE(node_comparison_test6) {
- key_type node(site_type(1, 1).index(1), site_type(0, 0).index(0));
- key_type nodes [] = {
- key_type(site_type(2, -3).index(2)),
- key_type(site_type(2, -2).index(2)),
- key_type(site_type(2, 0).index(2)),
- key_type(site_type(2, 1).index(2)),
- key_type(site_type(2, 2).index(2)),
- key_type(site_type(2, 3).index(2)),
- key_type(site_type(2, 5).index(2)),
- };
- bool res[] = {false, false, false, false, false, false, true};
- CHECK_NODE_COMPARISON(node, nodes, res, 7);
-}
-
-BOOST_AUTO_TEST_CASE(node_comparison_test7) {
- key_type node(site_type(0, 0).index(0), site_type(0, 2).index(1));
- key_type nodes[] = {
- key_type(site_type(1, 0).index(2)),
- key_type(site_type(1, 1).index(2)),
- key_type(site_type(1, 2).index(2)),
- };
- bool res[] = {false, false, true};
- CHECK_NODE_COMPARISON(node, nodes, res, 3);
-}
-
-BOOST_AUTO_TEST_CASE(node_comparison_test8) {
- key_type node(site_type(0, 0).index(0), site_type(1, 1).index(2));
- key_type nodes[] = {
- key_type(site_type(1, 0).index(1)),
- key_type(site_type(1, 1).index(2)),
- key_type(site_type(1, 2).index(3)),
- key_type(site_type(1, 1).index(2), site_type(0, 0).index(0)),
- };
- 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);
- site_type site2(-8, 0);
- site_type site3(0, 6);
- 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);
- site_type site2(min_int, max_int);
- site_type site3(max_int-1, max_int-1);
- site_type site4(max_int, max_int);
- 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);
- site_type site2(0, 4);
- site_type site3(site1.point0(), site2.point0());
- CHECK_CIRCLE_EXISTENCE(site1, site3, site2, false);
- site_type site4(-2, 0);
- site_type site5(0, 2);
- 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);
- site_type site2(-2, 10);
- site_type site3(4, 10);
- 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);
- site1.inverse();
- site_type site2(-2, 4, 10, 4);
- site_type site3(6, 2);
- site_type site4(1, 0);
- 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);
- site1.inverse();
- site_type site2(-1, 0, 8, 12);
- site_type site3(1, 1);
- 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);
- site1.inverse();
- site_type site2(-6, 4, 0, 12);
- site_type site3(1, 0);
- 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);
- site1.inverse();
- site_type site2(0, 12, 8, 6);
- site_type site3(1, 0);
- CHECK_CIRCLE_FORMATION_PREDICATE(site3, site2, site1, 1.0, 5.0, 6.0);
-}
-
-BOOST_AUTO_TEST_CASE(circle_formation_predicate_test10) {
- site_type site1(0, 0, 4, 0);
- site_type site2(0, 0, 0, 4);
- site_type site3(0, 4, 4, 4);
- 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(1, 0, 41, 30);
- site_type site2(-39, 30, 1, 60);
- site_type site3(1, 60, 41, 30);
- site1.inverse();
- CHECK_CIRCLE_FORMATION_PREDICATE(site1, site2, site3, 1.0, 30.0, 25.0);
-}
Copied: sandbox/gtl/libs/polygon/test/voronoi_calc_utils_test.cpp (from r75182, /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_utils_test.cpp 2011-11-08 09:16:57 EST (Tue, 08 Nov 2011)
@@ -1,4 +1,4 @@
-// Boost.Polygon library voronoi_calc_kernel_test.cpp file
+// Boost.Polygon library voronoi_calc_utils_test.cpp file
// Copyright Andrii Sydorchuk 2010-2011.
// Distributed under the Boost Software License, Version 1.0.
@@ -9,40 +9,40 @@
#include <map>
-#define BOOST_TEST_MODULE voronoi_calc_kernel_test
+#define BOOST_TEST_MODULE voronoi_calc_utils_test
#include <boost/test/test_case_template.hpp>
#include "boost/polygon/detail/voronoi_structures.hpp"
-#include "boost/polygon/detail/voronoi_calc_kernel.hpp"
+#include "boost/polygon/detail/voronoi_calc_utils.hpp"
using namespace boost::polygon::detail;
-typedef voronoi_calc_kernel<int> VCK;
+typedef voronoi_calc_utils<int> VCU;
typedef point_2d<int> point_type;
typedef site_event<int> site_type;
typedef circle_event<double> circle_type;
-VCK::event_comparison_predicate<site_type, circle_type> event_comparison;
+VCU::event_comparison_predicate<site_type, circle_type> event_comparison;
typedef beach_line_node_key<site_type> key_type;
-typedef VCK::distance_predicate<site_type> distance_predicate_type;
-typedef VCK::node_comparison_predicate<key_type> node_comparison_type;
+typedef VCU::distance_predicate<site_type> distance_predicate_type;
+typedef VCU::node_comparison_predicate<key_type> node_comparison_type;
typedef std::map<key_type, int, node_comparison_type> beach_line_type;
typedef beach_line_type::iterator bieach_line_iterator;
distance_predicate_type distance_predicate;
node_comparison_type node_comparison;
-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;
+typedef VCU::circle_existence_predicate<site_type> CEP_type;
+typedef VCU::gmp_circle_formation_functor<site_type, circle_type> GMP_CFF_type;
+typedef VCU::lazy_circle_formation_functor<site_type, circle_type> lazy_CFF_type;
+VCU::circle_formation_predicate<site_type, circle_type, CEP_type, GMP_CFF_type> gmp_predicate;
+VCU::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); \
- BOOST_CHECK_EQUAL(VCK::get_orientation(P1, P3, P2) == R2, true); \
- BOOST_CHECK_EQUAL(VCK::get_orientation(P2, P1, P3) == R2, true); \
- BOOST_CHECK_EQUAL(VCK::get_orientation(P2, P3, P1) == R1, true); \
- BOOST_CHECK_EQUAL(VCK::get_orientation(P3, P1, P2) == R1, true); \
- BOOST_CHECK_EQUAL(VCK::get_orientation(P3, P2, P1) == R2, true)
+ BOOST_CHECK_EQUAL(VCU::get_orientation(P1, P2, P3) == R1, true); \
+ BOOST_CHECK_EQUAL(VCU::get_orientation(P1, P3, P2) == R2, true); \
+ BOOST_CHECK_EQUAL(VCU::get_orientation(P2, P1, P3) == R2, true); \
+ BOOST_CHECK_EQUAL(VCU::get_orientation(P2, P3, P1) == R1, true); \
+ BOOST_CHECK_EQUAL(VCU::get_orientation(P3, P1, P2) == R1, true); \
+ BOOST_CHECK_EQUAL(VCU::get_orientation(P3, P2, P1) == R2, true)
#define CHECK_EVENT_COMPARISON(A, B, R1, R2) \
BOOST_CHECK_EQUAL(event_comparison(A, B), R1); \
@@ -81,9 +81,9 @@
point_type point3(max_int, max_int);
point_type point4(min_int, max_int);
point_type point5(max_int-1, max_int);
- CHECK_ORIENTATION(point1, point2, point3, VCK::COLLINEAR, VCK::COLLINEAR);
- CHECK_ORIENTATION(point1, point4, point3, VCK::RIGHT, VCK::LEFT);
- CHECK_ORIENTATION(point1, point5, point3, VCK::RIGHT, VCK::LEFT);
+ CHECK_ORIENTATION(point1, point2, point3, VCU::COLLINEAR, VCU::COLLINEAR);
+ CHECK_ORIENTATION(point1, point4, point3, VCU::RIGHT, VCU::LEFT);
+ CHECK_ORIENTATION(point1, point5, point3, VCU::RIGHT, VCU::LEFT);
}
BOOST_AUTO_TEST_CASE(event_comparison_test1) {
Deleted: sandbox/gtl/libs/polygon/test/voronoi_fpt_kernel_test.cpp
==============================================================================
--- sandbox/gtl/libs/polygon/test/voronoi_fpt_kernel_test.cpp 2011-11-08 09:16:57 EST (Tue, 08 Nov 2011)
+++ (empty file)
@@ -1,276 +0,0 @@
-// Boost.Polygon library voronoi_fpt_kernel_test.cpp file
-
-// Copyright Andrii Sydorchuk 2010-2011.
-// 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 <cmath>
-#include <ctime>
-#include <vector>
-
-#define BOOST_TEST_MODULE voronoi_fpt_kernel_test
-#include <boost/mpl/list.hpp>
-#include <boost/random/mersenne_twister.hpp>
-#include <boost/test/test_case_template.hpp>
-
-#include "boost/polygon/detail/voronoi_fpt_kernel.hpp"
-#include "boost/polygon/detail/voronoi_calc_kernel.hpp"
-using namespace boost::polygon::detail;
-
-typedef robust_fpt<double> rfpt_type;
-
-BOOST_AUTO_TEST_CASE(robust_fpt_constructors_test1) {
- rfpt_type a = rfpt_type();
- BOOST_CHECK_EQUAL(a.fpv() == 0.0, true);
- BOOST_CHECK_EQUAL(a.re() == 0.0, true);
- BOOST_CHECK_EQUAL(a.ulp() == 0, true);
-}
-
-BOOST_AUTO_TEST_CASE(robust_fpt_constructors_test2) {
- rfpt_type a(10.0);
- BOOST_CHECK_EQUAL(a.fpv() == 10.0, true);
- BOOST_CHECK_EQUAL(a.re() == 1.0, true);
- BOOST_CHECK_EQUAL(a.ulp() == 1.0, true);
-}
-
-BOOST_AUTO_TEST_CASE(robust_fpt_constructors_test3) {
- rfpt_type a(10.0, false);
- BOOST_CHECK_EQUAL(a.fpv() == 10.0, true);
- BOOST_CHECK_EQUAL(a.re() == 0.0, true);
- BOOST_CHECK_EQUAL(a.ulp() == 0.0, true);
-}
-
-BOOST_AUTO_TEST_CASE(robust_fpt_constructors_test4) {
- rfpt_type a(10.0, 3.0);
- BOOST_CHECK_EQUAL(a.fpv() == 10.0, true);
- BOOST_CHECK_EQUAL(a.re() == 3.0, true);
- BOOST_CHECK_EQUAL(a.ulp() == 3.0, true);
-
- rfpt_type b(10.0, 2.75);
- BOOST_CHECK_EQUAL(b.fpv() == 10.0, true);
- BOOST_CHECK_EQUAL(b.re() == 2.75, true);
- BOOST_CHECK_EQUAL(b.ulp() == 2.75, true);
-}
-
-BOOST_AUTO_TEST_CASE(robust_fpt_sum_test1) {
- rfpt_type a(2.0, 5.0);
- rfpt_type b(3.0, 4.0);
- rfpt_type c = a + b;
- BOOST_CHECK_EQUAL(c.fpv() == 5.0, true);
- BOOST_CHECK_EQUAL(c.re() == 6.0, true);
- BOOST_CHECK_EQUAL(c.ulp() == 6.0, true);
-
- c += b;
- BOOST_CHECK_EQUAL(c.fpv() == 8.0, true);
- BOOST_CHECK_EQUAL(c.re() == 7.0, true);
- BOOST_CHECK_EQUAL(c.ulp() == 7.0, true);
-}
-
-BOOST_AUTO_TEST_CASE(robust_fpt_sum_test2) {
- rfpt_type a(3.0, 2.0);
- rfpt_type b(-2.0, 3.0);
- rfpt_type c = a + b;
- BOOST_CHECK_EQUAL(c.fpv() == 1.0, true);
- BOOST_CHECK_EQUAL(c.re() == 13.0, true);
- BOOST_CHECK_EQUAL(c.ulp() == 13.0, true);
-
- c += b;
- BOOST_CHECK_EQUAL(c.fpv() == -1.0, true);
- BOOST_CHECK_EQUAL(c.re() == 20.0, true);
- BOOST_CHECK_EQUAL(c.ulp() == 20.0, true);
-}
-
-BOOST_AUTO_TEST_CASE(robust_fpt_dif_test1) {
- rfpt_type a(2.0, 5.0);
- rfpt_type b(-3.0, 4.0);
- rfpt_type c = a - b;
- BOOST_CHECK_EQUAL(c.fpv() == 5.0, true);
- BOOST_CHECK_EQUAL(c.re() == 6.0, true);
- BOOST_CHECK_EQUAL(c.ulp() == 6.0, true);
-
- c -= b;
- BOOST_CHECK_EQUAL(c.fpv() == 8.0, true);
- BOOST_CHECK_EQUAL(c.re() == 7.0, true);
- BOOST_CHECK_EQUAL(c.ulp() == 7.0, true);
-}
-
-BOOST_AUTO_TEST_CASE(robust_fpt_dif_test2) {
- rfpt_type a(3.0, 2.0);
- rfpt_type b(2.0, 3.0);
- rfpt_type c = a - b;
- BOOST_CHECK_EQUAL(c.fpv() == 1.0, true);
- BOOST_CHECK_EQUAL(c.re() == 13.0, true);
- BOOST_CHECK_EQUAL(c.ulp() == 13.0, true);
-
- c -= b;
- BOOST_CHECK_EQUAL(c.fpv() == -1.0, true);
- BOOST_CHECK_EQUAL(c.re() == 20.0, true);
- BOOST_CHECK_EQUAL(c.ulp() == 20.0, true);
-}
-
-BOOST_AUTO_TEST_CASE(robust_fpt_mult_test3) {
- rfpt_type a(2.0, 3.0);
- rfpt_type b(4.0);
- rfpt_type c = a * b;
- BOOST_CHECK_EQUAL(c.fpv() == 8.0, true);
- BOOST_CHECK_EQUAL(c.re() == 5.0, true);
- BOOST_CHECK_EQUAL(c.ulp() == 5.0, true);
-
- c *= b;
- BOOST_CHECK_EQUAL(c.fpv() == 32.0, true);
- BOOST_CHECK_EQUAL(c.re() == 7.0, true);
- BOOST_CHECK_EQUAL(c.ulp() == 7.0, true);
-}
-
-BOOST_AUTO_TEST_CASE(robust_fpt_div_test1) {
- rfpt_type a(2.0, 3.0);
- rfpt_type b(4.0);
- rfpt_type c = a / b;
- BOOST_CHECK_EQUAL(c.fpv() == 0.5, true);
- BOOST_CHECK_EQUAL(c.re() == 5.0, true);
- BOOST_CHECK_EQUAL(c.ulp() == 5.0, true);
-
- c /= b;
- BOOST_CHECK_EQUAL(c.fpv() == 0.125, true);
- BOOST_CHECK_EQUAL(c.re() == 7.0, true);
- BOOST_CHECK_EQUAL(c.ulp() == 7.0, true);
-}
-
-BOOST_AUTO_TEST_CASE(robust_dif_constructors_test) {
- robust_dif<int> rd1;
- BOOST_CHECK_EQUAL(rd1.pos(), 0);
- BOOST_CHECK_EQUAL(rd1.neg(), 0);
- BOOST_CHECK_EQUAL(rd1.dif(), 0);
-
- robust_dif<int> rd2(1);
- BOOST_CHECK_EQUAL(rd2.pos(), 1);
- BOOST_CHECK_EQUAL(rd2.neg(), 0);
- BOOST_CHECK_EQUAL(rd2.dif(), 1);
-
- robust_dif<int> rd3(-1);
- BOOST_CHECK_EQUAL(rd3.pos(), 0);
- BOOST_CHECK_EQUAL(rd3.neg(), 1);
- BOOST_CHECK_EQUAL(rd3.dif(), -1);
-
- robust_dif<int> rd4(1, 2);
- BOOST_CHECK_EQUAL(rd4.pos(), 1);
- BOOST_CHECK_EQUAL(rd4.neg(), 2);
- BOOST_CHECK_EQUAL(rd4.dif(), -1);
-}
-
-BOOST_AUTO_TEST_CASE(robust_dif_operators_test1) {
- robust_dif<int> a(5, 2), b(1, 10);
- int dif_a = a.dif();
- int dif_b = b.dif();
- robust_dif<int> sum = a + b;
- robust_dif<int> dif = a - b;
- robust_dif<int> mult = a * b;
- robust_dif<int> umin = -a;
- BOOST_CHECK_EQUAL(sum.dif(), dif_a + dif_b);
- BOOST_CHECK_EQUAL(dif.dif(), dif_a - dif_b);
- BOOST_CHECK_EQUAL(mult.dif(), dif_a * dif_b);
- BOOST_CHECK_EQUAL(umin.dif(), -dif_a);
-}
-
-BOOST_AUTO_TEST_CASE(robust_dif_operators_test2) {
- robust_dif<int> a(5, 2);
- for (int b = -3; b <= 3; b += 6) {
- int dif_a = a.dif();
- int dif_b = b;
- robust_dif<int> sum = a + b;
- robust_dif<int> dif = a - b;
- robust_dif<int> mult = a * b;
- robust_dif<int> div = a / b;
- BOOST_CHECK_EQUAL(sum.dif(), dif_a + dif_b);
- BOOST_CHECK_EQUAL(dif.dif(), dif_a - dif_b);
- BOOST_CHECK_EQUAL(mult.dif(), dif_a * dif_b);
- BOOST_CHECK_EQUAL(div.dif(), dif_a / dif_b);
- }
-}
-
-BOOST_AUTO_TEST_CASE(robust_dif_operators_test3) {
- robust_dif<int> b(5, 2);
- for (int a = -3; a <= 3; a += 6) {
- int dif_a = a;
- int dif_b = b.dif();
- robust_dif<int> sum = a + b;
- robust_dif<int> dif = a - b;
- robust_dif<int> mult = a * b;
- BOOST_CHECK_EQUAL(sum.dif(), dif_a + dif_b);
- BOOST_CHECK_EQUAL(dif.dif(), dif_a - dif_b);
- BOOST_CHECK_EQUAL(mult.dif(), dif_a * dif_b);
- }
-}
-
-BOOST_AUTO_TEST_CASE(robust_dif_operators_test4) {
- std::vector< robust_dif<int> > a4(4, robust_dif<int>(5, 2));
- std::vector< robust_dif<int> > b4(4, robust_dif<int>(1, 2));
- std::vector< robust_dif<int> > c4 = a4;
- c4[0] += b4[0];
- c4[1] -= b4[1];
- c4[2] *= b4[2];
- BOOST_CHECK_EQUAL(c4[0].dif(), a4[0].dif() + b4[0].dif());
- BOOST_CHECK_EQUAL(c4[1].dif(), a4[1].dif() - b4[1].dif());
- BOOST_CHECK_EQUAL(c4[2].dif(), a4[2].dif() * b4[2].dif());
- a4[0] += b4[0].dif();
- a4[1] -= b4[1].dif();
- a4[2] *= b4[2].dif();
- a4[3] /= b4[3].dif();
- BOOST_CHECK_EQUAL(c4[0].dif(), a4[0].dif());
- BOOST_CHECK_EQUAL(c4[1].dif(), a4[1].dif());
- BOOST_CHECK_EQUAL(c4[2].dif(), a4[2].dif());
- BOOST_CHECK_EQUAL(c4[3].dif() / b4[3].dif(), a4[3].dif());
-}
-
-template <typename mpz, typename mpf>
-class sqrt_expr_tester {
-public:
- bool run() {
- static boost::mt19937 gen(static_cast<unsigned int>(time(NULL)));
- bool ret_val = true;
- for (int i = 0; i < MX_SQRTS; ++i) {
- a[i] = gen() % 1000000;
- double random_int = gen() % 1000000;
- b[i] = random_int * random_int;
- }
- int mask = (1 << MX_SQRTS);
- for (int i = 0; i < mask; i++) {
- double expected_val = 0.0;
- for (int j = 0; j < MX_SQRTS; 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 = get_d(sqrt_expr_.eval4(A, B));
- ret_val &= almost_equal(expected_val, received_val, 1);
- }
- return ret_val;
- }
-private:
- static const int MX_SQRTS = 4;
- robust_sqrt_expr<mpz, mpf> sqrt_expr_;
- mpz A[MX_SQRTS];
- mpz B[MX_SQRTS];
- double a[MX_SQRTS];
- double b[MX_SQRTS];
-};
-
-BOOST_AUTO_TEST_CASE(mpz_sqrt_evaluator_test) {
- typedef mpt_wrapper<mpz_class, 8> mpz_wrapper_type;
- typedef mpt_wrapper<mpf_class, 2> mpf_wrapper_type;
- sqrt_expr_tester<mpz_class, mpf_class> tester1;
- sqrt_expr_tester<mpz_wrapper_type, mpf_wrapper_type> tester2;
- for (int i = 0; i < 2000; ++i) {
- BOOST_CHECK(tester1.run());
- BOOST_CHECK(tester2.run());
- }
-}
Copied: sandbox/gtl/libs/polygon/test/voronoi_robust_fpt_test.cpp (from r75083, /sandbox/gtl/libs/polygon/test/voronoi_fpt_kernel_test.cpp)
==============================================================================
--- /sandbox/gtl/libs/polygon/test/voronoi_fpt_kernel_test.cpp (original)
+++ sandbox/gtl/libs/polygon/test/voronoi_robust_fpt_test.cpp 2011-11-08 09:16:57 EST (Tue, 08 Nov 2011)
@@ -1,4 +1,4 @@
-// Boost.Polygon library voronoi_fpt_kernel_test.cpp file
+// Boost.Polygon library voronoi_robust_fpt_test.cpp file
// Copyright Andrii Sydorchuk 2010-2011.
// Distributed under the Boost Software License, Version 1.0.
@@ -11,13 +11,13 @@
#include <ctime>
#include <vector>
-#define BOOST_TEST_MODULE voronoi_fpt_kernel_test
+#define BOOST_TEST_MODULE voronoi_robust_fpt_test
#include <boost/mpl/list.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <boost/test/test_case_template.hpp>
-#include "boost/polygon/detail/voronoi_fpt_kernel.hpp"
-#include "boost/polygon/detail/voronoi_calc_kernel.hpp"
+#include "boost/polygon/detail/voronoi_robust_fpt.hpp"
+#include "boost/polygon/detail/voronoi_calc_utils.hpp"
using namespace boost::polygon::detail;
typedef robust_fpt<double> rfpt_type;
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