discretization_helpers.hpp

Go to the documentation of this file.

#ifndef MEDUSA_BITS_DOMAINS_DISCRETIZATION_HELPERS_HPP_
#define MEDUSA_BITS_DOMAINS_DISCRETIZATION_HELPERS_HPP_
 
#include <medusa/Config.hpp>
#include <medusa/bits/types/Range.hpp>
#include <Eigen/Core>
#include <random>
 
namespace mm {
 
namespace discretization_helpers {
 
template <typename vec_t, typename func_t>
Range<vec_t> discretizeLineWithDensity(const vec_t& p, const vec_t& q, const func_t& delta_r) {
    typedef typename vec_t::scalar_t scalar_t;
    scalar_t dist = (q-p).norm();
    assert_msg(dist >= 1e-15, "Points %s and %s are equal.", p, q);
    vec_t dir = (q-p) / dist;
    scalar_t cur_dist = 0;
    Range<scalar_t> dists;
    while (cur_dist < dist) {
        cur_dist += delta_r(p+cur_dist*dir);
        dists.push_back(cur_dist);
    }
    dists.pop_back();
    if (dists.size() <= 2) return {};
    // decide whether to include the last point
    scalar_t spacing_with_last_point = dist - dists.back();
    scalar_t spacing_without_last_point = dist - dists[dists.size()-2];
    scalar_t desired_spacing = delta_r(q);
    scalar_t error_with = std::abs(1 - spacing_with_last_point / desired_spacing);
    scalar_t error_without = std::abs(1 - spacing_without_last_point / desired_spacing);
    if (error_without < error_with) {
        dists.pop_back();
    }
 
    Range<vec_t> points(dists.size());
    for (int i = 0; i < dists.size(); ++i) {
        points[i] = p + dists[i]*dir;
    }
    return points;
}
 
template <typename scalar_t, int dim>
struct SphereDiscretization {
    template <typename generator_t>
    static std::vector<Eigen::Matrix<scalar_t, dim, 1>,
                       Eigen::aligned_allocator<Eigen::Matrix<scalar_t, dim, 1>>> construct(
            scalar_t radius, int num_samples, generator_t& generator) {
        scalar_t dphi = 2 * Pi<scalar_t>::value / num_samples;
        std::uniform_real_distribution<scalar_t> distribution(0, Pi<scalar_t>::value);
        scalar_t offset = distribution(generator);
        std::vector<Eigen::Matrix<scalar_t, dim, 1>,
                Eigen::aligned_allocator<Eigen::Matrix<scalar_t, dim, 1>>> result;
        for (int i = 0; i < num_samples / 2; ++i) {
            scalar_t phi = i * dphi + offset;
            if (phi > Pi<scalar_t>::value) phi -= Pi<scalar_t>::value;
            int slice_n = static_cast<int>(std::ceil(num_samples * std::sin(phi)));
            if (slice_n == 0) continue;
            auto slice = SphereDiscretization<scalar_t, dim - 1>::construct(
                    radius * std::sin(phi), slice_n, generator);
            Eigen::Matrix<scalar_t, dim, 1> v;
            for (const auto& p : slice) {
                v[0] = radius * std::cos(phi);
                v.template tail<dim - 1>() = p;
                result.push_back(v);
            }
        }
        return result;
    }
    static std::vector<Eigen::Matrix<scalar_t, dim, 1>,
            Eigen::aligned_allocator<Eigen::Matrix<scalar_t, dim, 1>>> construct(
            scalar_t radius, int num_samples) {
        scalar_t dphi = 2 * Pi<scalar_t>::value / num_samples;
        std::vector<Eigen::Matrix<scalar_t, dim, 1>,
                Eigen::aligned_allocator<Eigen::Matrix<scalar_t, dim, 1>>> result;
        for (int i = 0; i < num_samples / 2; ++i) {
            scalar_t phi = i * dphi;
            if (phi > Pi<scalar_t>::value) phi -= Pi<scalar_t>::value;
            int slice_n = static_cast<int>(std::ceil(num_samples * std::sin(phi)));
            if (slice_n == 0) continue;
            auto slice = SphereDiscretization<scalar_t, dim - 1>::construct(
                    radius * std::sin(phi), slice_n);
            Eigen::Matrix<scalar_t, dim, 1> v;
            for (const auto& p : slice) {
                v[0] = radius * std::cos(phi);
                v.template tail<dim - 1>() = p;
                result.push_back(v);
            }
        }
        return result;
    }
};
 
template <typename scalar_t>
struct SphereDiscretization<scalar_t, 2> {
    template <typename generator_t>
    static std::vector<Eigen::Matrix<scalar_t, 2, 1>,
            Eigen::aligned_allocator<Eigen::Matrix<scalar_t, 2, 1>>> construct(
            scalar_t radius, int num_samples, generator_t& generator) {
        scalar_t dphi = 2 * Pi<scalar_t>::value / num_samples;
        std::uniform_real_distribution<scalar_t> distribution(0, 2 * Pi<scalar_t>::value);
        scalar_t offset = distribution(generator);
        std::vector<Eigen::Matrix<scalar_t, 2, 1>,
                Eigen::aligned_allocator<Eigen::Matrix<scalar_t, 2, 1>>> result;
        for (int i = 0; i < num_samples; ++i) {
            scalar_t phi = i * dphi + offset;
            result.emplace_back(radius * std::cos(phi), radius * std::sin(phi));
        }
        return result;
    }
    static std::vector<Eigen::Matrix<scalar_t, 2, 1>,
            Eigen::aligned_allocator<Eigen::Matrix<scalar_t, 2, 1>>> construct(
            scalar_t radius, int num_samples) {
        scalar_t dphi = 2 * Pi<scalar_t>::value / num_samples;
        std::vector<Eigen::Matrix<scalar_t, 2, 1>,
                Eigen::aligned_allocator<Eigen::Matrix<scalar_t, 2, 1>>> result;
        for (int i = 0; i < num_samples; ++i) {
            scalar_t phi = i * dphi;
            result.emplace_back(radius * std::cos(phi), radius * std::sin(phi));
        }
        return result;
    }
};
 
template <typename scalar_t>
struct SphereDiscretization<scalar_t, 1> {
    template <typename generator_t>
    static std::vector<Eigen::Matrix<scalar_t, 1, 1>,
            Eigen::aligned_allocator<Eigen::Matrix<scalar_t, 1, 1>>> construct(
            scalar_t radius, int, generator_t&) {
        return {Eigen::Matrix<scalar_t, 1, 1>(-radius), Eigen::Matrix<scalar_t, 1, 1>(radius)};
    }
    static std::vector<Eigen::Matrix<scalar_t, 1, 1>,
            Eigen::aligned_allocator<Eigen::Matrix<scalar_t, 1, 1>>> construct(
            scalar_t radius, int) {
        return {Eigen::Matrix<scalar_t, 1, 1>(-radius), Eigen::Matrix<scalar_t, 1, 1>(radius)};
    }
};
 
}  // namespace discretization_helpers
}  // namespace mm
 
#endif  // MEDUSA_BITS_DOMAINS_DISCRETIZATION_HELPERS_HPP_