ceph/src/boost/libs/geometry/include/boost/geometry/strategies/geographic/side_detail.hpp

   1 // Boost.Geometry
   2
   3 // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
   4
   5 // This file was modified by Oracle on 2014, 2015.
   6 // Modifications copyright (c) 2014-2015 Oracle and/or its affiliates.
   7
   8 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
   9
  10 // Use, modification and distribution is subject to the Boost Software License,
  11 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  12 // http://www.boost.org/LICENSE_1_0.txt)
  13
  14 #ifndef BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_SIDE_DETAIL_HPP
  15 #define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_SIDE_DETAIL_HPP
  16
  17 #include <boost/geometry/core/cs.hpp>
  18 #include <boost/geometry/core/access.hpp>
  19 #include <boost/geometry/core/radian_access.hpp>
  20 #include <boost/geometry/core/radius.hpp>
  21
  22 #include <boost/geometry/util/math.hpp>
  23 #include <boost/geometry/util/promote_floating_point.hpp>
  24 #include <boost/geometry/util/select_calculation_type.hpp>
  25
  26 #include <boost/geometry/strategies/side.hpp>
  27 //#include <boost/geometry/strategies/concepts/side_concept.hpp>
  28
  29
  30 namespace boost { namespace geometry
  31 {
  32
  33
  34 namespace strategy { namespace side
  35 {
  36
  37 #ifndef DOXYGEN_NO_DETAIL
  38 namespace detail
  39 {
  40
  41 /*!
  42 \brief Check at which side of a segment a point lies
  43          left of segment (> 0), right of segment (< 0), on segment (0)
  44 \ingroup strategies
  45 \tparam InverseFormula Geodesic inverse solution formula.
  46 \tparam Model Reference model of coordinate system.
  47 \tparam CalculationType \tparam_calculation
  48  */
  49 template <template<typename, bool, bool> class InverseFormula,
  50           typename Model,
  51           typename CalculationType = void>
  52 class by_azimuth
  53 {
  54 public:
  55     by_azimuth(Model const& model = Model())
  56         : m_model(model)
  57     {}
  58
  59     template <typename P1, typename P2, typename P>
  60     inline int apply(P1 const& p1, P2 const& p2, P const& p)
  61     {
  62         typedef typename promote_floating_point
  63             <
  64                 typename select_calculation_type_alt
  65                     <
  66                         CalculationType,
  67                         P1, P2, P
  68                     >::type
  69             >::type calc_t;
  70
  71         typedef InverseFormula<calc_t, false, true> inverse_formula;
  72
  73         calc_t a1p = azimuth<calc_t, inverse_formula>(p1, p, m_model);
  74         calc_t a12 = azimuth<calc_t, inverse_formula>(p1, p2, m_model);
  75
  76         calc_t const pi = math::pi<calc_t>();
  77
  78         // instead of the formula from XTD
  79         //calc_t a_diff = asin(sin(a1p - a12));
  80
  81         calc_t a_diff = a1p - a12;
  82         // normalize, angle in [-pi, pi]
  83         while ( a_diff > pi )
  84             a_diff -= calc_t(2) * pi;
  85         while ( a_diff < -pi )
  86             a_diff += calc_t(2) * pi;
  87
  88         // NOTE: in general it shouldn't be required to support the pi/-pi case
  89         // because in non-cartesian systems it makes sense to check the side
  90         // only "between" the endpoints.
  91         // However currently the winding strategy calls the side strategy
  92         // for vertical segments to check if the point is "between the endpoints.
  93         // This could be avoided since the side strategy is not required for that
  94         // because meridian is the shortest path. So a difference of
  95         // longitudes would be sufficient (of course normalized to [-pi, pi]).
  96
  97         // NOTE: with the above said, the pi/-pi check is temporary
  98         // however in case if this was required
  99         // the geodesics on ellipsoid aren't "symmetrical"
 100         // therefore instead of comparing a_diff to pi and -pi
 101         // one should probably use inverse azimuths and compare
 102         // the difference to 0 as well
 103
 104         // positive azimuth is on the right side
 105         return math::equals(a_diff, 0)
 106             || math::equals(a_diff, pi)
 107             || math::equals(a_diff, -pi) ? 0
 108              : a_diff > 0 ? -1 // right
 109              : 1; // left
 110     }
 111
 112 private:
 113     template <typename ResultType,
 114               typename InverseFormulaType,
 115               typename Point1,
 116               typename Point2,
 117               typename ModelT>
 118     static inline ResultType azimuth(Point1 const& point1, Point2 const& point2, ModelT const& model)
 119     {
 120         return InverseFormulaType::apply(get_as_radian<0>(point1),
 121                                          get_as_radian<1>(point1),
 122                                          get_as_radian<0>(point2),
 123                                          get_as_radian<1>(point2),
 124                                          model).azimuth;
 125     }
 126
 127     Model m_model;
 128 };
 129
 130 } // detail
 131 #endif // DOXYGEN_NO_DETAIL
 132
 133 }} // namespace strategy::side
 134
 135
 136 }} // namespace boost::geometry
 137
 138
 139 #endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_SIDE_DETAIL_HPP