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