libcruft-util/geom/sample/surface.hpp

/*
 * This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this
 * file, You can obtain one at http://mozilla.org/MPL/2.0/.
 *
 * Copyright 2015-2018 Danny Robson <danny@nerdcruft.net>
 */

#pragma once

#include "volume.hpp"

#include "../ops.hpp"

#include "../../coord/fwd.hpp"
#include "../../extent.hpp"
#include "../../region.hpp"
#include "../../random.hpp"

#include <cstddef>
#include <random>


///////////////////////////////////////////////////////////////////////////////
namespace cruft::geom::sample {
    /// Approximate a poisson disc sampling through the "Mitchell's Best
    /// Candidate" algorithm.
    ///
    /// We try to keep adding a new point to an output set.
    /// * Generate `candidates_count` new points
    /// * Remove any point that fails AcceptT
    /// * Order the points by minimum distance to the output set
    /// * Remove any point with a distance less than DistanceT
    /// * Select the point with the highest distance
    ///
    /// Keep iterating until no point is added to the output set.
    ///
    /// \tparam SamplerT A surface sampler
    /// \tparam DistanceT The type of point-to-point distances
    /// \tparam GeneratorT The random generator passed to the sampler
    /// \tparam AcceptT A unary bool function that returns true if a sampled
    ///     point is permissible. The point is counted in the candidate set
    ///     regardless.
    ///
    /// \return A vector of the computed points
    template <
        typename SamplerT,
        typename DistanceT,
        typename GeneratorT,
        typename AcceptT
    >
    auto
    poisson (
        SamplerT const &sampler,
        GeneratorT &&gen,
        AcceptT &&accept,
        DistanceT &&minimum_distance,
        std::size_t const candidates_count
    ) {
        using point_type = decltype (sampler.eval (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.
        for (std::size_t i = 0; i < candidates_count; ++i) {
            auto const p = sampler.eval (gen);
            if (!accept (p))
                continue;

            selected.push_back (p);
            break;
        }

        // We couldn't find an initial candidate so abort.
        if (selected.empty ())
            return selected;

        // 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 sampler.eval (gen);
                }
            );

            // Remove points that aren't acceptable
            std::erase_if (
                candidates,
                [&] (auto const &p) { return !accept (p); }
            );

            // 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 = cruft::min (d2, cruft::distance2 (p, q));

                // record if it's the furthest away
                if (d2 > best_distance2 && d2 > cruft::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;
    }


    /// A convenience function that forwards to poisson with unconditional
    /// point acceptance.
    template <
        typename SamplerT,
        typename DistanceT,
        typename GeneratorT
    >
    auto
    poisson (
        SamplerT &&sampler,
        GeneratorT &&gen,
        DistanceT &&minimum_distance,
        std::size_t candidates_count
    ) {
        return poisson (
            std::forward<SamplerT> (sampler),
            std::forward<GeneratorT> (gen),
            [] (auto&&...) { return true; },
            std::forward<DistanceT> (minimum_distance),
            candidates_count
        );
    }


    /// A surface sampler specialisation for 2d extents.
    ///
    /// The actual work is handed off to the volume sampler, as it's
    /// equivalent in 2 dimensions.
    template <typename T>
    class surface<extent<2,T>> :
        public volume<extent<2,T>>
    { using volume<extent<2,T>>::volume; };


    /// A surface sampler specialisation for 2d regions
    ///
    /// We encapsulate an extent sampler and an offset, then mainly pass off
    /// the work to the extent sampler.
    template <typename T>
    class surface<region<2,T>> {
    public:
        using shape_type = region<2,T>;

        surface (shape_type const &_shape)
            : m_extent (_shape.e)
            , m_offset (_shape.p.template as<vector> ())
        { ; }


        template <typename GeneratorT>
        decltype(auto)
        eval (GeneratorT &&gen)
        {
            return m_extent.eval (
                std::forward<GeneratorT> (gen)
            ) + m_offset;
        }


    private:
        surface<extent<2,T>> m_extent;
        cruft::vector<2,T> m_offset;
    };
}