ceph/src/boost/libs/geometry/include/boost/geometry/algorithms/detail/envelope/segment.hpp

   1 // Boost.Geometry (aka GGL, Generic Geometry Library)
   2
   3 // Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands.
   4 // Copyright (c) 2008-2015 Bruno Lalande, Paris, France.
   5 // Copyright (c) 2009-2015 Mateusz Loskot, London, UK.
   6
   7 // This file was modified by Oracle on 2015, 2016.
   8 // Modifications copyright (c) 2015-2016, Oracle and/or its affiliates.
   9
  10 // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
  11 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
  12
  13 // Distributed under the Boost Software License, Version 1.0.
  14 // (See accompanying file LICENSE_1_0.txt or copy at
  15 // http://www.boost.org/LICENSE_1_0.txt)
  16
  17 #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_ENVELOPE_SEGMENT_HPP
  18 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_ENVELOPE_SEGMENT_HPP
  19
  20 #include <cstddef>
  21 #include <utility>
  22
  23 #include <boost/numeric/conversion/cast.hpp>
  24
  25 #include <boost/geometry/core/assert.hpp>
  26 #include <boost/geometry/core/coordinate_system.hpp>
  27 #include <boost/geometry/core/coordinate_type.hpp>
  28 #include <boost/geometry/core/cs.hpp>
  29 #include <boost/geometry/core/point_type.hpp>
  30 #include <boost/geometry/core/radian_access.hpp>
  31 #include <boost/geometry/core/tags.hpp>
  32
  33 #include <boost/geometry/util/math.hpp>
  34
  35 #include <boost/geometry/geometries/helper_geometry.hpp>
  36
  37 #include <boost/geometry/strategies/compare.hpp>
  38
  39 #include <boost/geometry/algorithms/detail/assign_indexed_point.hpp>
  40 #include <boost/geometry/algorithms/detail/normalize.hpp>
  41
  42 #include <boost/geometry/algorithms/detail/envelope/point.hpp>
  43 #include <boost/geometry/algorithms/detail/envelope/transform_units.hpp>
  44
  45 #include <boost/geometry/algorithms/detail/expand/point.hpp>
  46
  47 #include <boost/geometry/algorithms/dispatch/envelope.hpp>
  48
  49
  50 namespace boost { namespace geometry
  51 {
  52
  53 #ifndef DOXYGEN_NO_DETAIL
  54 namespace detail { namespace envelope
  55 {
  56
  57
  58 template <std::size_t Dimension, std::size_t DimensionCount>
  59 struct envelope_one_segment
  60 {
  61     template<typename Point, typename Box>
  62     static inline void apply(Point const& p1, Point const& p2, Box& mbr)
  63     {
  64         envelope_one_point<Dimension, DimensionCount>::apply(p1, mbr);
  65         detail::expand::point_loop
  66             <
  67                 strategy::compare::default_strategy,
  68                 strategy::compare::default_strategy,
  69                 Dimension,
  70                 DimensionCount
  71             >::apply(mbr, p2);
  72     }
  73 };
  74
  75
  76 // Computes the MBR of a segment given by (lon1, lat1) and (lon2,
  77 // lat2), and with azimuths a1 and a2 at the two endpoints of the
  78 // segment.
  79 // It is assumed that the spherical coordinates of the segment are
  80 // normalized and in radians.
  81 // The longitudes and latitudes of the endpoints are overridden by
  82 // those of the box.
  83 class compute_mbr_of_segment
  84 {
  85 private:
  86     // computes the azimuths of the segment with endpoints (lon1, lat1)
  87     // and (lon2, lat2)
  88     // radians
  89     template <typename CalculationType>
  90     static inline void azimuths(CalculationType const& lon1,
  91                                 CalculationType const& lat1,
  92                                 CalculationType const& lon2,
  93                                 CalculationType const& lat2,
  94                                 CalculationType& a1,
  95                                 CalculationType& a2)
  96     {
  97         BOOST_GEOMETRY_ASSERT(lon1 <= lon2);
  98
  99         CalculationType dlon = lon2 - lon1;
 100         CalculationType sin_dlon = sin(dlon);
 101         CalculationType cos_dlon = cos(dlon);
 102         CalculationType cos_lat1 = cos(lat1);
 103         CalculationType cos_lat2 = cos(lat2);
 104         CalculationType sin_lat1 = sin(lat1);
 105         CalculationType sin_lat2 = sin(lat2);
 106
 107         a1 = atan2(sin_dlon * cos_lat2,
 108                    cos_lat1 * sin_lat2 - sin_lat1 * cos_lat2 * cos_dlon);
 109
 110         a2 = atan2(-sin_dlon * cos_lat1,
 111                    cos_lat2 * sin_lat1 - sin_lat2 * cos_lat1 * cos_dlon);
 112         a2 += math::pi<CalculationType>();
 113     }
 114
 115     // degrees or radians
 116     template <typename CalculationType>
 117     static inline void swap(CalculationType& lon1,
 118                             CalculationType& lat1,
 119                             CalculationType& lon2,
 120                             CalculationType& lat2)
 121     {
 122         std::swap(lon1, lon2);
 123         std::swap(lat1, lat2);
 124     }
 125
 126     // radians
 127     template <typename CalculationType>
 128     static inline bool contains_pi_half(CalculationType const& a1,
 129                                         CalculationType const& a2)
 130     {
 131         // azimuths a1 and a2 are assumed to be in radians
 132         BOOST_GEOMETRY_ASSERT(! math::equals(a1, a2));
 133
 134         static CalculationType const pi_half = math::half_pi<CalculationType>();
 135
 136         return (a1 < a2)
 137             ? (a1 < pi_half && pi_half < a2)
 138             : (a1 > pi_half && pi_half > a2);
 139     }
 140
 141     // radians or degrees
 142     template <typename Units, typename CoordinateType>
 143     static inline bool crosses_antimeridian(CoordinateType const& lon1,
 144                                             CoordinateType const& lon2)
 145     {
 146         typedef math::detail::constants_on_spheroid
 147             <
 148                 CoordinateType, Units
 149             > constants;
 150
 151         return math::abs(lon1 - lon2) > constants::half_period(); // > pi
 152     }
 153
 154     // radians
 155     template <typename CalculationType>
 156     static inline CalculationType max_latitude(CalculationType const& azimuth,
 157                                                CalculationType const& latitude)
 158     {
 159         // azimuth and latitude are assumed to be in radians
 160         return acos( math::abs(cos(latitude) * sin(azimuth)) );
 161     }
 162
 163     // degrees or radians
 164     template <typename Units, typename CalculationType>
 165     static inline void compute_box_corners(CalculationType& lon1,
 166                                            CalculationType& lat1,
 167                                            CalculationType& lon2,
 168                                            CalculationType& lat2)
 169     {
 170         // coordinates are assumed to be in radians
 171         BOOST_GEOMETRY_ASSERT(lon1 <= lon2);
 172
 173         CalculationType lon1_rad = math::as_radian<Units>(lon1);
 174         CalculationType lat1_rad = math::as_radian<Units>(lat1);
 175         CalculationType lon2_rad = math::as_radian<Units>(lon2);
 176         CalculationType lat2_rad = math::as_radian<Units>(lat2);
 177
 178         CalculationType a1 = 0, a2 = 0;
 179         azimuths(lon1_rad, lat1_rad, lon2_rad, lat2_rad, a1, a2);
 180
 181         if (lat1 > lat2)
 182         {
 183             std::swap(lat1, lat2);
 184             std::swap(lat1_rad, lat2_rad);
 185         }
 186
 187         if (math::equals(a1, a2))
 188         {
 189             // the segment must lie on the equator or is very short
 190             return;
 191         }
 192
 193         if (contains_pi_half(a1, a2))
 194         {
 195             CalculationType const mid_lat = lat1 + lat2;
 196             if (mid_lat < 0)
 197             {
 198                 // update using min latitude
 199                 CalculationType const lat_min_rad = -max_latitude(a1, lat1_rad);
 200                 CalculationType const lat_min = math::from_radian<Units>(lat_min_rad);
 201
 202                 if (lat1 > lat_min)
 203                 {
 204                     lat1 = lat_min;
 205                 }
 206             }
 207             else if (mid_lat > 0)
 208             {
 209                 // update using max latitude
 210                 CalculationType const lat_max_rad = max_latitude(a1, lat1_rad);
 211                 CalculationType const lat_max = math::from_radian<Units>(lat_max_rad);
 212
 213                 if (lat2 < lat_max)
 214                 {
 215                     lat2 = lat_max;
 216                 }
 217             }
 218         }
 219     }
 220
 221     template <typename Units, typename CalculationType>
 222     static inline void apply(CalculationType& lon1,
 223                              CalculationType& lat1,
 224                              CalculationType& lon2,
 225                              CalculationType& lat2)
 226     {
 227         typedef math::detail::constants_on_spheroid
 228             <
 229                 CalculationType, Units
 230             > constants;
 231
 232         bool is_pole1 = math::equals(math::abs(lat1), constants::max_latitude());
 233         bool is_pole2 = math::equals(math::abs(lat2), constants::max_latitude());
 234
 235         if (is_pole1 && is_pole2)
 236         {
 237             // both points are poles; nothing more to do:
 238             // longitudes are already normalized to 0
 239             // but just in case
 240             lon1 = 0;
 241             lon2 = 0;
 242         }
 243         else if (is_pole1 && !is_pole2)
 244         {
 245             // first point is a pole, second point is not:
 246             // make the longitude of the first point the same as that
 247             // of the second point
 248             lon1 = lon2;
 249         }
 250         else if (!is_pole1 && is_pole2)
 251         {
 252             // second point is a pole, first point is not:
 253             // make the longitude of the second point the same as that
 254             // of the first point
 255             lon2 = lon1;
 256         }
 257
 258         if (lon1 == lon2)
 259         {
 260             // segment lies on a meridian
 261             if (lat1 > lat2)
 262             {
 263                 std::swap(lat1, lat2);
 264             }
 265             return;
 266         }
 267
 268         BOOST_GEOMETRY_ASSERT(!is_pole1 && !is_pole2);
 269
 270         if (lon1 > lon2)
 271         {
 272             swap(lon1, lat1, lon2, lat2);
 273         }
 274
 275         if (crosses_antimeridian<Units>(lon1, lon2))
 276         {
 277             lon1 += constants::period();
 278             swap(lon1, lat1, lon2, lat2);
 279         }
 280
 281         compute_box_corners<Units>(lon1, lat1, lon2, lat2);
 282     }
 283
 284 public:
 285     template <typename Units, typename CalculationType, typename Box>
 286     static inline void apply(CalculationType lon1,
 287                              CalculationType lat1,
 288                              CalculationType lon2,
 289                              CalculationType lat2,
 290                              Box& mbr)
 291     {
 292         typedef typename coordinate_type<Box>::type box_coordinate_type;
 293
 294         typedef typename helper_geometry
 295             <
 296                 Box, box_coordinate_type, Units
 297             >::type helper_box_type;
 298
 299         helper_box_type radian_mbr;
 300
 301         apply<Units>(lon1, lat1, lon2, lat2);
 302
 303         geometry::set
 304             <
 305                 min_corner, 0
 306             >(radian_mbr, boost::numeric_cast<box_coordinate_type>(lon1));
 307
 308         geometry::set
 309             <
 310                 min_corner, 1
 311             >(radian_mbr, boost::numeric_cast<box_coordinate_type>(lat1));
 312
 313         geometry::set
 314             <
 315                 max_corner, 0
 316             >(radian_mbr, boost::numeric_cast<box_coordinate_type>(lon2));
 317
 318         geometry::set
 319             <
 320                 max_corner, 1
 321             >(radian_mbr, boost::numeric_cast<box_coordinate_type>(lat2));
 322
 323         transform_units(radian_mbr, mbr);
 324     }
 325 };
 326
 327
 328 template <std::size_t DimensionCount>
 329 struct envelope_segment_on_sphere
 330 {
 331     template <typename Point, typename Box>
 332     static inline void apply(Point const& p1, Point const& p2, Box& mbr)
 333     {
 334         // first compute the envelope range for the first two coordinates
 335         Point p1_normalized = detail::return_normalized<Point>(p1);
 336         Point p2_normalized = detail::return_normalized<Point>(p2);
 337
 338         typedef typename coordinate_system<Point>::type::units units_type;
 339
 340         compute_mbr_of_segment::template apply<units_type>(
 341                     geometry::get<0>(p1_normalized),
 342                     geometry::get<1>(p1_normalized),
 343                     geometry::get<0>(p2_normalized),
 344                     geometry::get<1>(p2_normalized),
 345                     mbr);
 346
 347         // now compute the envelope range for coordinates of
 348         // dimension 2 and higher
 349         envelope_one_segment<2, DimensionCount>::apply(p1, p2, mbr);
 350     }
 351
 352     template <typename Segment, typename Box>
 353     static inline void apply(Segment const& segment, Box& mbr)
 354     {
 355         typename point_type<Segment>::type p[2];
 356         detail::assign_point_from_index<0>(segment, p[0]);
 357         detail::assign_point_from_index<1>(segment, p[1]);
 358         apply(p[0], p[1], mbr);
 359     }
 360 };
 361
 362
 363
 364 template <std::size_t DimensionCount, typename CS_Tag>
 365 struct envelope_segment
 366     : envelope_one_segment<0, DimensionCount>
 367 {};
 368
 369
 370 template <std::size_t DimensionCount>
 371 struct envelope_segment<DimensionCount, spherical_equatorial_tag>
 372     : envelope_segment_on_sphere<DimensionCount>
 373 {};
 374
 375
 376
 377 }} // namespace detail::envelope
 378 #endif // DOXYGEN_NO_DETAIL
 379
 380
 381 #ifndef DOXYGEN_NO_DISPATCH
 382 namespace dispatch
 383 {
 384
 385
 386 template <typename Segment, typename CS_Tag>
 387 struct envelope<Segment, segment_tag, CS_Tag>
 388 {
 389     template <typename Box>
 390     static inline void apply(Segment const& segment, Box& mbr)
 391     {
 392         typename point_type<Segment>::type p[2];
 393         detail::assign_point_from_index<0>(segment, p[0]);
 394         detail::assign_point_from_index<1>(segment, p[1]);
 395         detail::envelope::envelope_segment
 396             <
 397                 dimension<Segment>::value, CS_Tag
 398             >::apply(p[0], p[1], mbr);
 399     }
 400 };
 401
 402
 403 } // namespace dispatch
 404 #endif // DOXYGEN_NO_DISPATCH
 405
 406 }} // namespace boost::geometry
 407
 408 #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_ENVELOPE_SEGMENT_HPP