libs/polygon/example/voronoi_advanced_tutorial.cpp

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