orangeinp/detail/PolygonUtils.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/orangeinp/detail/PolygonUtils.hh
0006 //! \brief Utility standalone functions for polygons in 2D or 3D space.
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include <algorithm>
0011 #include <iostream>
0012 #include <vector>
0013
0014 #include "corecel/cont/Array.hh"
0015 #include "corecel/cont/Range.hh"
0016 #include "corecel/cont/Span.hh"
0017 #include "corecel/math/ArrayOperators.hh"
0018 #include "corecel/math/ArrayUtils.hh"
0019 #include "corecel/math/SoftEqual.hh"
0020 #include "orange/OrangeTypes.hh"
0021
0022 namespace celeritas
0023 {
0024 namespace orangeinp
0025 {
0026 namespace detail
0027 {
0028 //---------------------------------------------------------------------------//
0029 /*!
0030  *  Polygon orientation based on ordering of vertices.
0031  */
0032 enum class Orientation
0033 {
0034     clockwise = -1,
0035     collinear = 0,
0036     counterclockwise = 1,
0037 };
0038
0039 //---------------------------------------------------------------------------//
0040 // FREE FUNCTIONS
0041 //---------------------------------------------------------------------------//
0042 /*!
0043  * Find orientation of ordered vertices in 2D coordinates.
0044  */
0045 template<class T>
0046 inline Orientation calc_orientation(celeritas::Array<T, 2> const& a,
0047                                     celeritas::Array<T, 2> const& b,
0048                                     celeritas::Array<T, 2> const& c)
0049 {
0050     auto crossp = (b[0] - a[0]) * (c[1] - b[1]) - (b[1] - a[1]) * (c[0] - b[0]);
0051     return crossp < 0   ? Orientation::clockwise
0052            : crossp > 0 ? Orientation::counterclockwise
0053                         : Orientation::collinear;
0054 }
0055
0056 //---------------------------------------------------------------------------//
0057 /*!
0058  * Test whether a 2D polygon has the given orientation.
0059  *
0060  * The list of input corners must have at least 3 points to be a polygon.
0061  */
0062 inline bool has_orientation(Span<Real2 const> corners, Orientation o)
0063 {
0064     CELER_EXPECT(corners.size() > 2);
0065     for (auto i : range(corners.size()))
0066     {
0067         auto j = (i + 1) % corners.size();
0068         auto k = (i + 2) % corners.size();
0069         if (calc_orientation<Real2::value_type>(
0070                 corners[i], corners[j], corners[k])
0071             != o)
0072             return false;
0073     }
0074     return true;
0075 }
0076
0077 //---------------------------------------------------------------------------//
0078 /*!
0079  * Whether the orientation is the same or degenerate if allowed.
0080  */
0081 inline bool
0082 is_same_orientation(Orientation a, Orientation b, bool degen_ok = false)
0083 {
0084     if (a == Orientation::collinear || b == Orientation::collinear)
0085     {
0086         return degen_ok;
0087     }
0088     return (a == b);
0089 }
0090
0091 //---------------------------------------------------------------------------//
0092 /*!
0093  * Functor for calculating orientation with a tolerance for collinearity.
0094  *
0095  * Collinearity is based on a supplied absolute tolerance. For three ordered
0096  * points a, b, c, point b is collinear if the displacement, d, is less than
0097  * the absolute tolerance.
0098  * \verbatim
0099                 b
0100                . .
0101              .  .  .
0102            .    .    .
0103          .      . d    .
0104        .  t     .        .
0105      a . . . . . . . . . . c
0106    \endverbatim
0107  * The displacement is calculated as follows.
0108  *
0109  * Let:
0110  * \verbatim
0111    u = b - a
0112    v = c - a
0113    \endverbatim
0114  *
0115  * In 2D, the cross product can be written as,
0116  *
0117  * \verbatim
0118    u x v = |u| |v| sin(t),
0119    \endverbatim
0120  *
0121  * noting that this is a different cross product (different vectors) compared
0122  * to the cross product used for orientation determination. Geometrically, the
0123  * displacement can be calculated as,
0124  *
0125  *  \verbatim
0126     d = |u| sin(t).
0127     \endverbatim
0128  *
0129  * Therefore,
0130  *
0131  * \verbatim
0132    d = |u| (u x v) / (|u| |v|)
0133      = (u x v)/|v|.
0134    \endverbatim
0135  */
0136 template<class T>
0137 class SoftOrientation
0138 {
0139   public:
0140     using real_type = T;
0141     using Real2 = Array<real_type, 2>;
0142
0143     // Construct with default tolerance
0144     CELER_CONSTEXPR_FUNCTION SoftOrientation() : soft_zero_{} {}
0145
0146     // Construct with specified absolute tolerance
0147     explicit CELER_FUNCTION SoftOrientation(real_type abs_tol)
0148         : soft_zero_{abs_tol}
0149     {
0150     }
0151
0152     // Calculate orientation with tolerance for collinearity
0153     CELER_FUNCTION Orientation operator()(Real2 const& a,
0154                                           Real2 const& b,
0155                                           Real2 const& c) const
0156     {
0157         Real2 u{b[0] - a[0], b[1] - a[1]};
0158         Real2 v{c[0] - a[0], c[1] - a[1]};
0159
0160         auto cross_product = (u[0] * v[1] - v[0] * u[1]);
0161         auto magnitude = std::hypot(v[0], v[1]);
0162
0163         if (magnitude == 0 || soft_zero_(cross_product / magnitude))
0164         {
0165             return Orientation::collinear;
0166         }
0167         else
0168         {
0169             return calc_orientation<real_type>(a, b, c);
0170         }
0171     }
0172
0173   private:
0174     SoftZero<real_type> soft_zero_;
0175 };
0176
0177 //---------------------------------------------------------------------------//
0178 /*!
0179  * Check if a 2D polygon is convex.
0180  *
0181  * \param corners the vertices of the polygon
0182  * \param degen_ok allow consecutive collinear points
0183  */
0184 inline bool is_convex(Span<Real2 const> corners, bool degen_ok = false)
0185 {
0186     CELER_EXPECT(corners.size() > 2);
0187     auto ref = Orientation::collinear;
0188     for (auto i : range<size_type>(corners.size()))
0189     {
0190         auto j = (i + 1) % corners.size();
0191         auto k = (i + 2) % corners.size();
0192         auto cur = calc_orientation<Real2::value_type>(
0193             corners[i], corners[j], corners[k]);
0194         if (ref == Orientation::collinear)
0195         {
0196             // First non-collinear point
0197             ref = cur;
0198         }
0199         if (!is_same_orientation(cur, ref, degen_ok))
0200         {
0201             // Prohibited collinear orientation, or different orientation from
0202             // reference
0203             return false;
0204         }
0205     }
0206     return true;
0207 }
0208
0209 //---------------------------------------------------------------------------//
0210 /*!
0211  * Return the non-collinear subset of the supplied corners, as a copy.
0212  *
0213  * Points are checked for collinearity dynamically, i.e, if a point is found to
0214  * be collinear, it is not used for future collinearity checks.
0215  */
0216 inline std::vector<Real2>
0217 filter_collinear_points(std::vector<Real2> const& corners, double abs_tol)
0218 {
0219     CELER_EXPECT(corners.size() >= 3);
0220
0221     std::vector<Real2> result;
0222     result.reserve(corners.size());
0223
0224     SoftOrientation<Real2::value_type> soft_ori(abs_tol);
0225
0226     // Temporarily assume first point is not collinear. This is necessary to
0227     // handle the case where all points are locally collinear, but some points
0228     // are globally non-collinear, e.g., a many-sided regular polygon.
0229     result.push_back(corners[0]);
0230
0231     for (auto i : range<size_type>(1, corners.size()))
0232     {
0233         auto j = i + 1 < corners.size() ? i + 1 : 0;
0234
0235         if (soft_ori(result.back(), corners[i], corners[j])
0236             != Orientation::collinear)
0237         {
0238             result.push_back(corners[i]);
0239         }
0240     }
0241
0242     // Make sure there are enough filtered points to specify a polygon.
0243     CELER_ASSERT(result.size() >= 3);
0244
0245     // If it turns out that the first point is actually collinear, remove it.
0246     if (soft_ori(result.back(), result[0], result[1]) == Orientation::collinear)
0247     {
0248         result.erase(result.begin());
0249     }
0250
0251     // It shouldn't be possible for the potential removal of the first point
0252     // to leave us with fewer than 3 points, but check just in case.
0253     CELER_ENSURE(result.size() >= 3);
0254
0255     return result;
0256 }
0257
0258 //---------------------------------------------------------------------------//
0259 /*!
0260  * Calculate the min/max values of a polygon for a given dimension.
0261  */
0262 template<class T>
0263 inline std::pair<T, T>
0264 find_extrema(Span<Array<T, 2> const> polygon, size_type dim)
0265 {
0266     CELER_EXPECT(polygon.size() >= 3);
0267     CELER_EXPECT(dim < 2);
0268
0269     auto [poly_min_it, poly_max_it] = std::minmax_element(
0270         polygon.begin(), polygon.end(), [dim](auto const& a, auto const& b) {
0271             return a[dim] < b[dim];
0272         });
0273     return {(*poly_min_it)[dim], (*poly_max_it)[dim]};
0274 }
0275
0276 template<class T>
0277 inline auto find_extrema(T const& polygon, size_type dim)
0278 {
0279     return find_extrema(make_span(polygon), dim);
0280 }
0281
0282 //! Calculate the min/max values of a polygon for an x/y dimension
0283 template<class T>
0284 inline auto find_extrema(T const& polygon, Axis ax)
0285 {
0286     CELER_EXPECT(ax < Axis::z);
0287     return find_extrema(polygon, to_int(ax));
0288 }
0289
0290 //---------------------------------------------------------------------------//
0291 /*!
0292  * Construct a normal from three counterclockwise points.
0293  *
0294  * The direction of the normal is dictated by the right-hand rule, assuming
0295  * the normal vector, C, is given by:
0296  * \f[
0297  * \vec{A} = \vec{p_1} - \vec{p_0},
0298  * \vec{B} = \vec{p_2} - \vec{p_0},
0299  * \vec{C} = \vec{A} \times  \vec{B}.
0300  * \f]
0301  * \verbatim
0302              ^
0303              | C
0304              |
0305              |
0306    p1 _______| p0
0307         A   /
0308            / B
0309           /
0310          p2
0311    \endverbatim
0312  *
0313  * Passing in collinear or coincident points will result in an invalid normal.
0314  */
0315 inline Real3
0316 normal_from_triangle(Real3 const& p0, Real3 const& p1, Real3 const& p2)
0317 {
0318     return make_unit_vector(cross_product(p1 - p0, p2 - p0));
0319 }
0320
0321 //---------------------------------------------------------------------------//
0322 }  // namespace detail
0323 }  // namespace orangeinp
0324 }  // namespace celeritas