libcruft-util/geom/sample.hpp

195 lines
6.3 KiB
C++

/*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
* Copyright 2015-2018 Danny Robson <danny@nerdcruft.net>
*/
#ifndef __UTIL_GEOM_SAMPLE_HPP
#define __UTIL_GEOM_SAMPLE_HPP
#include "../coord/fwd.hpp"
#include "../extent.hpp"
#include "ops.hpp"
#include <cstddef>
#include <random>
///////////////////////////////////////////////////////////////////////////////
namespace util::geom {
/// a function object that selects a uniformly random point inside a shape
/// using a provided random generator. the point will lie within the shape,
/// inclusive of boundaries.
///
/// may be specialised for arbitrary shapes but uses rejection sampling
/// as a safe default. this implies that execution may not take a constant
/// time.
///
/// \tparam S coordinate type dimension
/// \tparam T value type for the shape/coordinate
/// \tparam K the shape type to test
/// \tparam G a UniformRandomBitGenerator, eg std::min19937
template <
size_t S,
typename T,
template <size_t,typename> class K,
typename G
>
struct sampler {
static point<S,T>
fn (K<S,T> k, G &g)
{
auto b = bounds (k);
while (true) {
auto p = sample (b, g);
if (intersects (k, p))
return p;
}
}
};
///////////////////////////////////////////////////////////////////////////
/// a convenience function that calls sample::fn to select a random point
/// in a provided shape.
template <
size_t S,
typename T,
template <size_t,typename> class K,
typename G // random generator
>
point<S,T>
sample (K<S,T> k, G &g)
{
return sampler<S,T,K,G>::fn (k, g);
}
std::vector<util::point2f>
poisson_sample (util::extent2i, float distance, int samples);
namespace surface {
// a generator of samples that lie on the surface of a shape
template <typename ShapeT>
class sampler;
template <typename ShapeT>
sampler (ShapeT const&) -> sampler<std::decay_t<ShapeT>>;
//---------------------------------------------------------------------
template <size_t S, typename T>
class sampler<util::extent<S,T>> {
public:
sampler (util::extent<S,T> _target):
target (_target)
{ ; }
template <typename GeneratorT>
util::point<S,T>
operator() (GeneratorT &&gen) const noexcept
{
util::point<S,T> p;
for (size_t i = 0; i < S; ++i) {
if constexpr (std::is_floating_point_v<T>) {
p[i] = std::uniform_real_distribution<T> (0, target[i]) (gen);
} else {
CHECK_GE (target[i], T{1});
p[i] = std::uniform_int_distribution<T> (0, target[i] - 1) (gen);
}
}
return p;
}
private:
util::extent<S,T> target;
};
/// approximate a poisson disc sampling through mitchell's best candidate.
///
/// try to keep adding a new point to a set. each new point is the
/// best of a set of candidates. the 'best' is the point that is
/// furthest from all selected points.
///
/// \return a vector of the computed points
template <typename SamplerT, typename DistanceT, typename GeneratorT>
auto
poisson (SamplerT const &target,
GeneratorT &&gen,
DistanceT &&minimum_distance,
size_t candidates_count)
{
using point_type = decltype (target (gen));
using value_type = typename point_type::value_type;
std::vector<point_type> selected;
std::vector<point_type> candidates;
// prime the found elements list with an initial point we can
// perform distance calculations on
selected.push_back (target (gen));
// keep trying to add one more new point
while (1) {
// generate a group of candidate points
candidates.clear ();
std::generate_n (
std::back_inserter (candidates),
candidates_count,
[&] (void) {
return target (gen);
}
);
// find the point whose minimum distance to the existing
// points is the greatest (while also being greater than the
// required minimum distance);
auto best_distance2 = std::numeric_limits<value_type>::lowest ();
size_t best_index = 0;
for (size_t i = 0; i < candidates.size (); ++i) {
auto const p = candidates[i];
auto d2 = std::numeric_limits<value_type>::max ();
// find the minimum distance from this candidate to the
// selected points
for (auto q: selected)
d2 = util::min (d2, util::distance2 (p, q));
// record if it's the furthest away
if (d2 > best_distance2 && d2 > util::pow (minimum_distance (p), 2)) {
best_distance2 = d2;
best_index = i;
}
}
// if we didn't find a suitable point then we give up and
// return the points we found, otherwise record the best point
if (best_distance2 <= 0)
break;
selected.push_back (candidates[best_index]);
}
return selected;
}
}
}
#endif