libcruft-util/geom/ellipse.hpp

198 lines
5.7 KiB
C++

/*
* 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 <cstdlib>
#include <random>
#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 = std::normal_distribution<float> (0, a2) (generator);
auto const y = std::normal_distribution<float> (0, b2) (generator);
auto const z = std::normal_distribution<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.hpp"
#include <cmath>
#include <random>
namespace cruft::geom {
/// Specialisation for uniform sampling of ellipses
template <typename T>
struct sampler<ellipse<2,T>>
{
/// Generate a random point within the ellipse.
template <typename GeneratorT>
static auto
eval (ellipse<2,T> shape, GeneratorT &&g)
{
// 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 * shape.radius;
return shape.origin + offset.template as<cruft::vector> ();
}
};
#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
}