ceph/src/boost/boost/geometry/strategy/geographic/area_box.hpp

   1 // Boost.Geometry
   2
   3 // Copyright (c) 2021, Oracle and/or its affiliates.
   4
   5 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
   6
   7 // Licensed under the Boost Software License version 1.0.
   8 // http://www.boost.org/users/license.html
   9
  10 #ifndef BOOST_GEOMETRY_STRATEGY_GEOGRAPHIC_AREA_BOX_HPP
  11 #define BOOST_GEOMETRY_STRATEGY_GEOGRAPHIC_AREA_BOX_HPP
  12
  13
  14 #include <boost/geometry/core/radian_access.hpp>
  15 #include <boost/geometry/srs/spheroid.hpp>
  16 #include <boost/geometry/strategies/spherical/get_radius.hpp>
  17 #include <boost/geometry/strategy/area.hpp>
  18 #include <boost/geometry/util/normalize_spheroidal_box_coordinates.hpp>
  19
  20
  21 namespace boost { namespace geometry
  22 {
  23
  24 namespace strategy { namespace area
  25 {
  26
  27 // Based on the approach for spherical coordinate system:
  28 // https://math.stackexchange.com/questions/131735/surface-element-in-spherical-coordinates
  29 // http://www.cs.cmu.edu/afs/cs/academic/class/16823-s16/www/pdfs/appearance-modeling-3.pdf
  30 // https://www.astronomyclub.xyz/celestial-sphere-2/solid-angle-on-the-celestial-sphere.html
  31 // https://mathworld.wolfram.com/SolidAngle.html
  32 // https://en.wikipedia.org/wiki/Spherical_coordinate_system
  33 // and equations for spheroid:
  34 // https://en.wikipedia.org/wiki/Geographic_coordinate_conversion
  35 // https://en.wikipedia.org/wiki/Meridian_arc
  36 // Note that the equations use geodetic latitudes so we do not have to convert them.
  37 // assume(y_max > y_min);
  38 // assume(x_max > x_min);
  39 // M: a*(1-e^2) / (1-e^2*sin(y)^2)^(3/2);
  40 // N: a / sqrt(1-e^2*sin(y)^2);
  41 // O: N*cos(y)*M;
  42 // tellsimp(log(abs(e*sin(y_min)+1)), p_min);
  43 // tellsimp(log(abs(e*sin(y_min)-1)), m_min);
  44 // tellsimp(log(abs(e*sin(y_max)+1)), p_max);
  45 // tellsimp(log(abs(e*sin(y_max)-1)), m_max);
  46 // S: integrate(integrate(O, y, y_min, y_max), x, x_min, x_max);
  47 // combine(S);
  48 //
  49 // An alternative solution to the above formula was suggested by Charles Karney
  50 // https://github.com/boostorg/geometry/pull/832
  51 // The following are formulas for area of a box defined by the equator and some latitude,
  52 // not arbitrary box.
  53 // For e^2 > 0
  54 // dlambda*b^2*sin(phi)/2*(1/(1-e^2*sin(phi)^2) + atanh(e*sin(phi))/(e*sin(phi)))
  55 // For e^2 < 0
  56 // dlambda*b^2*sin(phi)/2*(1/(1-e^2*sin(phi)^2) + atan(ea*sin(phi))/(ea*sin(phi)))
  57 // where ea = sqrt(-e^2)
  58 template
  59 <
  60     typename Spheroid = srs::spheroid<double>,
  61     typename CalculationType = void
  62 >
  63 class geographic_box
  64 {
  65 public:
  66     template <typename Box>
  67     struct result_type
  68         : strategy::area::detail::result_type
  69             <
  70                 Box,
  71                 CalculationType
  72             >
  73     {};
  74
  75     geographic_box() = default;
  76
  77     explicit geographic_box(Spheroid const& spheroid)
  78         : m_spheroid(spheroid)
  79     {}
  80
  81     template <typename Box>
  82     inline auto apply(Box const& box) const
  83     {
  84         typedef typename result_type<Box>::type return_type;
  85
  86         return_type const c0 = 0;
  87
  88         return_type x_min = get_as_radian<min_corner, 0>(box); // lon
  89         return_type y_min = get_as_radian<min_corner, 1>(box); // lat
  90         return_type x_max = get_as_radian<max_corner, 0>(box);
  91         return_type y_max = get_as_radian<max_corner, 1>(box);
  92
  93         math::normalize_spheroidal_box_coordinates<radian>(x_min, y_min, x_max, y_max);
  94
  95         if (x_min == x_max || y_max == y_min)
  96         {
  97             return c0;
  98         }
  99
 100         return_type const e2 = formula::eccentricity_sqr<return_type>(m_spheroid);
 101
 102         return_type const x_diff = x_max - x_min;
 103         return_type const sin_y_min = sin(y_min);
 104         return_type const sin_y_max = sin(y_max);
 105
 106         if (math::equals(e2, c0))
 107         {
 108             // spherical formula
 109             return_type const a = get_radius<0>(m_spheroid);
 110             return x_diff * (sin_y_max - sin_y_min) * a * a;
 111         }
 112
 113         return_type const c1 = 1;
 114         return_type const c2 = 2;
 115         return_type const b = get_radius<2>(m_spheroid);
 116
 117         /*
 118         return_type const c4 = 4;
 119         return_type const e = math::sqrt(e2);
 120
 121         return_type const p_min = log(math::abs(e * sin_y_min + c1));
 122         return_type const p_max = log(math::abs(e * sin_y_max + c1));
 123         return_type const m_min = log(math::abs(e * sin_y_min - c1));
 124         return_type const m_max = log(math::abs(e * sin_y_max - c1));
 125         return_type const n_min = e * sin_y_min * sin_y_min;
 126         return_type const n_max = e * sin_y_max * sin_y_max;
 127         return_type const d_min = e * n_min - c1;
 128         return_type const d_max = e * n_max - c1;
 129
 130         // NOTE: For equal latitudes the original formula generated by maxima may give negative
 131         //   result. It's caused by the order of operations, so here they're rearranged for
 132         //   symmetry.
 133         return_type const comp0 = (p_min - m_min) / (c4 * e * d_min);
 134         return_type const comp1 = sin_y_min / (c2 * d_min);
 135         return_type const comp2 = n_min * (m_min - p_min) / (c4 * d_min);
 136         return_type const comp3 = (p_max - m_max) / (c4 * e * d_max);
 137         return_type const comp4 = sin_y_max / (c2 * d_max);
 138         return_type const comp5 = n_max * (m_max - p_max) / (c4 * d_max);
 139         return_type const comp02 = comp0 + comp1 + comp2;
 140         return_type const comp35 = comp3 + comp4 + comp5;
 141
 142         return b * b * x_diff * (comp02 - comp35);
 143         */
 144
 145         return_type const comp0_min = c1 / (c1 - e2 * sin_y_min * sin_y_min);
 146         return_type const comp0_max = c1 / (c1 - e2 * sin_y_max * sin_y_max);
 147
 148         // NOTE: For latitudes equal to 0 the original formula returns NAN
 149         return_type comp1_min = 0, comp1_max = 0;
 150         if (e2 > c0)
 151         {
 152             return_type const e = math::sqrt(e2);
 153             return_type const e_sin_y_min = e * sin_y_min;
 154             return_type const e_sin_y_max = e * sin_y_max;
 155
 156             comp1_min = e_sin_y_min == c0 ? c1 : atanh(e_sin_y_min) / e_sin_y_min;
 157             comp1_max = e_sin_y_max == c0 ? c1 : atanh(e_sin_y_max) / e_sin_y_max;
 158         }
 159         else
 160         {
 161             return_type const ea = math::sqrt(-e2);
 162             return_type const ea_sin_y_min = ea * sin_y_min;
 163             return_type const ea_sin_y_max = ea * sin_y_max;
 164
 165             comp1_min = ea_sin_y_min == c0 ? c1 : atan(ea_sin_y_min) / ea_sin_y_min;
 166             comp1_max = ea_sin_y_max == c0 ? c1 : atan(ea_sin_y_max) / ea_sin_y_max;
 167         }
 168
 169         return_type const comp01_min = sin_y_min * (comp0_min + comp1_min);
 170         return_type const comp01_max = sin_y_max * (comp0_max + comp1_max);
 171
 172         return b * b * x_diff * (comp01_max - comp01_min) / c2;
 173     }
 174
 175     Spheroid model() const
 176     {
 177         return m_spheroid;
 178     }
 179
 180 private:
 181     Spheroid m_spheroid;
 182 };
 183
 184
 185 }} // namespace strategy::area
 186
 187
 188 }} // namespace boost::geometry
 189
 190
 191 #endif // BOOST_GEOMETRY_STRATEGY_GEOGRAPHIC_AREA_BOX_HPP