1 // Copyright (c) 2014
2 // Utrecht University (The Netherlands),
3 // ETH Zurich (Switzerland),
4 // INRIA Sophia-Antipolis (France),
5 // Max-Planck-Institute Saarbruecken (Germany),
6 // and Tel-Aviv University (Israel). All rights reserved.
7 //
8 // This file is part of CGAL (www.cgal.org)
9 //
10 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Generator/include/CGAL/random_convex_hull_in_disc_2.h $
11 // $Id: random_convex_hull_in_disc_2.h 3b1dfc7 2021-03-24T08:44:17+01:00 Simon Giraudot
12 // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
13 //
14 //
15 // Author(s) : Remy Thomasse <remy.thomasse@inria.fr>
16
17 #ifndef CGAL_RANDOM_CONVEX_HULL_DISC_H
18 #define CGAL_RANDOM_CONVEX_HULL_DISC_H 1
19 #include <boost/random.hpp>
20 #include <CGAL/Polygon_2_algorithms.h>
21 #include <CGAL/function_objects.h>
22 #include <CGAL/copy_n.h>
23 #include <CGAL/number_type_config.h>
24 #include <CGAL/double.h>
25 #include <list>
26
27 namespace CGAL {
28 namespace internal {
29 template <class P, class Traits>
30 struct compare_points_angle {
operatorcompare_points_angle31 bool operator()(const P& p, const P& q) {
32 P zero(0,0);
33 Traits traits;
34 typedef typename Traits::Orientation_2 Orientation_2;
35 Orientation_2 orientation_2=traits.orientation_2_object();
36 typedef typename Traits::Compare_y_2 Compare_y_2;
37 Compare_y_2 compare_y_2=traits.compare_y_2_object();
38
39 if (compare_y_2(p, zero) == LARGER ){
40 if (compare_y_2(q, zero)==LARGER)
41 return (orientation_2(zero,p,q)==LEFT_TURN);
42 else
43 return false;
44
45 } else {
46 if (compare_y_2(q,zero)==LARGER)
47 return true;
48 else
49 return (orientation_2(zero,p,q)==LEFT_TURN);
50 }
51 }
52 };
53
54 //////////////////////////////////////
55 template <class P, class GEN>
generate_points_annulus(long n,double a,double b,double small_radius,double big_radius,std::list<P> & l,GEN & gen)56 void generate_points_annulus(long n, double a, double b, double small_radius,
57 double big_radius, std::list<P>& l,
58 GEN& gen) { // generate n points between a and b
59 if (n > 1) {
60 boost::binomial_distribution<long> bin_distribution(n, .5);
61 boost::variate_generator<GEN&, boost::binomial_distribution<long> >
62 g(gen, bin_distribution);
63 long nb = g();
64 generate_points_annulus(nb, a, (a + b) / 2.0, small_radius, big_radius, l,
65 gen);
66 generate_points_annulus(n - nb, (a + b) / 2.0, b, small_radius, big_radius,
67 l, gen);
68 }
69 if (n == 1) // generation of a point
70 {
71
72 #if BOOST_VERSION < 104700
73
74 boost::uniform_real<double> random_squared_radius_distribution(
75 small_radius * small_radius / (big_radius * big_radius), 1);
76 boost::uniform_real<double> random_angle_distribution(a, b);
77 boost::variate_generator<
78 GEN&, boost::uniform_real<double> > random_angle(gen, random_angle_distribution);
79 boost::variate_generator<
80 GEN&, boost::uniform_real<double> > random_squared_radius(gen, random_squared_radius_distribution);
81
82 #else
83
84 boost::random::uniform_real_distribution<double> random_squared_radius_distribution(
85 small_radius * small_radius / (big_radius * big_radius), 1);
86
87 boost::random::uniform_real_distribution<double> random_angle_distribution(a, b);
88 boost::random::variate_generator<
89 GEN&, boost::random::uniform_real_distribution<double> > random_angle(gen, random_angle_distribution);
90 boost::random::variate_generator<
91 GEN&, boost::random::uniform_real_distribution<double> > random_squared_radius(gen, random_squared_radius_distribution);
92
93 #endif
94
95 double alpha = random_angle();
96 double r = big_radius * std::sqrt(random_squared_radius());
97 typedef Creator_uniform_2<double, P> Creator;
98 Creator creator;
99 typedef typename Creator::argument_type T;
100 l.push_back(creator(T(r * cos(alpha)), T(r * std::sin(alpha))));
101 }
102 }
103
104 template <class P>
Cyclic_increment(typename std::list<P>::iterator & it,std::list<P> & l)105 void Cyclic_increment(typename std::list<P>::iterator& it,
106 std::list<P>& l) {
107 ++it;
108 if (it == l.end()) {
109 it = l.begin();
110 }
111 }
112 //////////////////////////////////////////////////////////////////////////////
113 template <class P, class Traits>
Graham_without_sort_2(std::list<P> & l,const Traits & traits)114 void Graham_without_sort_2(std::list<P>& l, const Traits& traits) {
115 if (l.size() > 3) {
116 //typedef typename Traits::Left_turn_2 Left_turn;
117 //Left_turn left_turn = traits.left_turn_2_object();
118 typedef typename Traits::Orientation_2 Orientation_2;
119 typedef typename std::list<P>::iterator Iterator;
120 Orientation_2 orientation_2=traits.orientation_2_object();
121
122 typedef typename Traits::Compare_x_2 Compare_x_2;
123 Compare_x_2 compare_x_2=traits.compare_x_2_object();
124
125 Iterator pmin = l.begin();
126 for (Iterator it = l.begin(); it != l.end(); ++it) {
127 if (compare_x_2(*pmin, *it) == LARGER){
128 pmin = it;
129 }
130 } //*pmin is the extremal point on the left
131 Iterator u = pmin;
132 Iterator u_next = u;
133 Cyclic_increment(u_next, l);
134
135 Iterator u_next_next = u_next;
136 Cyclic_increment(u_next_next, l);
137
138 while (u_next != pmin) {
139
140 if (orientation_2(*u,*u_next,*u_next_next)==LEFT_TURN){
141 Cyclic_increment(u, l);
142 Cyclic_increment(u_next, l);
143 Cyclic_increment(u_next_next, l);
144 } else {
145 u_next = l.erase(u_next);
146 if (u_next == l.end()) u_next = l.begin();
147 if (u != pmin) {
148 u_next = u;
149 if (u == l.begin()) {
150 u = l.end();
151 }
152 --u;
153 } else {
154 u_next_next = u_next;
155 Cyclic_increment(u_next_next, l);
156 }
157 }
158 }
159 }
160 }
161
162 //////////////////////////////////////////////////////////////////////////////
163 template <class GEN, class Traits>
164
165 void random_convex_hull_in_disc_2(std::size_t n, double radius, std::list<typename Traits::Point_2>& l,
166 GEN& gen, const Traits& traits,
167 bool fast = true) {
168 CGAL_precondition(n >= 3);
169 typedef typename Traits::Point_2 P;
170 typedef typename Traits::FT FT;
171 std::size_t simulated_points = 0;
172 std::size_t generated_points = 0;
173 P zero(0, 0);
174 FT squared_radius = radius * radius;
175 typedef typename Traits::Compare_y_2 Compare_y_2;
176 Compare_y_2 compare_y_2=traits.compare_y_2_object();
177
178
179
180 do { // Initialization
181 long init =
182 static_cast<long>((std::min)(static_cast<std::size_t>(100), n - simulated_points));
183
184 generate_points_annulus(init, -CGAL_PI, CGAL_PI, 0, radius, l,
185 gen);
186
187 simulated_points += init;
188 generated_points += init;
189
190 Graham_without_sort_2(l, traits);
191 } while ((bounded_side_2(l.begin(), l.end(), zero, traits) !=
192 ON_BOUNDED_SIDE) &&
193 (simulated_points < n)); // initialization such that 0 in P_n
194
195 std::size_t T = n;
196 if (!fast) T = static_cast<std::size_t>(std::floor(n / CGAL::square(std::log(static_cast<double>(n)))));
197
198 while (simulated_points < n) {
199 // l is a list coming from a convex hull operation. we are moving the
200 // points s.t the angles are from -pi to pi.
201 {
202 typename std::list<P>::iterator it = l.begin();
203 while (compare_y_2(*it,zero) == LARGER){
204 l.push_back(*it);
205 l.pop_front();
206 it = l.begin();
207 }
208 it = l.end();
209 --it; // last element
210 while (compare_y_2(*it,zero) == SMALLER){
211 l.push_front(*it);
212 l.pop_back();
213 it = l.end();
214 --it; // last element
215 }
216 }
217 FT squared_small_radius = squared_radius;
218
219 {
220 typename std::list<P>::iterator it = l.begin();
221 typename std::list<P>::iterator it2=++l.begin();
222 for (; it != l.end();
223 ++it, Cyclic_increment(it2, l)) { // computation of annulus
224 typename Traits::Segment_2 s(*it, *it2);
225 FT temp=squared_distance(s,zero);
226 if ( compare(squared_small_radius,temp) == LARGER ) squared_small_radius=temp;
227 }
228 } // squared_small_radius=squared small radius of the annulus
229
230 FT p_disc = squared_small_radius / squared_radius;
231 long nb;
232 if (simulated_points < T) {
233 nb = static_cast<long>((std::min)(simulated_points, n - simulated_points));
234 } else {
235 nb = static_cast<long>((std::min)(T, n - simulated_points));
236 }
237 boost::binomial_distribution<long> dbin(nb, to_double(p_disc));
238 boost::variate_generator<
239 GEN&, boost::binomial_distribution<long> > bin(gen, dbin);
240
241 // How many points are falling in the small disc and wont be generated:
242 long k_disc = bin();
243 simulated_points += k_disc;
244
245 std::list<P> m;
246 generate_points_annulus(nb - k_disc, -CGAL_PI, CGAL_PI,
247 std::sqrt(to_double(squared_small_radius)), radius,
248 m, gen);
249 l.merge(m, compare_points_angle<P, Traits>());
250 generated_points += nb - k_disc;
251 simulated_points += nb - k_disc;
252 m.clear();
253 Graham_without_sort_2(l, traits);
254 }
255 }
256
257 } // namespace CGAL::internal
258
259 ///
260
261 template <class OutputIterator, class Traits, class Generator>
262 void random_convex_hull_in_disc_2(std::size_t n, double radius, Generator& gen,
263 OutputIterator it, const Traits& traits,
264 bool fast = true) {
265 typedef typename Traits::Point_2 Points;
266 std::list<Points> l;
267 internal::random_convex_hull_in_disc_2(n, radius, l, gen, traits, fast);
268 std::copy_n(l.begin(),l.size(),it);
269 }
270
271 } // namespace CGAL
272 #endif
273