1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
|
// Boost.Polygon library voronoi_advanced_tutorial.cpp file
// Copyright Andrii Sydorchuk 2010-2012.
// 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 <cstdio>
#include <ctime>
#include <string>
// This will work properly only with GCC compiler.
#include <ieee754.h>
typedef long double fpt80;
// Random generators and distributions.
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/uniform_int_distribution.hpp>
#include <boost/timer/timer.hpp>
#include <boost/polygon/voronoi.hpp>
using namespace boost::polygon;
struct my_ulp_comparison {
enum Result {
LESS = -1,
EQUAL = 0,
MORE = 1
};
Result operator()(fpt80 a, fpt80 b, unsigned int maxUlps) const {
if (a == b)
return EQUAL;
if (a > b) {
Result res = operator()(b, a, maxUlps);
if (res == EQUAL) return res;
return (res == LESS) ? MORE : LESS;
}
ieee854_long_double lhs, rhs;
lhs.d = a;
rhs.d = b;
if (lhs.ieee.negative ^ rhs.ieee.negative)
return lhs.ieee.negative ? LESS : MORE;
boost::uint64_t le = lhs.ieee.exponent; le =
(le << 32) + lhs.ieee.mantissa0;
boost::uint64_t re = rhs.ieee.exponent; re =
(re << 32) + rhs.ieee.mantissa0;
if (lhs.ieee.negative) {
if (le - 1 > re)
return LESS;
le = (le == re) ? 0 : 1;
le = (le << 32) + lhs.ieee.mantissa1;
re = rhs.ieee.mantissa1;
return (re + maxUlps < le) ? LESS : EQUAL;
} else {
if (le + 1 < re)
return LESS;
le = lhs.ieee.mantissa0;
re = (le == re) ? 0 : 1;
re = (re << 32) + rhs.ieee.mantissa1;
return (le + maxUlps < re) ? LESS : EQUAL;
}
}
};
struct my_fpt_converter {
template <typename T>
fpt80 operator()(const T& that) const {
return static_cast<fpt80>(that);
}
template <std::size_t N>
fpt80 operator()(const typename detail::extended_int<N> &that) const {
fpt80 result = 0.0;
for (std::size_t i = 1; i <= (std::min)((std::size_t)3, that.size()); ++i) {
if (i != 1)
result *= static_cast<fpt80>(0x100000000ULL);
result += that.chunks()[that.size() - i];
}
return (that.count() < 0) ? -result : result;
}
};
// Voronoi diagram traits.
struct my_voronoi_diagram_traits {
typedef fpt80 coordinate_type;
typedef voronoi_cell<coordinate_type> cell_type;
typedef voronoi_vertex<coordinate_type> vertex_type;
typedef voronoi_edge<coordinate_type> edge_type;
class vertex_equality_predicate_type {
public:
enum { ULPS = 128 };
bool operator()(const vertex_type &v1, const vertex_type &v2) const {
return (ulp_cmp(v1.x(), v2.x(), ULPS) == my_ulp_comparison::EQUAL &&
ulp_cmp(v1.y(), v2.y(), ULPS) == my_ulp_comparison::EQUAL);
}
private:
my_ulp_comparison ulp_cmp;
};
};
// Voronoi ctype traits for 48-bit signed integer input coordinates.
struct my_voronoi_ctype_traits {
typedef boost::int64_t int_type;
typedef detail::extended_int<3> int_x2_type;
typedef detail::extended_int<3> uint_x2_type;
typedef detail::extended_int<128> big_int_type;
typedef fpt80 fpt_type;
typedef fpt80 efpt_type;
typedef my_ulp_comparison ulp_cmp_type;
typedef my_fpt_converter to_fpt_converter_type;
typedef my_fpt_converter to_efpt_converter_type;
};
const unsigned int GENERATED_POINTS = 100;
const boost::int64_t MAX = 0x1000000000000LL;
int main() {
boost::mt19937_64 gen(std::time(0));
boost::random::uniform_int_distribution<boost::int64_t> distr(-MAX, MAX-1);
voronoi_builder<boost::int64_t, my_voronoi_ctype_traits> vb;
for (size_t i = 0; i < GENERATED_POINTS; ++i) {
boost::int64_t x = distr(gen);
boost::int64_t y = distr(gen);
vb.insert_point(x, y);
}
printf("Constructing Voronoi diagram of %d points...\n", GENERATED_POINTS);
boost::timer::cpu_timer t;
voronoi_diagram<fpt80, my_voronoi_diagram_traits> vd;
t.start();
vb.construct(&vd);
boost::timer::cpu_times times = t.elapsed();
std::string ftime = boost::timer::format(times, 5, "%w");
printf("Construction done in: %s seconds.\n", ftime.c_str());
printf("Resulting Voronoi graph has the following stats:\n");
printf("Number of Voronoi cells: %lu.\n", vd.num_cells());
printf("Number of Voronoi vertices: %lu.\n", vd.num_vertices());
printf("Number of Voronoi edges: %lu.\n", vd.num_edges());
return 0;
}
|
