|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r68306 - in sandbox/SOC/2010/sweepline: boost/sweepline boost/sweepline/detail libs/sweepline/example
From: sydorchuk.andriy_at_[hidden]
Date: 2011-01-19 14:03:47
Author: asydorchuk
Date: 2011-01-19 14:03:45 EST (Wed, 19 Jan 2011)
New Revision: 68306
URL: http://svn.boost.org/trac/boost/changeset/68306
Log:
Changed function that interpolates parabolic arcs with line segments.
Changed unsigned long long to BOOST_POLYGON_UNSIGNED_LONG_LONG.
Text files modified:
sandbox/SOC/2010/sweepline/boost/sweepline/detail/voronoi_formation.hpp | 54 ++++++------
sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_output.hpp | 166 ++++++++++++++++++++++-----------------
sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_sweepline.hpp | 3
sandbox/SOC/2010/sweepline/libs/sweepline/example/voronoi_visualizer.cpp | 2
4 files changed, 123 insertions(+), 102 deletions(-)
Modified: sandbox/SOC/2010/sweepline/boost/sweepline/detail/voronoi_formation.hpp
==============================================================================
--- sandbox/SOC/2010/sweepline/boost/sweepline/detail/voronoi_formation.hpp (original)
+++ sandbox/SOC/2010/sweepline/boost/sweepline/detail/voronoi_formation.hpp 2011-01-19 14:03:45 EST (Wed, 19 Jan 2011)
@@ -631,12 +631,13 @@
// Convert value to 64-bit unsigned integer.
// Return true if the value is positive, else false.
template <typename T>
- static inline bool convert_to_65_bit(T value, unsigned long long &res) {
+ static inline bool convert_to_65_bit(
+ T value, BOOST_POLYGON_UNSIGNED_LONG_LONG &res) {
if (value >= 0) {
- res = static_cast<unsigned long long>(value);
+ res = static_cast<BOOST_POLYGON_UNSIGNED_LONG_LONG>(value);
return true;
} else {
- res = static_cast<unsigned long long>(-value);
+ res = static_cast<BOOST_POLYGON_UNSIGNED_LONG_LONG>(-value);
return false;
}
}
@@ -645,8 +646,9 @@
// 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 reinterpretatoins differ in not more than maxUlps units.
- static inline bool almost_equal(double a, double b, int maxUlps) {
- long long ll_a, ll_b;
+ static inline bool almost_equal(double a, double b,
+ unsigned int maxUlps) {
+ BOOST_POLYGON_UNSIGNED_LONG_LONG ll_a, ll_b;
// Reinterpret double bits as 64-bit signed integer.
memcpy(&ll_a, &a, sizeof(double));
@@ -656,16 +658,17 @@
// 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 < 0)
- ll_a = 0x8000000000000000LL - ll_a;
- if (ll_b < 0)
- ll_b = 0x8000000000000000LL - ll_b;
+ 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.
- long long dif = ll_a - ll_b;
- return (dif <= maxUlps) && (dif >= -maxUlps);
+ if (ll_a > ll_b)
+ return ll_a - ll_b <= maxUlps;
+ return ll_b - ll_a <= maxUlps;
}
// Robust orientation test. Works correctly for any input type that
@@ -677,16 +680,15 @@
template <typename T>
static kOrientation orientation_test(T dif_x1_, T dif_y1_,
T dif_x2_, T dif_y2_) {
- typedef unsigned long long ull;
- ull dif_x1, dif_y1, dif_x2, dif_y2;
+ BOOST_POLYGON_UNSIGNED_LONG_LONG dif_x1, dif_y1, dif_x2, dif_y2;
bool dif_x1_plus, dif_x2_plus, dif_y1_plus, dif_y2_plus;
dif_x1_plus = convert_to_65_bit(dif_x1_, dif_x1);
dif_y1_plus = convert_to_65_bit(dif_y1_, dif_y1);
dif_x2_plus = convert_to_65_bit(dif_x2_, dif_x2);
dif_y2_plus = convert_to_65_bit(dif_y2_, dif_y2);
- ull expr_l = dif_x1 * dif_y2;
- ull expr_r = dif_x2 * dif_y1;
+ BOOST_POLYGON_UNSIGNED_LONG_LONG expr_l = dif_x1 * dif_y2;
+ BOOST_POLYGON_UNSIGNED_LONG_LONG expr_r = dif_x2 * dif_y1;
bool expr_l_plus = (dif_x1_plus == dif_y2_plus) ? true : false;
bool expr_r_plus = (dif_x2_plus == dif_y1_plus) ? true : false;
@@ -744,17 +746,16 @@
// The result is correct with epsilon relative error equal to 1EPS.
template <typename T>
static double robust_cross_product(T a1_, T b1_, T a2_, T b2_) {
- typedef unsigned long long ull;
- ull a1, b1, a2, b2;
+ BOOST_POLYGON_UNSIGNED_LONG_LONG 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);
- ull expr_l = a1 * b2;
+ BOOST_POLYGON_UNSIGNED_LONG_LONG expr_l = a1 * b2;
bool expr_l_plus = (a1_plus == b2_plus) ? true : false;
- ull expr_r = b1 * a2;
+ BOOST_POLYGON_UNSIGNED_LONG_LONG expr_r = b1 * a2;
bool expr_r_plus = (a2_plus == b1_plus) ? true : false;
if (expr_l == 0)
@@ -1062,9 +1063,8 @@
static bool robust_65bit_less_predicate(const point_2d<T> &left_point,
const point_2d<T> &right_point,
const point_2d<T> &new_point) {
- typedef long long ll;
- typedef unsigned long long ull;
- ull a1, a2, b1, b2, b1_sqr, b2_sqr, l_expr, r_expr;
+ BOOST_POLYGON_UNSIGNED_LONG_LONG a1, a2, b1, b2;
+ BOOST_POLYGON_UNSIGNED_LONG_LONG b1_sqr, b2_sqr, l_expr, r_expr;
bool l_expr_plus, r_expr_plus;
// Compute a1 = (x0-x1), a2 = (x0-x2), b1 = (y0 - y1), b2 = (y0 - y2).
@@ -1076,8 +1076,8 @@
// Compute b1_sqr = (y0-y1)^2 and b2_sqr = (y0-y2)^2.
b1_sqr = b1 * b1;
b2_sqr = b2 * b2;
- ull b1_sqr_mod = b1_sqr % a1;
- ull b2_sqr_mod = b2_sqr % a2;
+ BOOST_POLYGON_UNSIGNED_LONG_LONG b1_sqr_mod = b1_sqr % a1;
+ BOOST_POLYGON_UNSIGNED_LONG_LONG b2_sqr_mod = b2_sqr % a2;
// Compute l_expr = (x1 - x2).
l_expr_plus = convert_to_65_bit(left_point.x() - right_point.x(), l_expr);
@@ -1095,8 +1095,8 @@
}
// Compute right expression.
- ull value1 = b1_sqr / a1;
- ull value2 = b2_sqr / a2;
+ BOOST_POLYGON_UNSIGNED_LONG_LONG value1 = b1_sqr / a1;
+ BOOST_POLYGON_UNSIGNED_LONG_LONG value2 = b2_sqr / a2;
if (value1 >= value2) {
r_expr = value1 - value2;
r_expr_plus = true;
@@ -1157,8 +1157,6 @@
template <typename T>
static kPredicateResult fast_less_predicate(point_2d<T> site_point, site_event<T> segment_site,
point_2d<T> new_point, bool reverse_order) {
- typedef long long ll;
- typedef unsigned long long ull;
if (orientation_test(segment_site.point0(true), segment_site.point1(true),
new_point) != RIGHT_ORIENTATION) {
return (!segment_site.is_inverse()) ? LESS : MORE;
Modified: sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_output.hpp
==============================================================================
--- sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_output.hpp (original)
+++ sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_output.hpp 2011-01-19 14:03:45 EST (Wed, 19 Jan 2011)
@@ -172,6 +172,14 @@
return y_max_;
}
+ coordinate_type min_len() const {
+ return (std::min)(x_max_ - x_min_, y_max_ - y_min_);
+ }
+
+ coordinate_type max_len() const {
+ return (std::max)(x_max_ - x_min_, y_max_ - y_min_);
+ }
+
private:
coordinate_type x_min_;
coordinate_type y_min_;
@@ -218,9 +226,13 @@
// Compute intermediate points for the voronoi edge within the given
// bounding rectangle. The assumption is made that voronoi rectangle
// contains all the finite endpoints of the edge.
+ // Max_error is the maximum distance (relative to the minor of two
+ // rectangle sides) allowed between the parabola and line segments
+ // that interpolate it.
static std::vector<point_2d_type> get_intermediate_points(
const voronoi_edge<coordinate_type> *edge,
- const BRect<coordinate_type> &brect) {
+ const BRect<coordinate_type> &brect,
+ double max_error) {
// Retrieve the cell to the left of the current edge.
const voronoi_cell<coordinate_type> *cell1 = edge->cell();
@@ -272,8 +284,9 @@
// Edge is a segment or parabolic arc.
edge_points.push_back(edge->vertex0()->vertex());
edge_points.push_back(edge->vertex1()->vertex());
+ double max_dist = max_error * brect.min_len();
fill_intermediate_points(point1, point2, point3,
- edge_points);
+ edge_points, max_dist);
} else {
// Edge is a ray or line.
// Clip it with the bounding rectangle.
@@ -372,31 +385,26 @@
private:
voronoi_helper();
- // Find intermediate points of the parabola.
+ // Find intermediate points of the parabola. Number of points
+ // is defined by the value of max_dist parameter - maximum distance
+ // between parabola and line segments that interpolate it.
// Parabola is a locus of points equidistant from the point and segment
// sites. intermediate_points should contain two initial endpoints
// of the edge (voronoi vertices). Intermediate points are inserted
// between the given two endpoints.
+ // Max_dist is the maximum distance allowed between parabola and line
+ // segments that interpolate it.
static void fill_intermediate_points(
point_2d_type point_site,
point_2d_type segment_site_start,
point_2d_type segment_site_end,
- std::vector<point_2d_type> &intermediate_points) {
- // Get the number of intermediate points to represent parabola.
- int num_inter_points = get_intermediate_points_count(
- intermediate_points[0], intermediate_points[1]);
-
- // Check if there are any intermediate points.
- // In case bisector is a segment there won't be any.
- if (num_inter_points <= 0 ||
- point_site == segment_site_start ||
- point_site == segment_site_end) {
+ std::vector<point_2d_type> &intermediate_points,
+ double max_dist) {
+ // Check if bisector is a segment.
+ if (point_site == segment_site_start ||
+ point_site == segment_site_end)
return;
- }
-
- // Reserve space to avoid vector reallocations.
- intermediate_points.reserve(2 + num_inter_points);
-
+
// Apply the linear transformation to move start point of the
// segment to the point with coordinates (0, 0) and the direction
// of the segment to coincide the positive direction of the x-axis.
@@ -407,62 +415,86 @@
coordinate_type sqr_segment_length = segm_vec_x * segm_vec_x +
segm_vec_y * segm_vec_y;
- // Compute coordinates of the endpoints of the edge
+ // Compute x-coordinates of the endpoints of the edge
// in the transformed space.
- coordinate_type projection_start = get_point_projection(
- intermediate_points[0], segment_site_start, segment_site_end);
- coordinate_type projection_end = get_point_projection(
- intermediate_points[1], segment_site_start, segment_site_end);
- coordinate_type step = (projection_end - projection_start) *
- sqr_segment_length / (num_inter_points + 1);
+ coordinate_type projection_start = sqr_segment_length *
+ get_point_projection(intermediate_points[0],
+ segment_site_start,
+ segment_site_end);
+ coordinate_type projection_end = sqr_segment_length *
+ get_point_projection(intermediate_points[1],
+ segment_site_start,
+ segment_site_end);
// Compute parabola parameterers in the transformed space.
// Parabola has next representation:
- // f(x) = ((x-point_rot_x)^2 + point_rot_y^2) / (2.0*point_rot_y).
+ // f(x) = ((x-rot_x)^2 + rot_y^2) / (2.0*rot_y).
coordinate_type point_vec_x = point_site.x() -
segment_site_start.x();
coordinate_type point_vec_y = point_site.y() -
segment_site_start.y();
- coordinate_type point_rot_x = segm_vec_x * point_vec_x +
- segm_vec_y * point_vec_y;
- coordinate_type point_rot_y = segm_vec_x * point_vec_y -
- segm_vec_y * point_vec_x;
-
- // Compute the x-coordinate of the first intermediate point.
- coordinate_type projection_cur_x = projection_start *
- sqr_segment_length + step;
+ coordinate_type rot_x = segm_vec_x * point_vec_x +
+ segm_vec_y * point_vec_y;
+ coordinate_type rot_y = segm_vec_x * point_vec_y -
+ segm_vec_y * point_vec_x;
- // Temporary remove the last point.
- point_2d_type last_point = intermediate_points.back();
+ // Save the last point.
+ point_2d_type last_point = intermediate_points[1];
intermediate_points.pop_back();
- // Generate intermediate points.
- for (int i = 0; i < num_inter_points;
- i++, projection_cur_x += step) {
- // Compute coordinates of the intermediate points
- // in the transformed space.
- coordinate_type inter_rot_x = projection_cur_x;
- coordinate_type inter_rot_y =
- ((inter_rot_x - point_rot_x) *
- (inter_rot_x - point_rot_x) +
- point_rot_y * point_rot_y) / (2.0 * point_rot_y);
-
- // Compute coordinates of the intermediates points
- // in the initial space.
- coordinate_type new_point_x =
- (segm_vec_x * inter_rot_x - segm_vec_y * inter_rot_y) /
- sqr_segment_length + segment_site_start.x();
- coordinate_type new_point_y =
- (segm_vec_x * inter_rot_y + segm_vec_y * inter_rot_x) /
- sqr_segment_length + segment_site_start.y();
-
- // Append a new intermediate point to the end.
- intermediate_points.push_back(
- make_point_2d(new_point_x, new_point_y));
+ // Use stack to avoid recursion.
+ std::stack< coordinate_type > point_stack;
+ point_stack.push(projection_end);
+ coordinate_type cur_x = projection_start;
+ coordinate_type cur_y = parabola_y(cur_x, rot_x, rot_y);
+
+ // Adjust max_dist parameter in the transformed space.
+ max_dist *= max_dist * sqr_segment_length;
+
+ while (!point_stack.empty()) {
+ coordinate_type new_x = point_stack.top();
+ coordinate_type new_y = parabola_y(new_x, rot_x, rot_y);
+
+ // Compute coordinates of the point of the parabola that is
+ // furthest from the current line segment.
+ coordinate_type mid_x = (new_y - cur_y) / (new_x - cur_x) *
+ rot_y + rot_x;
+ coordinate_type mid_y = parabola_y(mid_x, rot_x, rot_y);
+
+ // Compute maximum distance between the given parabolic arc
+ // and line segment that interpolates it.
+ double dist = (new_y - cur_y) * (mid_x - cur_x) -
+ (new_x - cur_x) * (mid_y - cur_y);
+ dist = dist * dist / ((new_y - cur_y) * (new_y - cur_y) +
+ (new_x - cur_x) * (new_x - cur_x));
+ if (dist <= max_dist) {
+ // Distance between parabola and line segment is
+ // not greater than max_dist.
+ point_stack.pop();
+ coordinate_type inter_x =
+ (segm_vec_x * new_x - segm_vec_y * new_y) /
+ sqr_segment_length + segment_site_start.x();
+ coordinate_type inter_y =
+ (segm_vec_x * new_y + segm_vec_y * new_x) /
+ sqr_segment_length + segment_site_start.y();
+ intermediate_points.push_back(
+ make_point_2d(inter_x, inter_y));
+ cur_x = new_x;
+ cur_y = new_y;
+ } else {
+ point_stack.push(mid_x);
+ }
}
- // Append last point to the end.
- intermediate_points.push_back(last_point);
+ // Update last point.
+ intermediate_points.back() = last_point;
+ }
+
+ // Compute y(x) = ((x - a) * (x - a) + b * b) / (2 * b).
+ static coordinate_type parabola_y(coordinate_type x,
+ coordinate_type a,
+ coordinate_type b) {
+ return ((x - a) * (x - a) + b * b) / (2.0 * b);
}
// Check whether extent is compatible with the edge type.
@@ -529,18 +561,6 @@
segment_vec_y * point_vec_y;
return vec_dot / sqr_segment_length;
}
-
- // Return the number of the intermediate points to fill parabolic
- // arc with. Arc is represented by it's two endpoints. This doesn't
- // have any strong math background. Parameters were chosen
- // experimentally.
- static int get_intermediate_points_count(const point_2d_type &point1,
- const point_2d_type &point2) {
- coordinate_type vec_x = point1.x() - point2.x();
- coordinate_type vec_y = point1.y() - point2.y();
- coordinate_type sqr_len = vec_x * vec_x + vec_y * vec_y;
- return static_cast<int>(log(sqr_len) * 3.4 + 1E-6);
- }
};
// Represents voronoi cell.
Modified: sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_sweepline.hpp
==============================================================================
--- sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_sweepline.hpp (original)
+++ sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_sweepline.hpp 2011-01-19 14:03:45 EST (Wed, 19 Jan 2011)
@@ -16,6 +16,7 @@
#include <list>
#include <map>
#include <queue>
+#include <stack>
#include <vector>
#ifdef USE_MULTI_PRECISION_LIBRARY
@@ -23,6 +24,8 @@
#include "gmpxx.h"
#endif
+#define BOOST_POLYGON_UNSIGNED_LONG_LONG unsigned long long
+
#include "voronoi_output.hpp"
#include "detail/voronoi_formation.hpp"
Modified: sandbox/SOC/2010/sweepline/libs/sweepline/example/voronoi_visualizer.cpp
==============================================================================
--- sandbox/SOC/2010/sweepline/libs/sweepline/example/voronoi_visualizer.cpp (original)
+++ sandbox/SOC/2010/sweepline/libs/sweepline/example/voronoi_visualizer.cpp 2011-01-19 14:03:45 EST (Wed, 19 Jan 2011)
@@ -118,7 +118,7 @@
if (!it->is_primary() && primary_edges_only_)
continue;
std::vector<point_2d_type> temp_v =
- voronoi_helper<coordinate_type>::get_intermediate_points(&(*it), brect_);
+ voronoi_helper<coordinate_type>::get_intermediate_points(&(*it), brect_, 1E-3);
for (int i = 0; i < static_cast<int>(temp_v.size()) - 1; i++) {
glVertex2f(temp_v[i].x(), temp_v[i].y());
glVertex2f(temp_v[i+1].x(), temp_v[i+1].y());
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