random_point_generators.h
1/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
2 * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
3 * Author(s): Clement Jamin
4 *
5 * Copyright (C) 2016 Inria
6 *
7 * Modification(s):
8 * - 2019/08 Vincent Rouvreau: Fix issue #10 for CGAL
9 * - YYYY/MM Author: Description of the modification
10 */
11
12#ifndef RANDOM_POINT_GENERATORS_H_
13#define RANDOM_POINT_GENERATORS_H_
14
15#include <CGAL/number_utils.h>
16#include <CGAL/Random.h>
17#include <CGAL/point_generators_d.h>
18#include <CGAL/version.h> // for CGAL_VERSION_NR
19
20#include <vector> // for vector<>
21#include <boost/math/constants/constants.hpp> // for pi constant
22
23// Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10
24#if CGAL_VERSION_NR < 1041101000
25# error random_point_generators is only available for CGAL >= 4.11
26#endif
27
28namespace Gudhi {
29
31// Note: All these functions have been tested with the CGAL::Epick_d kernel
33
34// construct_point: dim 2
35
36template <typename Kernel>
37typename Kernel::Point_d construct_point(const Kernel &k,
38 typename Kernel::FT x1, typename Kernel::FT x2) {
39 typename Kernel::FT tab[2];
40 tab[0] = x1;
41 tab[1] = x2;
42 return k.construct_point_d_object()(2, &tab[0], &tab[2]);
43}
44
45// construct_point: dim 3
46
47template <typename Kernel>
48typename Kernel::Point_d construct_point(const Kernel &k,
49 typename Kernel::FT x1, typename Kernel::FT x2, typename Kernel::FT x3) {
50 typename Kernel::FT tab[3];
51 tab[0] = x1;
52 tab[1] = x2;
53 tab[2] = x3;
54 return k.construct_point_d_object()(3, &tab[0], &tab[3]);
55}
56
57// construct_point: dim 4
58
59template <typename Kernel>
60typename Kernel::Point_d construct_point(const Kernel &k,
61 typename Kernel::FT x1, typename Kernel::FT x2, typename Kernel::FT x3,
62 typename Kernel::FT x4) {
63 typename Kernel::FT tab[4];
64 tab[0] = x1;
65 tab[1] = x2;
66 tab[2] = x3;
67 tab[3] = x4;
68 return k.construct_point_d_object()(4, &tab[0], &tab[4]);
69}
70
71// construct_point: dim 5
72
73template <typename Kernel>
74typename Kernel::Point_d construct_point(const Kernel &k,
75 typename Kernel::FT x1, typename Kernel::FT x2, typename Kernel::FT x3,
76 typename Kernel::FT x4, typename Kernel::FT x5) {
77 typename Kernel::FT tab[5];
78 tab[0] = x1;
79 tab[1] = x2;
80 tab[2] = x3;
81 tab[3] = x4;
82 tab[4] = x5;
83 return k.construct_point_d_object()(5, &tab[0], &tab[5]);
84}
85
86// construct_point: dim 6
87
88template <typename Kernel>
89typename Kernel::Point_d construct_point(const Kernel &k,
90 typename Kernel::FT x1, typename Kernel::FT x2, typename Kernel::FT x3,
91 typename Kernel::FT x4, typename Kernel::FT x5, typename Kernel::FT x6) {
92 typename Kernel::FT tab[6];
93 tab[0] = x1;
94 tab[1] = x2;
95 tab[2] = x3;
96 tab[3] = x4;
97 tab[4] = x5;
98 tab[5] = x6;
99 return k.construct_point_d_object()(6, &tab[0], &tab[6]);
100}
101
102template <typename Kernel>
103std::vector<typename Kernel::Point_d> generate_points_on_plane(std::size_t num_points, int intrinsic_dim,
104 int ambient_dim,
105 double coord_min = -5., double coord_max = 5.) {
106 typedef typename Kernel::Point_d Point;
107 typedef typename Kernel::FT FT;
108 Kernel k;
109 CGAL::Random rng;
110 std::vector<Point> points;
111 points.reserve(num_points);
112 for (std::size_t i = 0; i < num_points;) {
113 std::vector<FT> pt(ambient_dim, FT(0));
114 for (int j = 0; j < intrinsic_dim; ++j)
115 pt[j] = rng.get_double(coord_min, coord_max);
116
117 Point p = k.construct_point_d_object()(ambient_dim, pt.begin(), pt.end());
118 points.push_back(p);
119 ++i;
120 }
121 return points;
122}
123
124template <typename Kernel>
125std::vector<typename Kernel::Point_d> generate_points_on_moment_curve(std::size_t num_points, int dim,
126 typename Kernel::FT min_x,
127 typename Kernel::FT max_x) {
128 typedef typename Kernel::Point_d Point;
129 typedef typename Kernel::FT FT;
130 Kernel k;
131 CGAL::Random rng;
132 std::vector<Point> points;
133 points.reserve(num_points);
134 for (std::size_t i = 0; i < num_points;) {
135 FT x = rng.get_double(min_x, max_x);
136 std::vector<FT> coords;
137 coords.reserve(dim);
138 for (int p = 1; p <= dim; ++p)
139 coords.push_back(std::pow(CGAL::to_double(x), p));
140 Point p = k.construct_point_d_object()(
141 dim, coords.begin(), coords.end());
142 points.push_back(p);
143 ++i;
144 }
145 return points;
146}
147
148
149// R = big radius, r = small radius
150template <typename Kernel/*, typename TC_basis*/>
151std::vector<typename Kernel::Point_d> generate_points_on_torus_3D(std::size_t num_points, double R, double r,
152 bool uniform = false) {
153 using namespace boost::math::double_constants;
154
155 typedef typename Kernel::Point_d Point;
156 typedef typename Kernel::FT FT;
157 Kernel k;
158 CGAL::Random rng;
159
160 // if uniform
161 std::size_t num_lines = (std::size_t)sqrt(num_points);
162
163 std::vector<Point> points;
164 points.reserve(num_points);
165 for (std::size_t i = 0; i < num_points;) {
166 FT u, v;
167 if (uniform) {
168 std::size_t k1 = i / num_lines;
169 std::size_t k2 = i % num_lines;
170 u = two_pi * k1 / num_lines;
171 v = two_pi * k2 / num_lines;
172 } else {
173 u = rng.get_double(0, two_pi);
174 v = rng.get_double(0, two_pi);
175 }
176 Point p = construct_point(k,
177 (R + r * std::cos(u)) * std::cos(v),
178 (R + r * std::cos(u)) * std::sin(v),
179 r * std::sin(u));
180 points.push_back(p);
181 ++i;
182 }
183 return points;
184}
185
186// "Private" function used by generate_points_on_torus_d
187template <typename Kernel, typename OutputIterator>
188static void generate_grid_points_on_torus_d(const Kernel &k, int dim, std::size_t num_slices,
189 OutputIterator out,
190 double radius_noise_percentage = 0.,
191 std::vector<typename Kernel::FT> current_point =
192 std::vector<typename Kernel::FT>()) {
193 using namespace boost::math::double_constants;
194
195 CGAL::Random rng;
196 int point_size = static_cast<int>(current_point.size());
197 if (point_size == 2 * dim) {
198 *out++ = k.construct_point_d_object()(point_size, current_point.begin(), current_point.end());
199 } else {
200 for (std::size_t slice_idx = 0; slice_idx < num_slices; ++slice_idx) {
201 double radius_noise_ratio = 1.;
202 if (radius_noise_percentage > 0.) {
203 radius_noise_ratio = rng.get_double(
204 (100. - radius_noise_percentage) / 100.,
205 (100. + radius_noise_percentage) / 100.);
206 }
207 std::vector<typename Kernel::FT> cp2 = current_point;
208 double alpha = two_pi * slice_idx / num_slices;
209 cp2.push_back(radius_noise_ratio * std::cos(alpha));
210 cp2.push_back(radius_noise_ratio * std::sin(alpha));
211 generate_grid_points_on_torus_d(
212 k, dim, num_slices, out, radius_noise_percentage, cp2);
213 }
214 }
215}
216
217template <typename Kernel>
218std::vector<typename Kernel::Point_d> generate_points_on_torus_d(std::size_t num_points, int dim, std::string sample = "random",
219 double radius_noise_percentage = 0.) {
220 using namespace boost::math::double_constants;
221
222 typedef typename Kernel::Point_d Point;
223 typedef typename Kernel::FT FT;
224 Kernel k;
225 CGAL::Random rng;
226
227 std::vector<Point> points;
228 points.reserve(num_points);
229 if (sample == "grid") {
230 std::size_t num_slices = (std::size_t)std::pow(num_points + .5, 1. / dim); // add .5 to avoid rounding down with numerical approximations
231 generate_grid_points_on_torus_d(
232 k, dim, num_slices, std::back_inserter(points), radius_noise_percentage);
233 } else {
234 for (std::size_t i = 0; i < num_points;) {
235 double radius_noise_ratio = 1.;
236 if (radius_noise_percentage > 0.) {
237 radius_noise_ratio = rng.get_double(
238 (100. - radius_noise_percentage) / 100.,
239 (100. + radius_noise_percentage) / 100.);
240 }
241 std::vector<typename Kernel::FT> pt;
242 pt.reserve(dim * 2);
243 for (int curdim = 0; curdim < dim; ++curdim) {
244 FT alpha = rng.get_double(0, two_pi);
245 pt.push_back(radius_noise_ratio * std::cos(alpha));
246 pt.push_back(radius_noise_ratio * std::sin(alpha));
247 }
248
249 Point p = k.construct_point_d_object()(pt.begin(), pt.end());
250 points.push_back(p);
251 ++i;
252 }
253 }
254 return points;
255}
256
257template <typename Kernel>
258std::vector<typename Kernel::Point_d> generate_points_on_sphere_d(std::size_t num_points, int dim, double radius,
259 double radius_noise_percentage = 0.) {
260 typedef typename Kernel::Point_d Point;
261 Kernel k;
262 CGAL::Random rng;
263 CGAL::Random_points_on_sphere_d<Point> generator(dim, radius);
264 std::vector<Point> points;
265 points.reserve(num_points);
266 for (std::size_t i = 0; i < num_points;) {
267 Point p = *generator++;
268 if (radius_noise_percentage > 0.) {
269 double radius_noise_ratio = rng.get_double(
270 (100. - radius_noise_percentage) / 100.,
271 (100. + radius_noise_percentage) / 100.);
272
273 typename Kernel::Point_to_vector_d k_pt_to_vec =
274 k.point_to_vector_d_object();
275 typename Kernel::Vector_to_point_d k_vec_to_pt =
276 k.vector_to_point_d_object();
277 typename Kernel::Scaled_vector_d k_scaled_vec =
278 k.scaled_vector_d_object();
279 p = k_vec_to_pt(k_scaled_vec(k_pt_to_vec(p), radius_noise_ratio));
280 }
281 points.push_back(p);
282 ++i;
283 }
284 return points;
285}
286
287template <typename Kernel>
288std::vector<typename Kernel::Point_d> generate_points_in_ball_d(std::size_t num_points, int dim, double radius) {
289 typedef typename Kernel::Point_d Point;
290 Kernel k;
291 CGAL::Random rng;
292 CGAL::Random_points_in_ball_d<Point> generator(dim, radius);
293 std::vector<Point> points;
294 points.reserve(num_points);
295 for (std::size_t i = 0; i < num_points;) {
296 Point p = *generator++;
297 points.push_back(p);
298 ++i;
299 }
300 return points;
301}
302
303template <typename Kernel>
304std::vector<typename Kernel::Point_d> generate_points_in_cube_d(std::size_t num_points, int dim, double radius) {
305 typedef typename Kernel::Point_d Point;
306 Kernel k;
307 CGAL::Random rng;
308 CGAL::Random_points_in_cube_d<Point> generator(dim, radius);
309 std::vector<Point> points;
310 points.reserve(num_points);
311 for (std::size_t i = 0; i < num_points;) {
312 Point p = *generator++;
313 points.push_back(p);
314 ++i;
315 }
316 return points;
317}
318
319template <typename Kernel>
320std::vector<typename Kernel::Point_d> generate_points_on_two_spheres_d(std::size_t num_points, int dim, double radius,
321 double distance_between_centers,
322 double radius_noise_percentage = 0.) {
323 typedef typename Kernel::FT FT;
324 typedef typename Kernel::Point_d Point;
325 typedef typename Kernel::Vector_d Vector;
326 Kernel k;
327 CGAL::Random rng;
328 CGAL::Random_points_on_sphere_d<Point> generator(dim, radius);
329 std::vector<Point> points;
330 points.reserve(num_points);
331
332 std::vector<FT> t(dim, FT(0));
333 t[0] = distance_between_centers;
334 Vector c1_to_c2(t.begin(), t.end());
335
336 for (std::size_t i = 0; i < num_points;) {
337 Point p = *generator++;
338 if (radius_noise_percentage > 0.) {
339 double radius_noise_ratio = rng.get_double(
340 (100. - radius_noise_percentage) / 100.,
341 (100. + radius_noise_percentage) / 100.);
342
343 typename Kernel::Point_to_vector_d k_pt_to_vec =
344 k.point_to_vector_d_object();
345 typename Kernel::Vector_to_point_d k_vec_to_pt =
346 k.vector_to_point_d_object();
347 typename Kernel::Scaled_vector_d k_scaled_vec =
348 k.scaled_vector_d_object();
349 p = k_vec_to_pt(k_scaled_vec(k_pt_to_vec(p), radius_noise_ratio));
350 }
351
352 typename Kernel::Translated_point_d k_transl =
353 k.translated_point_d_object();
354 Point p2 = k_transl(p, c1_to_c2);
355 points.push_back(p);
356 points.push_back(p2);
357 i += 2;
358 }
359 return points;
360}
361
362// Product of a 3-sphere and a circle => d = 3 / D = 5
363
364template <typename Kernel>
365std::vector<typename Kernel::Point_d> generate_points_on_3sphere_and_circle(std::size_t num_points,
366 double sphere_radius) {
367 using namespace boost::math::double_constants;
368
369 typedef typename Kernel::FT FT;
370 typedef typename Kernel::Point_d Point;
371 Kernel k;
372 CGAL::Random rng;
373 CGAL::Random_points_on_sphere_d<Point> generator(3, sphere_radius);
374 std::vector<Point> points;
375 points.reserve(num_points);
376
377 typename Kernel::Compute_coordinate_d k_coord =
378 k.compute_coordinate_d_object();
379 for (std::size_t i = 0; i < num_points;) {
380 Point p_sphere = *generator++; // First 3 coords
381
382 FT alpha = rng.get_double(0, two_pi);
383 std::vector<FT> pt(5);
384 pt[0] = k_coord(p_sphere, 0);
385 pt[1] = k_coord(p_sphere, 1);
386 pt[2] = k_coord(p_sphere, 2);
387 pt[3] = std::cos(alpha);
388 pt[4] = std::sin(alpha);
389 Point p(pt.begin(), pt.end());
390 points.push_back(p);
391 ++i;
392 }
393 return points;
394}
395
396// a = big radius, b = small radius
397template <typename Kernel>
398std::vector<typename Kernel::Point_d> generate_points_on_klein_bottle_3D(std::size_t num_points, double a, double b,
399 bool uniform = false) {
400 using namespace boost::math::double_constants;
401
402 typedef typename Kernel::Point_d Point;
403 typedef typename Kernel::FT FT;
404 Kernel k;
405 CGAL::Random rng;
406
407 // if uniform
408 std::size_t num_lines = (std::size_t)sqrt(num_points);
409
410 std::vector<Point> points;
411 points.reserve(num_points);
412 for (std::size_t i = 0; i < num_points;) {
413 FT u, v;
414 if (uniform) {
415 std::size_t k1 = i / num_lines;
416 std::size_t k2 = i % num_lines;
417 u = two_pi * k1 / num_lines;
418 v = two_pi * k2 / num_lines;
419 } else {
420 u = rng.get_double(0, two_pi);
421 v = rng.get_double(0, two_pi);
422 }
423 double tmp = cos(u / 2) * sin(v) - sin(u / 2) * sin(2. * v);
424 Point p = construct_point(k,
425 (a + b * tmp) * cos(u),
426 (a + b * tmp) * sin(u),
427 b * (sin(u / 2) * sin(v) + cos(u / 2) * sin(2. * v)));
428 points.push_back(p);
429 ++i;
430 }
431 return points;
432}
433
434// a = big radius, b = small radius
435template <typename Kernel>
436std::vector<typename Kernel::Point_d> generate_points_on_klein_bottle_4D(std::size_t num_points, double a, double b,
437 double noise = 0., bool uniform = false) {
438 using namespace boost::math::double_constants;
439
440 typedef typename Kernel::Point_d Point;
441 typedef typename Kernel::FT FT;
442 Kernel k;
443 CGAL::Random rng;
444
445 // if uniform
446 std::size_t num_lines = (std::size_t)sqrt(num_points);
447
448 std::vector<Point> points;
449 points.reserve(num_points);
450 for (std::size_t i = 0; i < num_points;) {
451 FT u, v;
452 if (uniform) {
453 std::size_t k1 = i / num_lines;
454 std::size_t k2 = i % num_lines;
455 u = two_pi * k1 / num_lines;
456 v = two_pi * k2 / num_lines;
457 } else {
458 u = rng.get_double(0, two_pi);
459 v = rng.get_double(0, two_pi);
460 }
461 Point p = construct_point(k,
462 (a + b * cos(v)) * cos(u) + (noise == 0. ? 0. : rng.get_double(0, noise)),
463 (a + b * cos(v)) * sin(u) + (noise == 0. ? 0. : rng.get_double(0, noise)),
464 b * sin(v) * cos(u / 2) + (noise == 0. ? 0. : rng.get_double(0, noise)),
465 b * sin(v) * sin(u / 2) + (noise == 0. ? 0. : rng.get_double(0, noise)));
466 points.push_back(p);
467 ++i;
468 }
469 return points;
470}
471
472
473// a = big radius, b = small radius
474
475template <typename Kernel>
476std::vector<typename Kernel::Point_d>
477generate_points_on_klein_bottle_variant_5D(
478 std::size_t num_points, double a, double b, bool uniform = false) {
479 using namespace boost::math::double_constants;
480
481 typedef typename Kernel::Point_d Point;
482 typedef typename Kernel::FT FT;
483 Kernel k;
484 CGAL::Random rng;
485
486 // if uniform
487 std::size_t num_lines = (std::size_t)sqrt(num_points);
488
489 std::vector<Point> points;
490 points.reserve(num_points);
491 for (std::size_t i = 0; i < num_points;) {
492 FT u, v;
493 if (uniform) {
494 std::size_t k1 = i / num_lines;
495 std::size_t k2 = i % num_lines;
496 u = two_pi * k1 / num_lines;
497 v = two_pi * k2 / num_lines;
498 } else {
499 u = rng.get_double(0, two_pi);
500 v = rng.get_double(0, two_pi);
501 }
502 FT x1 = (a + b * cos(v)) * cos(u);
503 FT x2 = (a + b * cos(v)) * sin(u);
504 FT x3 = b * sin(v) * cos(u / 2);
505 FT x4 = b * sin(v) * sin(u / 2);
506 FT x5 = x1 + x2 + x3 + x4;
507
508 Point p = construct_point(k, x1, x2, x3, x4, x5);
509 points.push_back(p);
510 ++i;
511 }
512 return points;
513}
514
515} // namespace Gudhi
516
517#endif // RANDOM_POINT_GENERATORS_H_