ellipsoid.hpp

 
#ifndef UFO_GEOMETRY_ELLIPSOID_HPP
#define UFO_GEOMETRY_ELLIPSOID_HPP
 
// UFO
#include <ufo/geometry/aabb.hpp>
#include <ufo/numeric/vec.hpp>
 
// STL
#include <concepts>
#include <cstddef>
#include <format>
#include <numbers>
#include <ostream>
 
namespace ufo
{
template <std::size_t Dim = 3, std::floating_point T = float>
struct Ellipsoid {
    using value_type = T;
 
    Vec<Dim, T> center_;
 
    Vec<Dim, T> radii;
 
    constexpr Ellipsoid() noexcept = default;
 
    constexpr Ellipsoid(Vec<Dim, T> center, Vec<Dim, T> radii) noexcept
        : center_(center), radii(radii)
    {
    }
 
    constexpr Ellipsoid(Ellipsoid const&) noexcept = default;
 
    template <std::convertible_to<T> U>
    constexpr explicit Ellipsoid(Ellipsoid<Dim, U> const& other) noexcept
        : center_(Vec<Dim, T>(other.center())), radii(Vec<Dim, T>(other.radii))
    {
    }
 
    [[nodiscard]] constexpr Vec<Dim, T>& center() noexcept { return center_; }
 
    [[nodiscard]] constexpr Vec<Dim, T> const& center() const noexcept { return center_; }
 
    [[nodiscard]] static constexpr std::size_t dimension() noexcept { return Dim; }
 
    [[nodiscard]] constexpr Vec<Dim, T> min() const noexcept { return center_ - radii; }
 
    [[nodiscard]] constexpr Vec<Dim, T> max() const noexcept { return center_ + radii; }
 
    [[nodiscard]] constexpr AABB<Dim, T> aabb() const noexcept
    {
        return AABB<Dim, T>(center_, radii);
    }
 
    [[nodiscard]] constexpr bool isDegenerate() const noexcept
    {
        return any(lessThan(radii, Vec<Dim, T>()));
    }
 
    [[nodiscard]] constexpr T diameter() const noexcept { return T(2) * ufo::max(radii); }
 
    [[nodiscard]] constexpr T volume() const noexcept
    {
        if constexpr (1 == Dim) {
            return radii[0] * T(2);
        } else if constexpr (2 == Dim) {
            return std::numbers::pi_v<T> * radii[0] * radii[1];
        } else if constexpr (3 == Dim) {
            return (T(4) / T(3)) * std::numbers::pi_v<T> * radii[0] * radii[1] * radii[2];
        } else if constexpr (4 == Dim) {
            return (std::numbers::pi_v<T> * std::numbers::pi_v<T> * radii[0] * radii[1] *
                    radii[2] * radii[3]) /
                   T(2);
        } else {
            // For higher dimensions, it would need a more general formula (Gamma function)
            return T(0);
        }
    }
 
    [[nodiscard]] bool operator==(Ellipsoid const&) const = default;
};
 
//
// Deduction guides
//
 
template <std::size_t Dim, std::floating_point T>
Ellipsoid(Vec<Dim, T>, Vec<Dim, T>) -> Ellipsoid<Dim, T>;
 
template <std::size_t Dim, std::floating_point T>
std::ostream& operator<<(std::ostream& out, Ellipsoid<Dim, T> const& ellipsoid)
{
    return out << "Center: [" << ellipsoid.center() << "], radii: [" << ellipsoid.radii
               << "]";
}
 
template <std::floating_point T>
using Ellipsoid1 = Ellipsoid<1, T>;
template <std::floating_point T>
using Ellipsoid2 = Ellipsoid<2, T>;
template <std::floating_point T>
using Ellipsoid3 = Ellipsoid<3, T>;
template <std::floating_point T>
using Ellipsoid4 = Ellipsoid<4, T>;
 
using Ellipsoid1f = Ellipsoid<1, float>;
using Ellipsoid2f = Ellipsoid<2, float>;
using Ellipsoid3f = Ellipsoid<3, float>;
using Ellipsoid4f = Ellipsoid<4, float>;
 
using Ellipsoid1d = Ellipsoid<1, double>;
using Ellipsoid2d = Ellipsoid<2, double>;
using Ellipsoid3d = Ellipsoid<3, double>;
using Ellipsoid4d = Ellipsoid<4, double>;
 
}  // namespace ufo
 
template <std::size_t Dim, std::floating_point T>
  requires std::formattable<T, char>
struct std::formatter<ufo::Ellipsoid<Dim, T>> {
    constexpr auto parse(std::format_parse_context& ctx) { return ctx.begin(); }
 
    auto format(ufo::Ellipsoid<Dim, T> const& ellipsoid, std::format_context& ctx) const
    {
        return std::format_to(ctx.out(), "Center: [{}], radii: [{}]", ellipsoid.center(),
                              ellipsoid.radii);
    }
};
 
#endif  // UFO_GEOMETRY_ELLIPSOID_HPP