libcruft-util/cruft/util/geom/ellipse.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 "fwd.hpp"

#include "../view.hpp"
#include "../point.hpp"
#include "../vector.hpp"

#include <cruft/util/rand/distribution/normal.hpp>

#include <cstdlib>
#include <iosfwd>


namespace cruft::geom {
    ///////////////////////////////////////////////////////////////////////////
    template <size_t S, typename ValueT>
    struct ellipse {
        // the centre point of the ellipsoid
        cruft::point<S,ValueT> origin;

        // the distance from the centre along each axis to the shape's edge
        cruft::vector<S,ValueT> radius;

        // the orientation of up for the shape
        cruft::vector<S,ValueT> up;
    };

    using ellipse2f = ellipse<2,float>;
    using ellipse3f = ellipse<3,float>;


    template <size_t S, typename T>
    bool
    intersects (ellipse<S,T>, point<S,T>);


    /// returns the approximate surface area of the ellipse
    ///
    /// the general form involves _substantially_ more expensive and
    /// complicated maths which is prohibitive right now.
    ///
    /// the relative error should be at most 1.061%
    inline float
    area (ellipse3f self)
    {
        auto const semiprod = self.radius * self.radius.indices<1,2,0> ();
        auto const semipow  = pow (semiprod, 1.6f);

        return 4 * pi<float> * std::pow (sum (semipow) / 3, 1/1.6f);
    }


    template <typename T>
    T
    area (ellipse<2,T> self)
    {
        return pi<T> * product (self.radius);
    }


    template <size_t S, typename T>
    T
    volume (ellipse<S,T> self)
    {
        return 4 / T{3} * pi<T> * product (self.radius);
    }


    template <size_t DimensionV, typename ValueT>
    point<DimensionV,ValueT>
    project (
        ray<DimensionV,ValueT>,
        ellipse<DimensionV,ValueT>
    );


    /// returns the distance along a ray to the surface of an ellipse.
    ///
    /// returns infinity if there is no intersection
    template <size_t DimensionV, typename ValueT>
    ValueT
    distance (
        ray<DimensionV,ValueT>,
        ellipse<DimensionV,ValueT>
    );


    // generate a covering ellipsoid for an arbitrary set of points
    //
    // this isn't guaranteed to be optimal in any specific sense. but it
    // ought not be outrageously inefficient...
    ellipse3f
    cover (cruft::view<const point3f*>);


    /// returns a point that is uniformly distributed about the ellipse
    /// surface.
    ///
    /// NOTE: I don't have a strong proof that the below is in fact properly
    /// uniformly distributed, so if you need a strong guarantee for your work
    /// then it might not be the best option. But visual inspection appears to
    /// confirm there aren't obvious patterns.
    ///
    /// the concept was taken from: https://math.stackexchange.com/a/2514182
    template <typename RandomT>
    cruft::point3f
    sample_surface (ellipse3f self, RandomT &generator)
    {
        const auto [a, b, c] = self.radius;
        const auto a2 = a * a;
        const auto b2 = b * b;
        const auto c2 = c * c;

        // generate a direction vector from a normally distributed random variable
        auto const x = cruft::rand::distribution::normal<float> (0, a2) (generator);
        auto const y = cruft::rand::distribution::normal<float> (0, b2) (generator);
        auto const z = cruft::rand::distribution::normal<float> (0, c2) (generator);

        // find the distance to the surface along the direction vector
        auto const d = std::sqrt (x * x / a2 + y * y / b2 + z * z / c2);

        return self.origin + cruft::vector3f {x,y,z} / d;
    }


    template <size_t S, typename T>
    std::ostream&
    operator<< (std::ostream&, ellipse<S,T>);
}


///////////////////////////////////////////////////////////////////////////////
#include "sample/fwd.hpp"

#include "../random.hpp"

#include <cmath>
#include <random>

namespace cruft::geom::sample {
    /// Specialisation for uniform sampling of ellipses
    template <typename T>
    class volume<ellipse<2,T>> {
    public:
        using shape_type = ellipse<2,T>;

        volume (shape_type &&) = delete;
        volume (shape_type const &target):
            m_target (target)
        { ; }


        /// Generate a random point within the ellipse.
        template <typename GeneratorT>
        auto
        eval (GeneratorT &&g) const noexcept
        {
            // We use a two step process:
            //   * Generate a point within a unit sphere
            //   * Transform the point to an ellipse.

            // TODO: We assume floating point for the time being because it
            //       simplifies interaction with trig routines. There's no
            //       intrinsic reason for this limitation though.
            static_assert (
                std::is_floating_point_v<T>,
                    "The current implementation assumes floating point."
            );

            // Choose a direction and a distance within the unit circle.
            //
            // `sqrt` of the distance is used to ensure a uniform
            // distribution.
            T phi = random::uniform<T> (g) * 2 * pi<T>;
            T rho = std::sqrt (random::uniform<T> (g));

            cruft::point2<T> const circle_pos {
                std::cos (phi),
                std::sin (phi)
            };

            auto const offset = circle_pos * rho * m_target.radius;
            return m_target.origin + offset.template as<cruft::vector> ();
        }


    private:
        shape_type const &m_target;
    };


#if 0
    // TODO: We should implement a higher dimensional ellipsoid sampler for
    //       efficiency gains over rejection sampling that we currently use.
    template <typename T>
    struct sampler<ellipse<3,T>>
    {
        template <typename GeneratorT>
        static cruft::point<3,T>
        eval (ellipse<3,T>, GeneratorT&&);
    };
#endif
}