surf/detail/QuadricSphereConverter.hh

0001 //------------------------------- -*- C++ -*- -------------------------------//
0002 // Copyright Celeritas contributors: see top-level COPYRIGHT file for details
0003 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0004 //---------------------------------------------------------------------------//
0005 //! \file orange/surf/detail/QuadricSphereConverter.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <optional>
0010
0011 #include "corecel/math/SoftEqual.hh"
0012
0013 #include "../Plane.hh"
0014 #include "../SimpleQuadric.hh"
0015
0016 namespace celeritas
0017 {
0018 namespace detail
0019 {
0020 //---------------------------------------------------------------------------//
0021 /*!
0022  * Try to convert a simple quadric to a sphere.
0023  *
0024  * The simple quadric must have already undergone normalization (zero or one
0025  * negative signs in the second order terms).
0026  */
0027 class QuadricSphereConverter
0028 {
0029   public:
0030     // Construct with tolerance
0031     inline QuadricSphereConverter(real_type tol);
0032
0033     // Try converting to a sphere
0034     std::optional<Sphere> operator()(SimpleQuadric const& sq) const;
0035
0036   private:
0037     SoftEqual<> soft_equal_;
0038 };
0039
0040 //---------------------------------------------------------------------------//
0041 // INLINE DEFINITIONS
0042 //---------------------------------------------------------------------------//
0043 /*!
0044  * Construct with tolerance.
0045  */
0046 QuadricSphereConverter::QuadricSphereConverter(real_type tol)
0047     : soft_equal_{tol}
0048 {
0049 }
0050
0051 //---------------------------------------------------------------------------//
0052 /*!
0053  * Try converting to a sphere with this orientation.
0054  *
0055  * \verbatim
0056    x^2 + y^2 + z^2
0057        - 2x_0 - 2y_0 - 2z_0
0058        + x_0^2 + y_0^2 + z_0^2 - r^2 = 0
0059  * \endverbatim
0060  * SQ:
0061  * \verbatim
0062    a x^2 + a y^2 + a z^2 + e x + f y + g z  + h = 0
0063  * \endverbatim
0064  * so
0065  * \verbatim
0066     -2 x_0 = e/b --> x_0 = e / (-2 * a)
0067     -2 y_0 = f/b --> y_0 = f / (-2 * a)
0068     -2 z_0 = g/b --> z_0 = g / (-2 * a)
0069    x_0^2 + y_0^2 + z_0^2 - r^2 = h/b --> r^2 = (origin . origin - h/b)
0070  * \endverbatim
0071  *
0072  * All the coefficients should be positive or nearly so.
0073  */
0074 std::optional<Sphere>
0075 QuadricSphereConverter::operator()(SimpleQuadric const& sq) const
0076 {
0077     CELER_EXPECT(std::all_of(
0078         sq.second().begin(), sq.second().end(), [this](real_type v) {
0079             return v >= 0 || soft_equal_(v, 0);
0080         }));
0081     constexpr auto X = to_int(Axis::x);
0082     constexpr auto Y = to_int(Axis::y);
0083     constexpr auto Z = to_int(Axis::z);
0084
0085     auto second = sq.second();
0086     if (!soft_equal_(second[X], second[Y])
0087         || !soft_equal_(second[X], second[Z]))
0088     {
0089         // Coefficients aren't equal: it's some sort of ellipsoid
0090         return {};
0091     }
0092
0093     real_type const inv_norm = 3 / (second[0] + second[1] + second[2]);
0094     CELER_ASSERT(inv_norm > 0);
0095
0096     Real3 origin = make_array(sq.first());
0097     origin *= real_type{-0.5} * inv_norm;
0098
0099     real_type radius_sq = dot_product(origin, origin) - sq.zeroth() * inv_norm;
0100     if (radius_sq <= 0)
0101     {
0102         // No real solution
0103         return {};
0104     }
0105
0106     // Clear potential signed zeros before returning
0107     origin += real_type{0};
0108     return Sphere::from_radius_sq(origin, radius_sq);
0109 }
0110
0111 //---------------------------------------------------------------------------//
0112 }  // namespace detail
0113 }  // namespace celeritas