ceph/src/boost/boost/geometry/strategies/geographic/intersection.hpp

   1 // Boost.Geometry
   2
   3 // Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland.
   4
   5 // Copyright (c) 2016-2021, Oracle and/or its affiliates.
   6 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
   7
   8 // Use, modification and distribution is subject to the Boost Software License,
   9 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  10 // http://www.boost.org/LICENSE_1_0.txt)
  11
  12 #ifndef BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_INTERSECTION_HPP
  13 #define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_INTERSECTION_HPP
  14
  15 #include <algorithm>
  16 #include <type_traits>
  17
  18 #include <boost/geometry/core/cs.hpp>
  19 #include <boost/geometry/core/access.hpp>
  20 #include <boost/geometry/core/radian_access.hpp>
  21 #include <boost/geometry/core/tags.hpp>
  22
  23 #include <boost/geometry/algorithms/detail/assign_values.hpp>
  24 #include <boost/geometry/algorithms/detail/assign_indexed_point.hpp>
  25 #include <boost/geometry/algorithms/detail/equals/point_point.hpp>
  26 #include <boost/geometry/algorithms/detail/recalculate.hpp>
  27
  28 #include <boost/geometry/formulas/andoyer_inverse.hpp>
  29 #include <boost/geometry/formulas/sjoberg_intersection.hpp>
  30 #include <boost/geometry/formulas/spherical.hpp>
  31 #include <boost/geometry/formulas/unit_spheroid.hpp>
  32
  33 #include <boost/geometry/geometries/concepts/point_concept.hpp>
  34 #include <boost/geometry/geometries/concepts/segment_concept.hpp>
  35 #include <boost/geometry/geometries/segment.hpp>
  36
  37 #include <boost/geometry/policies/robustness/segment_ratio.hpp>
  38
  39 #include <boost/geometry/srs/spheroid.hpp>
  40
  41 #include <boost/geometry/strategy/geographic/area.hpp>
  42 #include <boost/geometry/strategy/geographic/envelope.hpp>
  43 #include <boost/geometry/strategy/geographic/expand_segment.hpp>
  44 #include <boost/geometry/strategy/spherical/expand_box.hpp>
  45
  46 #include <boost/geometry/strategies/geographic/disjoint_segment_box.hpp>
  47 #include <boost/geometry/strategies/geographic/distance.hpp>
  48 #include <boost/geometry/strategies/geographic/parameters.hpp>
  49 #include <boost/geometry/strategies/geographic/point_in_poly_winding.hpp>
  50 #include <boost/geometry/strategies/geographic/side.hpp>
  51 #include <boost/geometry/strategies/spherical/disjoint_box_box.hpp>
  52 #include <boost/geometry/strategies/spherical/point_in_point.hpp>
  53 #include <boost/geometry/strategies/intersection.hpp>
  54 #include <boost/geometry/strategies/intersection_result.hpp>
  55 #include <boost/geometry/strategies/side_info.hpp>
  56
  57 #include <boost/geometry/util/math.hpp>
  58 #include <boost/geometry/util/select_calculation_type.hpp>
  59
  60
  61 namespace boost { namespace geometry
  62 {
  63
  64 namespace strategy { namespace intersection
  65 {
  66
  67 // CONSIDER: Improvement of the robustness/accuracy/repeatability by
  68 // moving all segments to 0 longitude
  69 // picking latitudes closer to 0
  70 // etc.
  71
  72 template
  73 <
  74     typename FormulaPolicy = strategy::andoyer,
  75     std::size_t Order = strategy::default_order<FormulaPolicy>::value,
  76     typename Spheroid = srs::spheroid<double>,
  77     typename CalculationType = void
  78 >
  79 struct geographic_segments
  80 {
  81     typedef geographic_tag cs_tag;
  82
  83     enum intersection_point_flag { ipi_inters = 0, ipi_at_a1, ipi_at_a2, ipi_at_b1, ipi_at_b2 };
  84
  85     template <typename CoordinateType, typename SegmentRatio>
  86     struct segment_intersection_info
  87     {
  88         template <typename Point, typename Segment1, typename Segment2>
  89         void calculate(Point& point, Segment1 const& a, Segment2 const& b) const
  90         {
  91             if (ip_flag == ipi_inters)
  92             {
  93                 // TODO: assign the rest of coordinates
  94                 set_from_radian<0>(point, lon);
  95                 set_from_radian<1>(point, lat);
  96             }
  97             else if (ip_flag == ipi_at_a1)
  98             {
  99                 detail::assign_point_from_index<0>(a, point);
 100             }
 101             else if (ip_flag == ipi_at_a2)
 102             {
 103                 detail::assign_point_from_index<1>(a, point);
 104             }
 105             else if (ip_flag == ipi_at_b1)
 106             {
 107                 detail::assign_point_from_index<0>(b, point);
 108             }
 109             else // ip_flag == ipi_at_b2
 110             {
 111                 detail::assign_point_from_index<1>(b, point);
 112             }
 113         }
 114
 115         CoordinateType lon;
 116         CoordinateType lat;
 117         SegmentRatio robust_ra;
 118         SegmentRatio robust_rb;
 119         intersection_point_flag ip_flag;
 120     };
 121
 122     explicit geographic_segments(Spheroid const& spheroid = Spheroid())
 123         : m_spheroid(spheroid)
 124     {}
 125
 126     Spheroid model() const
 127     {
 128         return m_spheroid;
 129     }
 130
 131     // Relate segments a and b
 132     template
 133     <
 134         typename UniqueSubRange1,
 135         typename UniqueSubRange2,
 136         typename Policy
 137     >
 138     inline typename Policy::return_type apply(UniqueSubRange1 const& range_p,
 139                                               UniqueSubRange2 const& range_q,
 140                                               Policy const&) const
 141     {
 142         typedef typename UniqueSubRange1::point_type point1_type;
 143         typedef typename UniqueSubRange2::point_type point2_type;
 144         typedef model::referring_segment<point1_type const> segment_type1;
 145         typedef model::referring_segment<point2_type const> segment_type2;
 146
 147         BOOST_CONCEPT_ASSERT( (concepts::ConstPoint<point1_type>) );
 148         BOOST_CONCEPT_ASSERT( (concepts::ConstPoint<point2_type>) );
 149
 150         /*
 151         typename coordinate_type<Point1>::type
 152             const a1_lon = get<0>(a1),
 153             const a2_lon = get<0>(a2);
 154         typename coordinate_type<Point2>::type
 155             const b1_lon = get<0>(b1),
 156             const b2_lon = get<0>(b2);
 157         bool is_a_reversed = a1_lon > a2_lon || a1_lon == a2_lon && get<1>(a1) > get<1>(a2);
 158         bool is_b_reversed = b1_lon > b2_lon || b1_lon == b2_lon && get<1>(b1) > get<1>(b2);
 159         */
 160
 161         point1_type const& p0 = range_p.at(0);
 162         point1_type const& p1 = range_p.at(1);
 163         point2_type const& q0 = range_q.at(0);
 164         point2_type const& q1 = range_q.at(1);
 165
 166         bool const is_p_reversed = get<1>(p0) > get<1>(p1);
 167         bool const is_q_reversed = get<1>(q0) > get<1>(q1);
 168
 169         // Call apply with original segments and ordered points
 170         return apply<Policy>(segment_type1(p0, p1),
 171                              segment_type2(q0, q1),
 172                              (is_p_reversed ? p1 : p0),
 173                              (is_p_reversed ? p0 : p1),
 174                              (is_q_reversed ? q1 : q0),
 175                              (is_q_reversed ? q0 : q1),
 176                              is_p_reversed, is_q_reversed);
 177     }
 178
 179 private:
 180     // Relate segments a and b
 181     template
 182     <
 183         typename Policy,
 184         typename Segment1,
 185         typename Segment2,
 186         typename Point1,
 187         typename Point2
 188     >
 189     inline typename Policy::return_type apply(Segment1 const& a, Segment2 const& b,
 190                                               Point1 const& a1, Point1 const& a2,
 191                                               Point2 const& b1, Point2 const& b2,
 192                                               bool is_a_reversed, bool is_b_reversed) const
 193     {
 194         BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment1>) );
 195         BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment2>) );
 196
 197         typedef typename select_calculation_type
 198             <Segment1, Segment2, CalculationType>::type calc_t;
 199
 200         typedef srs::spheroid<calc_t> spheroid_type;
 201
 202         static const calc_t c0 = 0;
 203
 204         // normalized spheroid
 205         spheroid_type spheroid = formula::unit_spheroid<spheroid_type>(m_spheroid);
 206
 207         // TODO: check only 2 first coordinates here?
 208         bool a_is_point = equals_point_point(a1, a2);
 209         bool b_is_point = equals_point_point(b1, b2);
 210
 211         if(a_is_point && b_is_point)
 212         {
 213             return equals_point_point(a1, b2)
 214                 ? Policy::degenerate(a, true)
 215                 : Policy::disjoint()
 216                 ;
 217         }
 218
 219         calc_t const a1_lon = get_as_radian<0>(a1);
 220         calc_t const a1_lat = get_as_radian<1>(a1);
 221         calc_t const a2_lon = get_as_radian<0>(a2);
 222         calc_t const a2_lat = get_as_radian<1>(a2);
 223         calc_t const b1_lon = get_as_radian<0>(b1);
 224         calc_t const b1_lat = get_as_radian<1>(b1);
 225         calc_t const b2_lon = get_as_radian<0>(b2);
 226         calc_t const b2_lat = get_as_radian<1>(b2);
 227
 228         side_info sides;
 229
 230         // NOTE: potential optimization, don't calculate distance at this point
 231         // this would require to reimplement inverse strategy to allow
 232         // calculation of distance if needed, probably also storing intermediate
 233         // results somehow inside an object.
 234         typedef typename FormulaPolicy::template inverse<calc_t, true, true, false, false, false> inverse_dist_azi;
 235         typedef typename inverse_dist_azi::result_type inverse_result;
 236
 237         // TODO: no need to call inverse formula if we know that the points are equal
 238         // distance can be set to 0 in this case and azimuth may be not calculated
 239         bool is_equal_a1_b1 = equals_point_point(a1, b1);
 240         bool is_equal_a2_b1 = equals_point_point(a2, b1);
 241         bool degen_neq_coords = false;
 242
 243         inverse_result res_b1_b2, res_b1_a1, res_b1_a2;
 244         if (! b_is_point)
 245         {
 246             res_b1_b2 = inverse_dist_azi::apply(b1_lon, b1_lat, b2_lon, b2_lat, spheroid);
 247             if (math::equals(res_b1_b2.distance, c0))
 248             {
 249                 b_is_point = true;
 250                 degen_neq_coords = true;
 251             }
 252             else
 253             {
 254                 res_b1_a1 = inverse_dist_azi::apply(b1_lon, b1_lat, a1_lon, a1_lat, spheroid);
 255                 if (math::equals(res_b1_a1.distance, c0))
 256                 {
 257                     is_equal_a1_b1 = true;
 258                 }
 259                 res_b1_a2 = inverse_dist_azi::apply(b1_lon, b1_lat, a2_lon, a2_lat, spheroid);
 260                 if (math::equals(res_b1_a2.distance, c0))
 261                 {
 262                     is_equal_a2_b1 = true;
 263                 }
 264                 sides.set<0>(is_equal_a1_b1 ? 0 : formula::azimuth_side_value(res_b1_a1.azimuth, res_b1_b2.azimuth),
 265                              is_equal_a2_b1 ? 0 : formula::azimuth_side_value(res_b1_a2.azimuth, res_b1_b2.azimuth));
 266                 if (sides.same<0>())
 267                 {
 268                     // Both points are at the same side of other segment, we can leave
 269                     return Policy::disjoint();
 270                 }
 271             }
 272         }
 273
 274         bool is_equal_a1_b2 = equals_point_point(a1, b2);
 275
 276         inverse_result res_a1_a2, res_a1_b1, res_a1_b2;
 277         if (! a_is_point)
 278         {
 279             res_a1_a2 = inverse_dist_azi::apply(a1_lon, a1_lat, a2_lon, a2_lat, spheroid);
 280             if (math::equals(res_a1_a2.distance, c0))
 281             {
 282                 a_is_point = true;
 283                 degen_neq_coords = true;
 284             }
 285             else
 286             {
 287                 res_a1_b1 = inverse_dist_azi::apply(a1_lon, a1_lat, b1_lon, b1_lat, spheroid);
 288                 if (math::equals(res_a1_b1.distance, c0))
 289                 {
 290                     is_equal_a1_b1 = true;
 291                 }
 292                 res_a1_b2 = inverse_dist_azi::apply(a1_lon, a1_lat, b2_lon, b2_lat, spheroid);
 293                 if (math::equals(res_a1_b2.distance, c0))
 294                 {
 295                     is_equal_a1_b2 = true;
 296                 }
 297                 sides.set<1>(is_equal_a1_b1 ? 0 : formula::azimuth_side_value(res_a1_b1.azimuth, res_a1_a2.azimuth),
 298                              is_equal_a1_b2 ? 0 : formula::azimuth_side_value(res_a1_b2.azimuth, res_a1_a2.azimuth));
 299                 if (sides.same<1>())
 300                 {
 301                     // Both points are at the same side of other segment, we can leave
 302                     return Policy::disjoint();
 303                 }
 304             }
 305         }
 306
 307         if(a_is_point && b_is_point)
 308         {
 309             return is_equal_a1_b2
 310                 ? Policy::degenerate(a, true)
 311                 : Policy::disjoint()
 312                 ;
 313         }
 314
 315         // NOTE: at this point the segments may still be disjoint
 316         // NOTE: at this point one of the segments may be degenerated
 317
 318         bool collinear = sides.collinear();
 319
 320         if (! collinear)
 321         {
 322             // WARNING: the side strategy doesn't have the info about the other
 323             // segment so it may return results inconsistent with this intersection
 324             // strategy, as it checks both segments for consistency
 325
 326             if (sides.get<0, 0>() == 0 && sides.get<0, 1>() == 0)
 327             {
 328                 collinear = true;
 329                 sides.set<1>(0, 0);
 330             }
 331             else if (sides.get<1, 0>() == 0 && sides.get<1, 1>() == 0)
 332             {
 333                 collinear = true;
 334                 sides.set<0>(0, 0);
 335             }
 336         }
 337
 338         if (collinear)
 339         {
 340             if (a_is_point)
 341             {
 342                 return collinear_one_degenerated<Policy, calc_t>(a, true, b1, b2, a1, a2, res_b1_b2, res_b1_a1, res_b1_a2, is_b_reversed, degen_neq_coords);
 343             }
 344             else if (b_is_point)
 345             {
 346                 return collinear_one_degenerated<Policy, calc_t>(b, false, a1, a2, b1, b2, res_a1_a2, res_a1_b1, res_a1_b2, is_a_reversed, degen_neq_coords);
 347             }
 348             else
 349             {
 350                 calc_t dist_a1_a2, dist_a1_b1, dist_a1_b2;
 351                 calc_t dist_b1_b2, dist_b1_a1, dist_b1_a2;
 352                 // use shorter segment
 353                 if (res_a1_a2.distance <= res_b1_b2.distance)
 354                 {
 355                     calculate_collinear_data(a1, a2, b1, b2, res_a1_a2, res_a1_b1, res_a1_b2, dist_a1_a2, dist_a1_b1);
 356                     calculate_collinear_data(a1, a2, b2, b1, res_a1_a2, res_a1_b2, res_a1_b1, dist_a1_a2, dist_a1_b2);
 357                     dist_b1_b2 = dist_a1_b2 - dist_a1_b1;
 358                     dist_b1_a1 = -dist_a1_b1;
 359                     dist_b1_a2 = dist_a1_a2 - dist_a1_b1;
 360                 }
 361                 else
 362                 {
 363                     calculate_collinear_data(b1, b2, a1, a2, res_b1_b2, res_b1_a1, res_b1_a2, dist_b1_b2, dist_b1_a1);
 364                     calculate_collinear_data(b1, b2, a2, a1, res_b1_b2, res_b1_a2, res_b1_a1, dist_b1_b2, dist_b1_a2);
 365                     dist_a1_a2 = dist_b1_a2 - dist_b1_a1;
 366                     dist_a1_b1 = -dist_b1_a1;
 367                     dist_a1_b2 = dist_b1_b2 - dist_b1_a1;
 368                 }
 369
 370                 // NOTE: this is probably not needed
 371                 int a1_on_b = position_value(c0, dist_a1_b1, dist_a1_b2);
 372                 int a2_on_b = position_value(dist_a1_a2, dist_a1_b1, dist_a1_b2);
 373                 int b1_on_a = position_value(c0, dist_b1_a1, dist_b1_a2);
 374                 int b2_on_a = position_value(dist_b1_b2, dist_b1_a1, dist_b1_a2);
 375
 376                 if ((a1_on_b < 1 && a2_on_b < 1) || (a1_on_b > 3 && a2_on_b > 3))
 377                 {
 378                     return Policy::disjoint();
 379                 }
 380
 381                 if (a1_on_b == 1)
 382                 {
 383                     dist_b1_a1 = 0;
 384                     dist_a1_b1 = 0;
 385                 }
 386                 else if (a1_on_b == 3)
 387                 {
 388                     dist_b1_a1 = dist_b1_b2;
 389                     dist_a1_b2 = 0;
 390                 }
 391
 392                 if (a2_on_b == 1)
 393                 {
 394                     dist_b1_a2 = 0;
 395                     dist_a1_b1 = dist_a1_a2;
 396                 }
 397                 else if (a2_on_b == 3)
 398                 {
 399                     dist_b1_a2 = dist_b1_b2;
 400                     dist_a1_b2 = dist_a1_a2;
 401                 }
 402
 403                 bool opposite = ! same_direction(res_a1_a2.azimuth, res_b1_b2.azimuth);
 404
 405                 // NOTE: If segment was reversed opposite, positions and segment ratios has to be altered
 406                 if (is_a_reversed)
 407                 {
 408                     // opposite
 409                     opposite = ! opposite;
 410                     // positions
 411                     std::swap(a1_on_b, a2_on_b);
 412                     b1_on_a = 4 - b1_on_a;
 413                     b2_on_a = 4 - b2_on_a;
 414                     // distances for ratios
 415                     std::swap(dist_b1_a1, dist_b1_a2);
 416                     dist_a1_b1 = dist_a1_a2 - dist_a1_b1;
 417                     dist_a1_b2 = dist_a1_a2 - dist_a1_b2;
 418                 }
 419                 if (is_b_reversed)
 420                 {
 421                     // opposite
 422                     opposite = ! opposite;
 423                     // positions
 424                     a1_on_b = 4 - a1_on_b;
 425                     a2_on_b = 4 - a2_on_b;
 426                     std::swap(b1_on_a, b2_on_a);
 427                     // distances for ratios
 428                     dist_b1_a1 = dist_b1_b2 - dist_b1_a1;
 429                     dist_b1_a2 = dist_b1_b2 - dist_b1_a2;
 430                     std::swap(dist_a1_b1, dist_a1_b2);
 431                 }
 432
 433                 segment_ratio<calc_t> ra_from(dist_b1_a1, dist_b1_b2);
 434                 segment_ratio<calc_t> ra_to(dist_b1_a2, dist_b1_b2);
 435                 segment_ratio<calc_t> rb_from(dist_a1_b1, dist_a1_a2);
 436                 segment_ratio<calc_t> rb_to(dist_a1_b2, dist_a1_a2);
 437
 438                 return Policy::segments_collinear(a, b, opposite,
 439                     a1_on_b, a2_on_b, b1_on_a, b2_on_a,
 440                     ra_from, ra_to, rb_from, rb_to);
 441             }
 442         }
 443         else // crossing or touching
 444         {
 445             if (a_is_point || b_is_point)
 446             {
 447                 return Policy::disjoint();
 448             }
 449
 450             calc_t lon = 0, lat = 0;
 451             intersection_point_flag ip_flag;
 452             calc_t dist_a1_a2, dist_a1_i1, dist_b1_b2, dist_b1_i1;
 453             if (calculate_ip_data(a1, a2, b1, b2,
 454                                   a1_lon, a1_lat, a2_lon, a2_lat,
 455                                   b1_lon, b1_lat, b2_lon, b2_lat,
 456                                   res_a1_a2, res_a1_b1, res_a1_b2,
 457                                   res_b1_b2, res_b1_a1, res_b1_a2,
 458                                   sides, spheroid,
 459                                   lon, lat,
 460                                   dist_a1_a2, dist_a1_i1, dist_b1_b2, dist_b1_i1,
 461                                   ip_flag))
 462             {
 463                 // NOTE: If segment was reversed sides and segment ratios has to be altered
 464                 if (is_a_reversed)
 465                 {
 466                     // sides
 467                     sides_reverse_segment<0>(sides);
 468                     // distance for ratio
 469                     dist_a1_i1 = dist_a1_a2 - dist_a1_i1;
 470                     // ip flag
 471                     ip_flag_reverse_segment(ip_flag, ipi_at_a1, ipi_at_a2);
 472                 }
 473                 if (is_b_reversed)
 474                 {
 475                     // sides
 476                     sides_reverse_segment<1>(sides);
 477                     // distance for ratio
 478                     dist_b1_i1 = dist_b1_b2 - dist_b1_i1;
 479                     // ip flag
 480                     ip_flag_reverse_segment(ip_flag, ipi_at_b1, ipi_at_b2);
 481                 }
 482
 483                 // intersects
 484                 segment_intersection_info
 485                     <
 486                         calc_t,
 487                         segment_ratio<calc_t>
 488                     > sinfo;
 489
 490                 sinfo.lon = lon;
 491                 sinfo.lat = lat;
 492                 sinfo.robust_ra.assign(dist_a1_i1, dist_a1_a2);
 493                 sinfo.robust_rb.assign(dist_b1_i1, dist_b1_b2);
 494                 sinfo.ip_flag = ip_flag;
 495
 496                 return Policy::segments_crosses(sides, sinfo, a, b);
 497             }
 498             else
 499             {
 500                 return Policy::disjoint();
 501             }
 502         }
 503     }
 504
 505     template <typename Policy, typename CalcT, typename Segment, typename Point1, typename Point2, typename ResultInverse>
 506     static inline typename Policy::return_type
 507         collinear_one_degenerated(Segment const& segment, bool degenerated_a,
 508                                   Point1 const& a1, Point1 const& a2,
 509                                   Point2 const& b1, Point2 const& b2,
 510                                   ResultInverse const& res_a1_a2,
 511                                   ResultInverse const& res_a1_b1,
 512                                   ResultInverse const& res_a1_b2,
 513                                   bool is_other_reversed,
 514                                   bool degen_neq_coords)
 515     {
 516         CalcT dist_1_2, dist_1_o;
 517         if (! calculate_collinear_data(a1, a2, b1, b2, res_a1_a2, res_a1_b1, res_a1_b2, dist_1_2, dist_1_o, degen_neq_coords))
 518         {
 519             return Policy::disjoint();
 520         }
 521
 522         // NOTE: If segment was reversed segment ratio has to be altered
 523         if (is_other_reversed)
 524         {
 525             // distance for ratio
 526             dist_1_o = dist_1_2 - dist_1_o;
 527         }
 528
 529         return Policy::one_degenerate(segment, segment_ratio<CalcT>(dist_1_o, dist_1_2), degenerated_a);
 530     }
 531
 532     // TODO: instead of checks below test bi against a1 and a2 here?
 533     //       in order to make this independent from is_near()
 534     template <typename Point1, typename Point2, typename ResultInverse, typename CalcT>
 535     static inline bool calculate_collinear_data(Point1 const& a1, Point1 const& a2, // in
 536                                                 Point2 const& b1, Point2 const& /*b2*/, // in
 537                                                 ResultInverse const& res_a1_a2,     // in
 538                                                 ResultInverse const& res_a1_b1,     // in
 539                                                 ResultInverse const& res_a1_b2,     // in
 540                                                 CalcT& dist_a1_a2,                  // out
 541                                                 CalcT& dist_a1_b1,                  // out
 542                                                 bool degen_neq_coords = false)      // in
 543     {
 544         dist_a1_a2 = res_a1_a2.distance;
 545
 546         dist_a1_b1 = res_a1_b1.distance;
 547         if (! same_direction(res_a1_b1.azimuth, res_a1_a2.azimuth))
 548         {
 549             dist_a1_b1 = -dist_a1_b1;
 550         }
 551
 552         // if b1 is close a1
 553         if (is_endpoint_equal(dist_a1_b1, a1, b1))
 554         {
 555             dist_a1_b1 = 0;
 556             return true;
 557         }
 558         // if b1 is close a2
 559         else if (is_endpoint_equal(dist_a1_a2 - dist_a1_b1, a2, b1))
 560         {
 561             dist_a1_b1 = dist_a1_a2;
 562             return true;
 563         }
 564
 565         // check the other endpoint of degenerated segment near a pole
 566         if (degen_neq_coords)
 567         {
 568             static CalcT const c0 = 0;
 569             if (math::equals(res_a1_b2.distance, c0))
 570             {
 571                 dist_a1_b1 = 0;
 572                 return true;
 573             }
 574             else if (math::equals(dist_a1_a2 - res_a1_b2.distance, c0))
 575             {
 576                 dist_a1_b1 = dist_a1_a2;
 577                 return true;
 578             }
 579         }
 580
 581         // or i1 is on b
 582         return segment_ratio<CalcT>(dist_a1_b1, dist_a1_a2).on_segment();
 583     }
 584
 585     template <typename Point1, typename Point2, typename CalcT, typename ResultInverse, typename Spheroid_>
 586     static inline bool calculate_ip_data(Point1 const& a1, Point1 const& a2,       // in
 587                                          Point2 const& b1, Point2 const& b2,       // in
 588                                          CalcT const& a1_lon, CalcT const& a1_lat, // in
 589                                          CalcT const& a2_lon, CalcT const& a2_lat, // in
 590                                          CalcT const& b1_lon, CalcT const& b1_lat, // in
 591                                          CalcT const& b2_lon, CalcT const& b2_lat, // in
 592                                          ResultInverse const& res_a1_a2,           // in
 593                                          ResultInverse const& res_a1_b1,           // in
 594                                          ResultInverse const& res_a1_b2,           // in
 595                                          ResultInverse const& res_b1_b2,           // in
 596                                          ResultInverse const& res_b1_a1,           // in
 597                                          ResultInverse const& res_b1_a2,           // in
 598                                          side_info const& sides,                   // in
 599                                          Spheroid_ const& spheroid,                // in
 600                                          CalcT & lon, CalcT & lat,             // out
 601                                          CalcT& dist_a1_a2, CalcT& dist_a1_ip, // out
 602                                          CalcT& dist_b1_b2, CalcT& dist_b1_ip, // out
 603                                          intersection_point_flag& ip_flag)     // out
 604     {
 605         dist_a1_a2 = res_a1_a2.distance;
 606         dist_b1_b2 = res_b1_b2.distance;
 607
 608         // assign the IP if some endpoints overlap
 609         if (equals_point_point(a1, b1))
 610         {
 611             lon = a1_lon;
 612             lat = a1_lat;
 613             dist_a1_ip = 0;
 614             dist_b1_ip = 0;
 615             ip_flag = ipi_at_a1;
 616             return true;
 617         }
 618         else if (equals_point_point(a1, b2))
 619         {
 620             lon = a1_lon;
 621             lat = a1_lat;
 622             dist_a1_ip = 0;
 623             dist_b1_ip = dist_b1_b2;
 624             ip_flag = ipi_at_a1;
 625             return true;
 626         }
 627         else if (equals_point_point(a2, b1))
 628         {
 629             lon = a2_lon;
 630             lat = a2_lat;
 631             dist_a1_ip = dist_a1_a2;
 632             dist_b1_ip = 0;
 633             ip_flag = ipi_at_a2;
 634             return true;
 635         }
 636         else if (equals_point_point(a2, b2))
 637         {
 638             lon = a2_lon;
 639             lat = a2_lat;
 640             dist_a1_ip = dist_a1_a2;
 641             dist_b1_ip = dist_b1_b2;
 642             ip_flag = ipi_at_a2;
 643             return true;
 644         }
 645
 646         // at this point we know that the endpoints doesn't overlap
 647         // check cases when an endpoint lies on the other geodesic
 648         if (sides.template get<0, 0>() == 0) // a1 wrt b
 649         {
 650             if (res_b1_a1.distance <= res_b1_b2.distance
 651                 && same_direction(res_b1_a1.azimuth, res_b1_b2.azimuth))
 652             {
 653                 lon = a1_lon;
 654                 lat = a1_lat;
 655                 dist_a1_ip = 0;
 656                 dist_b1_ip = res_b1_a1.distance;
 657                 ip_flag = ipi_at_a1;
 658                 return true;
 659             }
 660             else
 661             {
 662                 return false;
 663             }
 664         }
 665         else if (sides.template get<0, 1>() == 0) // a2 wrt b
 666         {
 667             if (res_b1_a2.distance <= res_b1_b2.distance
 668                 && same_direction(res_b1_a2.azimuth, res_b1_b2.azimuth))
 669             {
 670                 lon = a2_lon;
 671                 lat = a2_lat;
 672                 dist_a1_ip = res_a1_a2.distance;
 673                 dist_b1_ip = res_b1_a2.distance;
 674                 ip_flag = ipi_at_a2;
 675                 return true;
 676             }
 677             else
 678             {
 679                 return false;
 680             }
 681         }
 682         else if (sides.template get<1, 0>() == 0) // b1 wrt a
 683         {
 684             if (res_a1_b1.distance <= res_a1_a2.distance
 685                 && same_direction(res_a1_b1.azimuth, res_a1_a2.azimuth))
 686             {
 687                 lon = b1_lon;
 688                 lat = b1_lat;
 689                 dist_a1_ip = res_a1_b1.distance;
 690                 dist_b1_ip = 0;
 691                 ip_flag = ipi_at_b1;
 692                 return true;
 693             }
 694             else
 695             {
 696                 return false;
 697             }
 698         }
 699         else if (sides.template get<1, 1>() == 0) // b2 wrt a
 700         {
 701             if (res_a1_b2.distance <= res_a1_a2.distance
 702                 && same_direction(res_a1_b2.azimuth, res_a1_a2.azimuth))
 703             {
 704                 lon = b2_lon;
 705                 lat = b2_lat;
 706                 dist_a1_ip = res_a1_b2.distance;
 707                 dist_b1_ip = res_b1_b2.distance;
 708                 ip_flag = ipi_at_b2;
 709                 return true;
 710             }
 711             else
 712             {
 713                 return false;
 714             }
 715         }
 716
 717         // At this point neither the endpoints overlaps
 718         // nor any andpoint lies on the other geodesic
 719         // So the endpoints should lie on the opposite sides of both geodesics
 720
 721         bool const ok = formula::sjoberg_intersection<CalcT, FormulaPolicy::template inverse, Order>
 722                         ::apply(a1_lon, a1_lat, a2_lon, a2_lat, res_a1_a2.azimuth,
 723                                 b1_lon, b1_lat, b2_lon, b2_lat, res_b1_b2.azimuth,
 724                                 lon, lat, spheroid);
 725
 726         if (! ok)
 727         {
 728             return false;
 729         }
 730
 731         typedef typename FormulaPolicy::template inverse<CalcT, true, true, false, false, false> inverse_dist_azi;
 732         typedef typename inverse_dist_azi::result_type inverse_result;
 733
 734         inverse_result const res_a1_ip = inverse_dist_azi::apply(a1_lon, a1_lat, lon, lat, spheroid);
 735         dist_a1_ip = res_a1_ip.distance;
 736         if (! same_direction(res_a1_ip.azimuth, res_a1_a2.azimuth))
 737         {
 738             dist_a1_ip = -dist_a1_ip;
 739         }
 740
 741         bool is_on_a = segment_ratio<CalcT>(dist_a1_ip, dist_a1_a2).on_segment();
 742         // NOTE: not fully consistent with equals_point_point() since radians are always used.
 743         bool is_on_a1 = math::equals(lon, a1_lon) && math::equals(lat, a1_lat);
 744         bool is_on_a2 = math::equals(lon, a2_lon) && math::equals(lat, a2_lat);
 745
 746         if (! (is_on_a || is_on_a1 || is_on_a2))
 747         {
 748             return false;
 749         }
 750
 751         inverse_result const res_b1_ip = inverse_dist_azi::apply(b1_lon, b1_lat, lon, lat, spheroid);
 752         dist_b1_ip = res_b1_ip.distance;
 753         if (! same_direction(res_b1_ip.azimuth, res_b1_b2.azimuth))
 754         {
 755             dist_b1_ip = -dist_b1_ip;
 756         }
 757
 758         bool is_on_b = segment_ratio<CalcT>(dist_b1_ip, dist_b1_b2).on_segment();
 759         // NOTE: not fully consistent with equals_point_point() since radians are always used.
 760         bool is_on_b1 = math::equals(lon, b1_lon) && math::equals(lat, b1_lat);
 761         bool is_on_b2 = math::equals(lon, b2_lon) && math::equals(lat, b2_lat);
 762
 763         if (! (is_on_b || is_on_b1 || is_on_b2))
 764         {
 765             return false;
 766         }
 767
 768         typedef typename FormulaPolicy::template inverse<CalcT, true, false, false, false, false> inverse_dist;
 769
 770         ip_flag = ipi_inters;
 771
 772         if (is_on_b1)
 773         {
 774             lon = b1_lon;
 775             lat = b1_lat;
 776             dist_a1_ip = inverse_dist::apply(a1_lon, a1_lat, lon, lat, spheroid).distance; // for consistency
 777             dist_b1_ip = 0;
 778             ip_flag = ipi_at_b1;
 779         }
 780         else if (is_on_b2)
 781         {
 782             lon = b2_lon;
 783             lat = b2_lat;
 784             dist_a1_ip = inverse_dist::apply(a1_lon, a1_lat, lon, lat, spheroid).distance; // for consistency
 785             dist_b1_ip = res_b1_b2.distance;
 786             ip_flag = ipi_at_b2;
 787         }
 788
 789         if (is_on_a1)
 790         {
 791             lon = a1_lon;
 792             lat = a1_lat;
 793             dist_a1_ip = 0;
 794             dist_b1_ip = inverse_dist::apply(b1_lon, b1_lat, lon, lat, spheroid).distance; // for consistency
 795             ip_flag = ipi_at_a1;
 796         }
 797         else if (is_on_a2)
 798         {
 799             lon = a2_lon;
 800             lat = a2_lat;
 801             dist_a1_ip = res_a1_a2.distance;
 802             dist_b1_ip = inverse_dist::apply(b1_lon, b1_lat, lon, lat, spheroid).distance; // for consistency
 803             ip_flag = ipi_at_a2;
 804         }
 805
 806         return true;
 807     }
 808
 809     template <typename CalcT, typename P1, typename P2>
 810     static inline bool is_endpoint_equal(CalcT const& dist,
 811                                          P1 const& ai, P2 const& b1)
 812     {
 813         static CalcT const c0 = 0;
 814         return is_near(dist) && (math::equals(dist, c0) || equals_point_point(ai, b1));
 815     }
 816
 817     template <typename CalcT>
 818     static inline bool is_near(CalcT const& dist)
 819     {
 820         // NOTE: This strongly depends on the Inverse method
 821         CalcT const small_number = CalcT(std::is_same<CalcT, float>::value ? 0.0001 : 0.00000001);
 822         return math::abs(dist) <= small_number;
 823     }
 824
 825     template <typename ProjCoord1, typename ProjCoord2>
 826     static inline int position_value(ProjCoord1 const& ca1,
 827                                      ProjCoord2 const& cb1,
 828                                      ProjCoord2 const& cb2)
 829     {
 830         // S1x  0   1    2     3   4
 831         // S2       |---------->
 832         return math::equals(ca1, cb1) ? 1
 833              : math::equals(ca1, cb2) ? 3
 834              : cb1 < cb2 ?
 835                 ( ca1 < cb1 ? 0
 836                 : ca1 > cb2 ? 4
 837                 : 2 )
 838               : ( ca1 > cb1 ? 0
 839                 : ca1 < cb2 ? 4
 840                 : 2 );
 841     }
 842
 843     template <typename CalcT>
 844     static inline bool same_direction(CalcT const& azimuth1, CalcT const& azimuth2)
 845     {
 846         // distance between two angles normalized to (-180, 180]
 847         CalcT const angle_diff = math::longitude_distance_signed<radian>(azimuth1, azimuth2);
 848         return math::abs(angle_diff) <= math::half_pi<CalcT>();
 849     }
 850
 851     template <int Which>
 852     static inline void sides_reverse_segment(side_info & sides)
 853     {
 854         // names assuming segment A is reversed (Which == 0)
 855         int a1_wrt_b = sides.template get<Which, 0>();
 856         int a2_wrt_b = sides.template get<Which, 1>();
 857         std::swap(a1_wrt_b, a2_wrt_b);
 858         sides.template set<Which>(a1_wrt_b, a2_wrt_b);
 859         int b1_wrt_a = sides.template get<1 - Which, 0>();
 860         int b2_wrt_a = sides.template get<1 - Which, 1>();
 861         sides.template set<1 - Which>(-b1_wrt_a, -b2_wrt_a);
 862     }
 863
 864     static inline void ip_flag_reverse_segment(intersection_point_flag & ip_flag,
 865                                                intersection_point_flag const& ipi_at_p1,
 866                                                intersection_point_flag const& ipi_at_p2)
 867     {
 868         ip_flag = ip_flag == ipi_at_p1 ? ipi_at_p2 :
 869                   ip_flag == ipi_at_p2 ? ipi_at_p1 :
 870                   ip_flag;
 871     }
 872
 873     template <typename Point1, typename Point2>
 874     static inline bool equals_point_point(Point1 const& point1, Point2 const& point2)
 875     {
 876         return strategy::within::spherical_point_point::apply(point1, point2);
 877     }
 878
 879 private:
 880     Spheroid m_spheroid;
 881 };
 882
 883
 884 }} // namespace strategy::intersection
 885
 886 }} // namespace boost::geometry
 887
 888
 889 #endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_INTERSECTION_HPP