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  * Find the orientation of a simple 2D polygon from its signed area.
0059  *
0060  * The input must have at least three ordered vertices and may be nonconvex.
0061  * The closing edge is included implicitly. A zero signed area returns
0062  * \c Orientation::collinear.
0063  *
0064  * Sum signed triangle areas relative to the first vertex to avoid cancellation
0065  * from large coordinate offsets. The first and closing edges contribute zero.
0066  */
0067 inline Orientation calc_orientation(Span<Real2 const> corners)
0068 {
0069     CELER_EXPECT(corners.size() >= 3);
0070
0071     auto const& origin = corners.front();
0072     real_type twice_area = 0;
0073     for (auto i : range<size_type>(1, corners.size() - 1))
0074     {
0075         auto const a = corners[i] - origin;
0076         auto const b = corners[i + 1] - origin;
0077         twice_area += a[0] * b[1] - b[0] * a[1];
0078     }
0079     return twice_area < 0   ? Orientation::clockwise
0080            : twice_area > 0 ? Orientation::counterclockwise
0081                             : Orientation::collinear;
0082 }
0083
0084 //---------------------------------------------------------------------------//
0085 /*!
0086  * Test whether a 2D polygon has the given orientation.
0087  *
0088  * The list of input corners must have at least 3 points to be a polygon.
0089  */
0090 inline bool has_orientation(Span<Real2 const> corners, Orientation o)
0091 {
0092     CELER_EXPECT(corners.size() > 2);
0093     for (auto i : range(corners.size()))
0094     {
0095         auto j = (i + 1) % corners.size();
0096         auto k = (i + 2) % corners.size();
0097         if (calc_orientation<Real2::value_type>(
0098                 corners[i], corners[j], corners[k])
0099             != o)
0100             return false;
0101     }
0102     return true;
0103 }
0104
0105 //---------------------------------------------------------------------------//
0106 /*!
0107  * Whether the orientation is the same or degenerate if allowed.
0108  */
0109 inline bool is_same_orientation(
0110     Orientation a, Orientation b, bool degen_ok = false)
0111 {
0112     if (a == Orientation::collinear || b == Orientation::collinear)
0113     {
0114         return degen_ok;
0115     }
0116     return (a == b);
0117 }
0118
0119 //---------------------------------------------------------------------------//
0120 /*!
0121  * Functor for calculating orientation with a tolerance for collinearity.
0122  *
0123  * Collinearity is based on a supplied absolute tolerance. For three ordered
0124  * points a, b, c, point b is collinear if the displacement, d, is less than
0125  * the absolute tolerance.
0126  * \verbatim
0127                 b
0128                . .
0129              .  .  .
0130            .    .    .
0131          .      . d    .
0132        .  t     .        .
0133      a . . . . . . . . . . c
0134    \endverbatim
0135  * The displacement is calculated as follows.
0136  *
0137  * Let:
0138  * \verbatim
0139    u = b - a
0140    v = c - a
0141    \endverbatim
0142  *
0143  * In 2D, the cross product can be written as,
0144  *
0145  * \verbatim
0146    u x v = |u| |v| sin(t),
0147    \endverbatim
0148  *
0149  * noting that this is a different cross product (different vectors) compared
0150  * to the cross product used for orientation determination. Geometrically, the
0151  * displacement can be calculated as,
0152  *
0153  *  \verbatim
0154     d = |u| sin(t).
0155     \endverbatim
0156  *
0157  * Therefore,
0158  *
0159  * \verbatim
0160    d = |u| (u x v) / (|u| |v|)
0161      = (u x v)/|v|.
0162    \endverbatim
0163  */
0164 template<class T>
0165 class SoftOrientation
0166 {
0167   public:
0168     using real_type = T;
0169     using Real2 = Array<real_type, 2>;
0170
0171     // Construct with default tolerance
0172     CELER_CONSTEXPR_FUNCTION SoftOrientation() : soft_zero_{} {}
0173
0174     // Construct with specified absolute tolerance
0175     explicit CELER_FUNCTION SoftOrientation(real_type abs_tol)
0176         : soft_zero_{abs_tol}
0177     {
0178     }
0179
0180     // Calculate orientation with tolerance for collinearity
0181     CELER_FUNCTION Orientation operator()(
0182         Real2 const& a, Real2 const& b, Real2 const& c) const
0183     {
0184         Real2 u{b[0] - a[0], b[1] - a[1]};
0185         Real2 v{c[0] - a[0], c[1] - a[1]};
0186
0187         auto cross_product = (u[0] * v[1] - v[0] * u[1]);
0188         auto magnitude = std::hypot(v[0], v[1]);
0189
0190         if (magnitude == 0 || soft_zero_(cross_product / magnitude))
0191         {
0192             return Orientation::collinear;
0193         }
0194         else
0195         {
0196             return calc_orientation<real_type>(a, b, c);
0197         }
0198     }
0199
0200   private:
0201     SoftZero<real_type> soft_zero_;
0202 };
0203
0204 //---------------------------------------------------------------------------//
0205 /*!
0206  * Check if a 2D polygon is convex.
0207  *
0208  * \param corners the vertices of the polygon
0209  * \param degen_ok allow consecutive collinear points
0210  */
0211 inline bool is_convex(Span<Real2 const> corners, bool degen_ok = false)
0212 {
0213     CELER_EXPECT(corners.size() > 2);
0214     auto ref = Orientation::collinear;
0215     for (auto i : range<size_type>(corners.size()))
0216     {
0217         auto j = (i + 1) % corners.size();
0218         auto k = (i + 2) % corners.size();
0219         auto cur = calc_orientation<Real2::value_type>(
0220             corners[i], corners[j], corners[k]);
0221         if (ref == Orientation::collinear)
0222         {
0223             // First non-collinear point
0224             ref = cur;
0225         }
0226         if (!is_same_orientation(cur, ref, degen_ok))
0227         {
0228             // Prohibited collinear orientation, or different orientation from
0229             // reference
0230             return false;
0231         }
0232     }
0233     return true;
0234 }
0235
0236 //---------------------------------------------------------------------------//
0237 /*!
0238  * Return the non-collinear subset of the supplied corners, as a copy.
0239  *
0240  * Points are checked for collinearity dynamically, i.e, if a point is found to
0241  * be collinear, it is not used for future collinearity checks.
0242  */
0243 inline std::vector<Real2> filter_collinear_points(
0244     std::vector<Real2> const& corners, double abs_tol)
0245 {
0246     CELER_EXPECT(corners.size() >= 3);
0247
0248     std::vector<Real2> result;
0249     result.reserve(corners.size());
0250
0251     SoftOrientation<Real2::value_type> soft_ori(abs_tol);
0252
0253     // Temporarily assume first point is not collinear. This is necessary to
0254     // handle the case where all points are locally collinear, but some points
0255     // are globally non-collinear, e.g., a many-sided regular polygon.
0256     result.push_back(corners[0]);
0257
0258     for (auto i : range<size_type>(1, corners.size()))
0259     {
0260         auto j = i + 1 < corners.size() ? i + 1 : 0;
0261
0262         if (soft_ori(result.back(), corners[i], corners[j])
0263             != Orientation::collinear)
0264         {
0265             result.push_back(corners[i]);
0266         }
0267     }
0268
0269     // Make sure there are enough filtered points to specify a polygon.
0270     CELER_ASSERT(result.size() >= 3);
0271
0272     // If it turns out that the first point is actually collinear, remove it.
0273     if (soft_ori(result.back(), result[0], result[1]) == Orientation::collinear)
0274     {
0275         result.erase(result.begin());
0276     }
0277
0278     // It shouldn't be possible for the potential removal of the first point
0279     // to leave us with fewer than 3 points, but check just in case.
0280     CELER_ENSURE(result.size() >= 3);
0281
0282     return result;
0283 }
0284
0285 //---------------------------------------------------------------------------//
0286 /*!
0287  * Calculate the min/max values of a polygon for a given dimension.
0288  */
0289 template<class T>
0290 inline std::pair<T, T> find_extrema(Span<Array<T, 2> const> polygon,
0291                                     size_type dim)
0292 {
0293     CELER_EXPECT(polygon.size() >= 3);
0294     CELER_EXPECT(dim < 2);
0295
0296     auto [poly_min_it, poly_max_it] = std::minmax_element(
0297         polygon.begin(), polygon.end(), [dim](auto const& a, auto const& b) {
0298             return a[dim] < b[dim];
0299         });
0300     return {(*poly_min_it)[dim], (*poly_max_it)[dim]};
0301 }
0302
0303 template<class T>
0304 inline auto find_extrema(T const& polygon, size_type dim)
0305 {
0306     return find_extrema(make_span(polygon), dim);
0307 }
0308
0309 //! Calculate the min/max values of a polygon for an x/y dimension
0310 template<class T>
0311 inline auto find_extrema(T const& polygon, Axis ax)
0312 {
0313     CELER_EXPECT(ax < Axis::z);
0314     return find_extrema(polygon, to_int(ax));
0315 }
0316
0317 //---------------------------------------------------------------------------//
0318 /*!
0319  * Construct a normal from three counterclockwise points.
0320  *
0321  * The direction of the normal is dictated by the right-hand rule, assuming
0322  * the normal vector, C, is given by:
0323  * \f[
0324  * \vec{A} = \vec{p_1} - \vec{p_0},
0325  * \vec{B} = \vec{p_2} - \vec{p_0},
0326  * \vec{C} = \vec{A} \times  \vec{B}.
0327  * \f]
0328  * \verbatim
0329              ^
0330              | C
0331              |
0332              |
0333    p1 _______| p0
0334         A   /
0335            / B
0336           /
0337          p2
0338    \endverbatim
0339  *
0340  * Passing in collinear or coincident points will result in an invalid normal.
0341  */
0342 inline Real3 normal_from_triangle(
0343     Real3 const& p0, Real3 const& p1, Real3 const& p2)
0344 {
0345     return make_unit_vector(cross_product(p1 - p0, p2 - p0));
0346 }
0347
0348 //---------------------------------------------------------------------------//
0349 }  // namespace detail
0350 }  // namespace orangeinp
0351 }  // namespace celeritas