libcruft-util/geom/ellipse.hpp

174 lines
4.8 KiB
C++
Raw Normal View History

2015-10-14 15:32:53 +11:00
/*
* 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>
2015-10-14 15:32:53 +11:00
*/
#ifndef __UTIL_GEOM_ELLIPSE_HPP
#define __UTIL_GEOM_ELLIPSE_HPP
2018-04-16 15:57:35 +10:00
#include "fwd.hpp"
#include "../view.hpp"
2015-10-14 15:32:53 +11:00
#include "../point.hpp"
#include "../vector.hpp"
#include <cstdlib>
#include <random>
#include <iosfwd>
2017-01-05 15:06:49 +11:00
namespace util::geom {
2015-10-14 15:32:53 +11:00
///////////////////////////////////////////////////////////////////////////
template <size_t S, typename ValueT>
2015-10-14 15:32:53 +11:00
struct ellipse {
2018-04-16 15:57:35 +10:00
// the centre point of the ellipsoid
util::point<S,ValueT> origin;
2018-04-16 15:57:35 +10:00
// the distance from the centre along each axis to the shape's edge
util::vector<S,ValueT> radius;
2018-04-16 15:57:35 +10:00
// the orientation of up for the shape
util::vector<S,ValueT> up;
2015-10-14 15:32:53 +11:00
};
2018-04-16 15:57:35 +10:00
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);
}
2018-04-16 15:57:35 +10:00
template <size_t DimensionV, typename ValueT>
point<DimensionV,ValueT>
project (
ray<DimensionV,ValueT>,
ellipse<DimensionV,ValueT>
);
2018-04-16 15:57:35 +10:00
/// 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 (util::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>
util::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 + util::vector3f {x,y,z} / d;
}
template <size_t S, typename T>
std::ostream&
operator<< (std::ostream&, ellipse<S,T>);
2017-01-05 15:06:49 +11:00
}
2015-10-14 15:32:53 +11:00
///////////////////////////////////////////////////////////////////////////////
#include "sample.hpp"
#include <cmath>
#include <random>
namespace util::geom {
template <typename T, template <size_t,typename> class K, typename G>
struct sampler<2,T,K,G>
{
static util::point<2,T> fn (K<2,T> k, G &g)
{
std::uniform_real_distribution<T> dist;
float phi = dist (g) * 2 * pi<T>;
float rho = std::sqrt (dist (g));
return util::point<2,T> {
std::cos (phi),
std::sin (phi)
} * rho * k.radius + k.origin.template as<util::vector> ();
}
};
}
2015-10-14 15:32:53 +11:00
#endif