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