ceph/src/boost/libs/geometry/include/boost/geometry/strategies/spherical/area_huiller.hpp

   1 // Boost.Geometry (aka GGL, Generic Geometry Library)
   2
   3 // Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands.
   4
   5 // This file was modified by Oracle on 2015.
   6 // Modifications copyright (c) 2015, Oracle and/or its affiliates.
   7
   8 // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
   9 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
  10
  11 // Use, modification and distribution is subject to the Boost Software License,
  12 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  13 // http://www.boost.org/LICENSE_1_0.txt)
  14
  15 #ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AREA_HUILLER_HPP
  16 #define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AREA_HUILLER_HPP
  17
  18
  19
  20 #include <boost/geometry/strategies/spherical/distance_haversine.hpp>
  21
  22 #include <boost/geometry/core/radian_access.hpp>
  23 #include <boost/geometry/util/math.hpp>
  24
  25
  26 namespace boost { namespace geometry
  27 {
  28
  29 namespace strategy { namespace area
  30 {
  31
  32
  33
  34 /*!
  35 \brief Area calculation by spherical excess / Huiller's formula
  36 \ingroup strategies
  37 \tparam PointOfSegment point type of segments of rings/polygons
  38 \tparam CalculationType \tparam_calculation
  39 \author Barend Gehrels. Adapted from:
  40 - http://webdocs.cs.ualberta.ca/~graphics/books/GraphicsGems/gemsiv/sph_poly.c
  41 - http://tog.acm.org/resources/GraphicsGems/gemsiv/sph_poly.c
  42 - http://williams.best.vwh.net/avform.htm
  43 \note The version in Graphics Gems IV (page 132-137) didn't account for
  44 polygons crossing the 0 and 180 meridians. The fix for this algorithm
  45 can be found in Graphics Gems V (pages 45-46). See:
  46 - http://kysmykseka.net/koti/wizardry/Game%20Development/Programming/Graphics%20Gems%204.pdf
  47 - http://kysmykseka.net/koti/wizardry/Game%20Development/Programming/Graphics%20Gems%205.pdf
  48 \note This version works for convex and non-convex polygons, for 180 meridian
  49 crossing polygons and for polygons with holes. However, some cases (especially
  50 180 meridian cases) must still be checked.
  51 \note The version which sums angles, which is often seen, doesn't handle non-convex
  52 polygons correctly.
  53 \note The version which sums longitudes, see http://hdl.handle.net/2014/40409,
  54 is simple and works well in most cases but not in 180 meridian crossing cases.
  55 This probably could be solved.
  56
  57 \note This version is made for spherical equatorial coordinate systems
  58
  59 \qbk{
  60
  61 [heading Example]
  62 [area_with_strategy]
  63 [area_with_strategy_output]
  64
  65
  66 [heading See also]
  67 [link geometry.reference.algorithms.area.area_2_with_strategy area (with strategy)]
  68 }
  69
  70 */
  71 template
  72 <
  73     typename PointOfSegment,
  74     typename CalculationType = void
  75 >
  76 class huiller
  77 {
  78 typedef typename boost::mpl::if_c
  79     <
  80         boost::is_void<CalculationType>::type::value,
  81         typename select_most_precise
  82             <
  83                 typename coordinate_type<PointOfSegment>::type,
  84                 double
  85             >::type,
  86         CalculationType
  87     >::type calculation_type;
  88
  89 protected :
  90     struct excess_sum
  91     {
  92         calculation_type sum;
  93
  94         // Distances are calculated on unit sphere here
  95         strategy::distance::haversine<calculation_type> distance_over_unit_sphere;
  96
  97
  98         inline excess_sum()
  99             : sum(0)
 100             , distance_over_unit_sphere(1)
 101         {}
 102         inline calculation_type area(calculation_type radius) const
 103         {
 104             return - sum * radius * radius;
 105         }
 106     };
 107
 108 public :
 109     typedef calculation_type return_type;
 110     typedef PointOfSegment segment_point_type;
 111     typedef excess_sum state_type;
 112
 113     inline huiller(calculation_type radius = 1.0)
 114         : m_radius(radius)
 115     {}
 116
 117     inline void apply(PointOfSegment const& p1,
 118                 PointOfSegment const& p2,
 119                 excess_sum& state) const
 120     {
 121         if (! geometry::math::equals(get<0>(p1), get<0>(p2)))
 122         {
 123             calculation_type const half = 0.5;
 124             calculation_type const two = 2.0;
 125             calculation_type const four = 4.0;
 126             calculation_type const pi
 127                 = geometry::math::pi<calculation_type>();
 128             calculation_type const two_pi
 129                 = geometry::math::two_pi<calculation_type>();
 130             calculation_type const half_pi
 131                 = geometry::math::half_pi<calculation_type>();
 132
 133             // Distance p1 p2
 134             calculation_type a = state.distance_over_unit_sphere.apply(p1, p2);
 135
 136             // Sides on unit sphere to south pole
 137             calculation_type b = half_pi - geometry::get_as_radian<1>(p2);
 138             calculation_type c = half_pi - geometry::get_as_radian<1>(p1);
 139
 140             // Semi parameter
 141             calculation_type s = half * (a + b + c);
 142
 143             // E: spherical excess, using l'Huiller's formula
 144             // [tg(e / 4)]2   =   tg[s / 2]  tg[(s-a) / 2]  tg[(s-b) / 2]  tg[(s-c) / 2]
 145             calculation_type excess = four
 146                 * atan(geometry::math::sqrt(geometry::math::abs(tan(s / two)
 147                     * tan((s - a) / two)
 148                     * tan((s - b) / two)
 149                     * tan((s - c) / two))));
 150
 151             excess = geometry::math::abs(excess);
 152
 153             // In right direction: positive, add area. In left direction: negative, subtract area.
 154             // Longitude comparisons are not so obvious. If one is negative and other is positive,
 155             // we have to take the dateline into account.
 156
 157             calculation_type lon_diff = geometry::get_as_radian<0>(p2)
 158                                       - geometry::get_as_radian<0>(p1);
 159             if (lon_diff <= 0)
 160             {
 161                 lon_diff += two_pi;
 162             }
 163
 164             if (lon_diff > pi)
 165             {
 166                 excess = -excess;
 167             }
 168
 169             state.sum += excess;
 170         }
 171     }
 172
 173     inline return_type result(excess_sum const& state) const
 174     {
 175         return state.area(m_radius);
 176     }
 177
 178 private :
 179     /// Radius of the sphere
 180     calculation_type m_radius;
 181 };
 182
 183 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
 184
 185 namespace services
 186 {
 187
 188
 189 template <typename Point>
 190 struct default_strategy<spherical_equatorial_tag, Point>
 191 {
 192     typedef strategy::area::huiller<Point> type;
 193 };
 194
 195 // Note: spherical polar coordinate system requires "get_as_radian_equatorial"
 196 /***template <typename Point>
 197 struct default_strategy<spherical_polar_tag, Point>
 198 {
 199     typedef strategy::area::huiller<Point> type;
 200 };***/
 201
 202 } // namespace services
 203
 204 #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
 205
 206
 207 }} // namespace strategy::area
 208
 209
 210
 211
 212 }} // namespace boost::geometry
 213
 214 #endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AREA_HUILLER_HPP