voronoi_advanced_tutorial.cpp 4.5 KB

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