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