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 <iostream>
0011 #include <vector>
0012
0013 #include "corecel/cont/Array.hh"
0014 #include "corecel/cont/Range.hh"
0015 #include "corecel/cont/Span.hh"
0016 #include "corecel/math/ArrayOperators.hh"
0017 #include "corecel/math/ArrayUtils.hh"
0018 #include "corecel/math/SoftEqual.hh"
0019 #include "orange/OrangeTypes.hh"
0020
0021 namespace celeritas
0022 {
0023 namespace orangeinp
0024 {
0025 namespace detail
0026 {
0027 //---------------------------------------------------------------------------//
0028 /*!
0029  *  Polygon orientation based on ordering of vertices.
0030  */
0031 enum class Orientation
0032 {
0033     clockwise = -1,
0034     collinear = 0,
0035     counterclockwise = 1,
0036 };
0037
0038 //---------------------------------------------------------------------------//
0039 /*!
0040  * Functor for calculating orientation with a tolerance for collinearity.
0041  */
0042 template<class T>
0043 class SoftOrientation
0044 {
0045   public:
0046     using real_type = T;
0047     using Real2 = Array<real_type, 2>;
0048
0049     // Construct with default tolerance
0050     CELER_CONSTEXPR_FUNCTION SoftOrientation() : soft_zero_{} {}
0051
0052     // Construct with specified absolute tolerance
0053     explicit CELER_FUNCTION SoftOrientation(real_type abs_tol)
0054         : soft_zero_{abs_tol}
0055     {
0056     }
0057
0058     // Calculate orientation with tolerance for collinearity
0059     CELER_FUNCTION Orientation operator()(Real2 const& a,
0060                                           Real2 const& b,
0061                                           Real2 const& c) const
0062     {
0063         auto crossp = (b[0] - a[0]) * (c[1] - b[1])
0064                       - (b[1] - a[1]) * (c[0] - b[0]);
0065         return soft_zero_(crossp) ? Orientation::collinear
0066                : crossp < 0       ? Orientation::clockwise
0067                                   : Orientation::counterclockwise;
0068     }
0069
0070   private:
0071     SoftZero<real_type> soft_zero_;
0072 };
0073
0074 //---------------------------------------------------------------------------//
0075 // FREE FUNCTIONS
0076 //---------------------------------------------------------------------------//
0077 /*!
0078  * Find orientation of ordered vertices in 2D coordinates.
0079  */
0080 inline Orientation
0081 calc_orientation(Real2 const& a, Real2 const& b, Real2 const& c)
0082 {
0083     auto crossp = (b[0] - a[0]) * (c[1] - b[1]) - (b[1] - a[1]) * (c[0] - b[0]);
0084     return crossp < 0   ? Orientation::clockwise
0085            : crossp > 0 ? Orientation::counterclockwise
0086                         : Orientation::collinear;
0087 }
0088
0089 //---------------------------------------------------------------------------//
0090 /*!
0091  * Test whether a 2D polygon has the given orientation.
0092  *
0093  * The list of input corners must have at least 3 points to be a polygon.
0094  */
0095 inline bool has_orientation(Span<Real2 const> corners, Orientation o)
0096 {
0097     CELER_EXPECT(corners.size() > 2);
0098     for (auto i : range(corners.size()))
0099     {
0100         auto j = (i + 1) % corners.size();
0101         auto k = (i + 2) % corners.size();
0102         if (calc_orientation(corners[i], corners[j], corners[k]) != o)
0103             return false;
0104     }
0105     return true;
0106 }
0107
0108 //---------------------------------------------------------------------------//
0109 /*!
0110  * Whether the orientation is the same or degenerate if allowed.
0111  */
0112 inline bool
0113 is_same_orientation(Orientation a, Orientation b, bool degen_ok = false)
0114 {
0115     if (a == Orientation::collinear || b == Orientation::collinear)
0116     {
0117         return degen_ok;
0118     }
0119     return (a == b);
0120 }
0121
0122 //---------------------------------------------------------------------------//
0123 /*!
0124  * Check if a 2D polygon is convex.
0125  *
0126  * \param corners the vertices of the polygon
0127  * \param degen_ok allow consecutive collinear points
0128  */
0129 inline bool is_convex(Span<Real2 const> corners, bool degen_ok = false)
0130 {
0131     CELER_EXPECT(corners.size() > 2);
0132     auto ref = Orientation::collinear;
0133     for (auto i : range<size_type>(corners.size()))
0134     {
0135         auto j = (i + 1) % corners.size();
0136         auto k = (i + 2) % corners.size();
0137         auto cur = calc_orientation(corners[i], corners[j], corners[k]);
0138         if (ref == Orientation::collinear)
0139         {
0140             // First non-collinear point
0141             ref = cur;
0142         }
0143         if (!is_same_orientation(cur, ref, degen_ok))
0144         {
0145             // Prohibited collinear orientation, or different orientation from
0146             // reference
0147             return false;
0148         }
0149     }
0150     return true;
0151 }
0152
0153 //---------------------------------------------------------------------------//
0154 }  // namespace detail
0155 }  // namespace orangeinp
0156 }  // namespace celeritas