Loading...
Searching...
No Matches
random_point_generators.h
/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
* See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
* Author(s): Clement Jamin
*
* Copyright (C) 2016 Inria
*
* Modification(s):
* - 2019/08 Vincent Rouvreau: Fix issue #10 for CGAL
* - YYYY/MM Author: Description of the modification
*/
#ifndef RANDOM_POINT_GENERATORS_H_
#define RANDOM_POINT_GENERATORS_H_
#include <CGAL/number_utils.h>
#include <CGAL/Random.h>
#include <CGAL/point_generators_d.h>
#include <CGAL/version.h> // for CGAL_VERSION_NR
#include <vector> // for vector<>
#include <boost/math/constants/constants.hpp> // for pi constant
// Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10
#if CGAL_VERSION_NR < 1041101000
# error random_point_generators is only available for CGAL >= 4.11
#endif
namespace Gudhi {
// Note: All these functions have been tested with the CGAL::Epick_d kernel
// construct_point: dim 2
template <typename Kernel>
typename Kernel::Point_d construct_point(const Kernel &k,
typename Kernel::FT x1, typename Kernel::FT x2) {
typename Kernel::FT tab[2];
tab[0] = x1;
tab[1] = x2;
return k.construct_point_d_object()(2, &tab[0], &tab[2]);
}
// construct_point: dim 3
template <typename Kernel>
typename Kernel::Point_d construct_point(const Kernel &k,
typename Kernel::FT x1, typename Kernel::FT x2, typename Kernel::FT x3) {
typename Kernel::FT tab[3];
tab[0] = x1;
tab[1] = x2;
tab[2] = x3;
return k.construct_point_d_object()(3, &tab[0], &tab[3]);
}
// construct_point: dim 4
template <typename Kernel>
typename Kernel::Point_d construct_point(const Kernel &k,
typename Kernel::FT x1, typename Kernel::FT x2, typename Kernel::FT x3,
typename Kernel::FT x4) {
typename Kernel::FT tab[4];
tab[0] = x1;
tab[1] = x2;
tab[2] = x3;
tab[3] = x4;
return k.construct_point_d_object()(4, &tab[0], &tab[4]);
}
// construct_point: dim 5
template <typename Kernel>
typename Kernel::Point_d construct_point(const Kernel &k,
typename Kernel::FT x1, typename Kernel::FT x2, typename Kernel::FT x3,
typename Kernel::FT x4, typename Kernel::FT x5) {
typename Kernel::FT tab[5];
tab[0] = x1;
tab[1] = x2;
tab[2] = x3;
tab[3] = x4;
tab[4] = x5;
return k.construct_point_d_object()(5, &tab[0], &tab[5]);
}
// construct_point: dim 6
template <typename Kernel>
typename Kernel::Point_d construct_point(const Kernel &k,
typename Kernel::FT x1, typename Kernel::FT x2, typename Kernel::FT x3,
typename Kernel::FT x4, typename Kernel::FT x5, typename Kernel::FT x6) {
typename Kernel::FT tab[6];
tab[0] = x1;
tab[1] = x2;
tab[2] = x3;
tab[3] = x4;
tab[4] = x5;
tab[5] = x6;
return k.construct_point_d_object()(6, &tab[0], &tab[6]);
}
template <typename Kernel>
std::vector<typename Kernel::Point_d> generate_points_on_plane(std::size_t num_points, int intrinsic_dim,
int ambient_dim,
double coord_min = -5., double coord_max = 5.) {
typedef typename Kernel::Point_d Point;
typedef typename Kernel::FT FT;
Kernel k;
CGAL::Random rng;
std::vector<Point> points;
points.reserve(num_points);
for (std::size_t i = 0; i < num_points;) {
std::vector<FT> pt(ambient_dim, FT(0));
for (int j = 0; j < intrinsic_dim; ++j)
pt[j] = rng.get_double(coord_min, coord_max);
Point p = k.construct_point_d_object()(ambient_dim, pt.begin(), pt.end());
points.push_back(p);
++i;
}
return points;
}
template <typename Kernel>
std::vector<typename Kernel::Point_d> generate_points_on_moment_curve(std::size_t num_points, int dim,
typename Kernel::FT min_x,
typename Kernel::FT max_x) {
typedef typename Kernel::Point_d Point;
typedef typename Kernel::FT FT;
Kernel k;
CGAL::Random rng;
std::vector<Point> points;
points.reserve(num_points);
for (std::size_t i = 0; i < num_points;) {
FT x = rng.get_double(min_x, max_x);
std::vector<FT> coords;
coords.reserve(dim);
for (int p = 1; p <= dim; ++p)
coords.push_back(std::pow(CGAL::to_double(x), p));
Point p = k.construct_point_d_object()(
dim, coords.begin(), coords.end());
points.push_back(p);
++i;
}
return points;
}
// R = big radius, r = small radius
template <typename Kernel/*, typename TC_basis*/>
std::vector<typename Kernel::Point_d> generate_points_on_torus_3D(std::size_t num_points, double R, double r,
bool uniform = false) {
using namespace boost::math::double_constants;
typedef typename Kernel::Point_d Point;
typedef typename Kernel::FT FT;
Kernel k;
CGAL::Random rng;
// if uniform
std::size_t num_lines = (std::size_t)sqrt(num_points);
std::vector<Point> points;
points.reserve(num_points);
for (std::size_t i = 0; i < num_points;) {
FT u, v;
if (uniform) {
std::size_t k1 = i / num_lines;
std::size_t k2 = i % num_lines;
u = two_pi * k1 / num_lines;
v = two_pi * k2 / num_lines;
} else {
u = rng.get_double(0, two_pi);
v = rng.get_double(0, two_pi);
}
Point p = construct_point(k,
(R + r * std::cos(u)) * std::cos(v),
(R + r * std::cos(u)) * std::sin(v),
r * std::sin(u));
points.push_back(p);
++i;
}
return points;
}
// "Private" function used by generate_points_on_torus_d
template <typename Kernel, typename OutputIterator>
static void generate_grid_points_on_torus_d(const Kernel &k, int dim, std::size_t num_slices,
OutputIterator out,
double radius_noise_percentage = 0.,
std::vector<typename Kernel::FT> current_point =
std::vector<typename Kernel::FT>()) {
using namespace boost::math::double_constants;
CGAL::Random rng;
int point_size = static_cast<int>(current_point.size());
if (point_size == 2 * dim) {
*out++ = k.construct_point_d_object()(point_size, current_point.begin(), current_point.end());
} else {
for (std::size_t slice_idx = 0; slice_idx < num_slices; ++slice_idx) {
double radius_noise_ratio = 1.;
if (radius_noise_percentage > 0.) {
radius_noise_ratio = rng.get_double(
(100. - radius_noise_percentage) / 100.,
(100. + radius_noise_percentage) / 100.);
}
std::vector<typename Kernel::FT> cp2 = current_point;
double alpha = two_pi * slice_idx / num_slices;
cp2.push_back(radius_noise_ratio * std::cos(alpha));
cp2.push_back(radius_noise_ratio * std::sin(alpha));
generate_grid_points_on_torus_d(
k, dim, num_slices, out, radius_noise_percentage, cp2);
}
}
}
template <typename Kernel>
std::vector<typename Kernel::Point_d> generate_points_on_torus_d(std::size_t num_points, int dim, std::string sample = "random",
double radius_noise_percentage = 0.) {
using namespace boost::math::double_constants;
typedef typename Kernel::Point_d Point;
typedef typename Kernel::FT FT;
Kernel k;
CGAL::Random rng;
std::vector<Point> points;
points.reserve(num_points);
if (sample == "grid") {
std::size_t num_slices = (std::size_t)std::pow(num_points + .5, 1. / dim); // add .5 to avoid rounding down with numerical approximations
generate_grid_points_on_torus_d(
k, dim, num_slices, std::back_inserter(points), radius_noise_percentage);
} else {
for (std::size_t i = 0; i < num_points;) {
double radius_noise_ratio = 1.;
if (radius_noise_percentage > 0.) {
radius_noise_ratio = rng.get_double(
(100. - radius_noise_percentage) / 100.,
(100. + radius_noise_percentage) / 100.);
}
std::vector<typename Kernel::FT> pt;
pt.reserve(dim * 2);
for (int curdim = 0; curdim < dim; ++curdim) {
FT alpha = rng.get_double(0, two_pi);
pt.push_back(radius_noise_ratio * std::cos(alpha));
pt.push_back(radius_noise_ratio * std::sin(alpha));
}
Point p = k.construct_point_d_object()(pt.begin(), pt.end());
points.push_back(p);
++i;
}
}
return points;
}
template <typename Kernel>
std::vector<typename Kernel::Point_d> generate_points_on_sphere_d(std::size_t num_points, int dim, double radius,
double radius_noise_percentage = 0.) {
typedef typename Kernel::Point_d Point;
Kernel k;
CGAL::Random rng;
CGAL::Random_points_on_sphere_d<Point> generator(dim, radius);
std::vector<Point> points;
points.reserve(num_points);
for (std::size_t i = 0; i < num_points;) {
Point p = *generator++;
if (radius_noise_percentage > 0.) {
double radius_noise_ratio = rng.get_double(
(100. - radius_noise_percentage) / 100.,
(100. + radius_noise_percentage) / 100.);
typename Kernel::Point_to_vector_d k_pt_to_vec =
k.point_to_vector_d_object();
typename Kernel::Vector_to_point_d k_vec_to_pt =
k.vector_to_point_d_object();
typename Kernel::Scaled_vector_d k_scaled_vec =
k.scaled_vector_d_object();
p = k_vec_to_pt(k_scaled_vec(k_pt_to_vec(p), radius_noise_ratio));
}
points.push_back(p);
++i;
}
return points;
}
template <typename Kernel>
std::vector<typename Kernel::Point_d> generate_points_in_ball_d(std::size_t num_points, int dim, double radius) {
typedef typename Kernel::Point_d Point;
Kernel k;
CGAL::Random rng;
CGAL::Random_points_in_ball_d<Point> generator(dim, radius);
std::vector<Point> points;
points.reserve(num_points);
for (std::size_t i = 0; i < num_points;) {
Point p = *generator++;
points.push_back(p);
++i;
}
return points;
}
template <typename Kernel>
std::vector<typename Kernel::Point_d> generate_points_in_cube_d(std::size_t num_points, int dim, double radius) {
typedef typename Kernel::Point_d Point;
Kernel k;
CGAL::Random rng;
CGAL::Random_points_in_cube_d<Point> generator(dim, radius);
std::vector<Point> points;
points.reserve(num_points);
for (std::size_t i = 0; i < num_points;) {
Point p = *generator++;
points.push_back(p);
++i;
}
return points;
}
template <typename Kernel>
std::vector<typename Kernel::Point_d> generate_points_on_two_spheres_d(std::size_t num_points, int dim, double radius,
double distance_between_centers,
double radius_noise_percentage = 0.) {
typedef typename Kernel::FT FT;
typedef typename Kernel::Point_d Point;
typedef typename Kernel::Vector_d Vector;
Kernel k;
CGAL::Random rng;
CGAL::Random_points_on_sphere_d<Point> generator(dim, radius);
std::vector<Point> points;
points.reserve(num_points);
std::vector<FT> t(dim, FT(0));
t[0] = distance_between_centers;
Vector c1_to_c2(t.begin(), t.end());
for (std::size_t i = 0; i < num_points;) {
Point p = *generator++;
if (radius_noise_percentage > 0.) {
double radius_noise_ratio = rng.get_double(
(100. - radius_noise_percentage) / 100.,
(100. + radius_noise_percentage) / 100.);
typename Kernel::Point_to_vector_d k_pt_to_vec =
k.point_to_vector_d_object();
typename Kernel::Vector_to_point_d k_vec_to_pt =
k.vector_to_point_d_object();
typename Kernel::Scaled_vector_d k_scaled_vec =
k.scaled_vector_d_object();
p = k_vec_to_pt(k_scaled_vec(k_pt_to_vec(p), radius_noise_ratio));
}
typename Kernel::Translated_point_d k_transl =
k.translated_point_d_object();
Point p2 = k_transl(p, c1_to_c2);
points.push_back(p);
points.push_back(p2);
i += 2;
}
return points;
}
// Product of a 3-sphere and a circle => d = 3 / D = 5
template <typename Kernel>
std::vector<typename Kernel::Point_d> generate_points_on_3sphere_and_circle(std::size_t num_points,
double sphere_radius) {
using namespace boost::math::double_constants;
typedef typename Kernel::FT FT;
typedef typename Kernel::Point_d Point;
Kernel k;
CGAL::Random rng;
CGAL::Random_points_on_sphere_d<Point> generator(3, sphere_radius);
std::vector<Point> points;
points.reserve(num_points);
typename Kernel::Compute_coordinate_d k_coord =
k.compute_coordinate_d_object();
for (std::size_t i = 0; i < num_points;) {
Point p_sphere = *generator++; // First 3 coords
FT alpha = rng.get_double(0, two_pi);
std::vector<FT> pt(5);
pt[0] = k_coord(p_sphere, 0);
pt[1] = k_coord(p_sphere, 1);
pt[2] = k_coord(p_sphere, 2);
pt[3] = std::cos(alpha);
pt[4] = std::sin(alpha);
Point p(pt.begin(), pt.end());
points.push_back(p);
++i;
}
return points;
}
// a = big radius, b = small radius
template <typename Kernel>
std::vector<typename Kernel::Point_d> generate_points_on_klein_bottle_3D(std::size_t num_points, double a, double b,
bool uniform = false) {
using namespace boost::math::double_constants;
typedef typename Kernel::Point_d Point;
typedef typename Kernel::FT FT;
Kernel k;
CGAL::Random rng;
// if uniform
std::size_t num_lines = (std::size_t)sqrt(num_points);
std::vector<Point> points;
points.reserve(num_points);
for (std::size_t i = 0; i < num_points;) {
FT u, v;
if (uniform) {
std::size_t k1 = i / num_lines;
std::size_t k2 = i % num_lines;
u = two_pi * k1 / num_lines;
v = two_pi * k2 / num_lines;
} else {
u = rng.get_double(0, two_pi);
v = rng.get_double(0, two_pi);
}
double tmp = cos(u / 2) * sin(v) - sin(u / 2) * sin(2. * v);
Point p = construct_point(k,
(a + b * tmp) * cos(u),
(a + b * tmp) * sin(u),
b * (sin(u / 2) * sin(v) + cos(u / 2) * sin(2. * v)));
points.push_back(p);
++i;
}
return points;
}
// a = big radius, b = small radius
template <typename Kernel>
std::vector<typename Kernel::Point_d> generate_points_on_klein_bottle_4D(std::size_t num_points, double a, double b,
double noise = 0., bool uniform = false) {
using namespace boost::math::double_constants;
typedef typename Kernel::Point_d Point;
typedef typename Kernel::FT FT;
Kernel k;
CGAL::Random rng;
// if uniform
std::size_t num_lines = (std::size_t)sqrt(num_points);
std::vector<Point> points;
points.reserve(num_points);
for (std::size_t i = 0; i < num_points;) {
FT u, v;
if (uniform) {
std::size_t k1 = i / num_lines;
std::size_t k2 = i % num_lines;
u = two_pi * k1 / num_lines;
v = two_pi * k2 / num_lines;
} else {
u = rng.get_double(0, two_pi);
v = rng.get_double(0, two_pi);
}
Point p = construct_point(k,
(a + b * cos(v)) * cos(u) + (noise == 0. ? 0. : rng.get_double(0, noise)),
(a + b * cos(v)) * sin(u) + (noise == 0. ? 0. : rng.get_double(0, noise)),
b * sin(v) * cos(u / 2) + (noise == 0. ? 0. : rng.get_double(0, noise)),
b * sin(v) * sin(u / 2) + (noise == 0. ? 0. : rng.get_double(0, noise)));
points.push_back(p);
++i;
}
return points;
}
// a = big radius, b = small radius
template <typename Kernel>
std::vector<typename Kernel::Point_d>
generate_points_on_klein_bottle_variant_5D(
std::size_t num_points, double a, double b, bool uniform = false) {
using namespace boost::math::double_constants;
typedef typename Kernel::Point_d Point;
typedef typename Kernel::FT FT;
Kernel k;
CGAL::Random rng;
// if uniform
std::size_t num_lines = (std::size_t)sqrt(num_points);
std::vector<Point> points;
points.reserve(num_points);
for (std::size_t i = 0; i < num_points;) {
FT u, v;
if (uniform) {
std::size_t k1 = i / num_lines;
std::size_t k2 = i % num_lines;
u = two_pi * k1 / num_lines;
v = two_pi * k2 / num_lines;
} else {
u = rng.get_double(0, two_pi);
v = rng.get_double(0, two_pi);
}
FT x1 = (a + b * cos(v)) * cos(u);
FT x2 = (a + b * cos(v)) * sin(u);
FT x3 = b * sin(v) * cos(u / 2);
FT x4 = b * sin(v) * sin(u / 2);
FT x5 = x1 + x2 + x3 + x4;
Point p = construct_point(k, x1, x2, x3, x4, x5);
points.push_back(p);
++i;
}
return points;
}
} // namespace Gudhi
#endif // RANDOM_POINT_GENERATORS_H_
Gudhi
Gudhi namespace.
Definition SimplicialComplexForAlpha.h:14