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-2020, 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     unsigned int 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     typedef side::geographic
  84         <
  85             FormulaPolicy, Spheroid, CalculationType
  86         > side_strategy_type;
  87
  88     inline side_strategy_type get_side_strategy() const
  89     {
  90         return side_strategy_type(m_spheroid);
  91     }
  92
  93     template <typename Geometry1, typename Geometry2>
  94     struct point_in_geometry_strategy
  95     {
  96         typedef strategy::within::geographic_winding
  97             <
  98                 typename point_type<Geometry1>::type,
  99                 typename point_type<Geometry2>::type,
 100                 FormulaPolicy,
 101                 Spheroid,
 102                 CalculationType
 103             > type;
 104     };
 105
 106     template <typename Geometry1, typename Geometry2>
 107     inline typename point_in_geometry_strategy<Geometry1, Geometry2>::type
 108         get_point_in_geometry_strategy() const
 109     {
 110         typedef typename point_in_geometry_strategy
 111             <
 112                 Geometry1, Geometry2
 113             >::type strategy_type;
 114         return strategy_type(m_spheroid);
 115     }
 116
 117     template <typename Geometry>
 118     struct area_strategy
 119     {
 120         typedef area::geographic
 121             <
 122                 FormulaPolicy,
 123                 Order,
 124                 Spheroid,
 125                 CalculationType
 126             > type;
 127     };
 128
 129     template <typename Geometry>
 130     inline typename area_strategy<Geometry>::type get_area_strategy() const
 131     {
 132         typedef typename area_strategy<Geometry>::type strategy_type;
 133         return strategy_type(m_spheroid);
 134     }
 135
 136     template <typename Geometry>
 137     struct distance_strategy
 138     {
 139         typedef distance::geographic
 140             <
 141                 FormulaPolicy,
 142                 Spheroid,
 143                 CalculationType
 144             > type;
 145     };
 146
 147     template <typename Geometry>
 148     inline typename distance_strategy<Geometry>::type get_distance_strategy() const
 149     {
 150         typedef typename distance_strategy<Geometry>::type strategy_type;
 151         return strategy_type(m_spheroid);
 152     }
 153
 154     typedef envelope::geographic<FormulaPolicy, Spheroid, CalculationType>
 155         envelope_strategy_type;
 156
 157     inline envelope_strategy_type get_envelope_strategy() const
 158     {
 159         return envelope_strategy_type(m_spheroid);
 160     }
 161
 162     typedef expand::geographic_segment<FormulaPolicy, Spheroid, CalculationType>
 163         expand_strategy_type;
 164
 165     inline expand_strategy_type get_expand_strategy() const
 166     {
 167         return expand_strategy_type(m_spheroid);
 168     }
 169
 170     typedef within::spherical_point_point point_in_point_strategy_type;
 171
 172     static inline point_in_point_strategy_type get_point_in_point_strategy()
 173     {
 174         return point_in_point_strategy_type();
 175     }
 176
 177     typedef within::spherical_point_point equals_point_point_strategy_type;
 178
 179     static inline equals_point_point_strategy_type get_equals_point_point_strategy()
 180     {
 181         return equals_point_point_strategy_type();
 182     }
 183
 184     typedef disjoint::spherical_box_box disjoint_box_box_strategy_type;
 185
 186     static inline disjoint_box_box_strategy_type get_disjoint_box_box_strategy()
 187     {
 188         return disjoint_box_box_strategy_type();
 189     }
 190
 191     typedef disjoint::segment_box_geographic
 192         <
 193             FormulaPolicy, Spheroid, CalculationType
 194         > disjoint_segment_box_strategy_type;
 195
 196     inline disjoint_segment_box_strategy_type get_disjoint_segment_box_strategy() const
 197     {
 198         return disjoint_segment_box_strategy_type(m_spheroid);
 199     }
 200
 201     typedef covered_by::spherical_point_box disjoint_point_box_strategy_type;
 202     typedef covered_by::spherical_point_box covered_by_point_box_strategy_type;
 203     typedef within::spherical_point_box within_point_box_strategy_type;
 204     typedef envelope::spherical_box envelope_box_strategy_type;
 205     typedef expand::spherical_box expand_box_strategy_type;
 206
 207     enum intersection_point_flag { ipi_inters = 0, ipi_at_a1, ipi_at_a2, ipi_at_b1, ipi_at_b2 };
 208
 209     template <typename CoordinateType, typename SegmentRatio>
 210     struct segment_intersection_info
 211     {
 212         template <typename Point, typename Segment1, typename Segment2>
 213         void calculate(Point& point, Segment1 const& a, Segment2 const& b) const
 214         {
 215             if (ip_flag == ipi_inters)
 216             {
 217                 // TODO: assign the rest of coordinates
 218                 set_from_radian<0>(point, lon);
 219                 set_from_radian<1>(point, lat);
 220             }
 221             else if (ip_flag == ipi_at_a1)
 222             {
 223                 detail::assign_point_from_index<0>(a, point);
 224             }
 225             else if (ip_flag == ipi_at_a2)
 226             {
 227                 detail::assign_point_from_index<1>(a, point);
 228             }
 229             else if (ip_flag == ipi_at_b1)
 230             {
 231                 detail::assign_point_from_index<0>(b, point);
 232             }
 233             else // ip_flag == ipi_at_b2
 234             {
 235                 detail::assign_point_from_index<1>(b, point);
 236             }
 237         }
 238
 239         CoordinateType lon;
 240         CoordinateType lat;
 241         SegmentRatio robust_ra;
 242         SegmentRatio robust_rb;
 243         intersection_point_flag ip_flag;
 244     };
 245
 246     explicit geographic_segments(Spheroid const& spheroid = Spheroid())
 247         : m_spheroid(spheroid)
 248     {}
 249
 250     // Relate segments a and b
 251     template
 252     <
 253         typename UniqueSubRange1,
 254         typename UniqueSubRange2,
 255         typename Policy
 256     >
 257     inline typename Policy::return_type apply(UniqueSubRange1 const& range_p,
 258                                               UniqueSubRange2 const& range_q,
 259                                               Policy const&) const
 260     {
 261         typedef typename UniqueSubRange1::point_type point1_type;
 262         typedef typename UniqueSubRange2::point_type point2_type;
 263         typedef model::referring_segment<point1_type const> segment_type1;
 264         typedef model::referring_segment<point2_type const> segment_type2;
 265
 266         BOOST_CONCEPT_ASSERT( (concepts::ConstPoint<point1_type>) );
 267         BOOST_CONCEPT_ASSERT( (concepts::ConstPoint<point2_type>) );
 268
 269         /*
 270         typename coordinate_type<Point1>::type
 271             const a1_lon = get<0>(a1),
 272             const a2_lon = get<0>(a2);
 273         typename coordinate_type<Point2>::type
 274             const b1_lon = get<0>(b1),
 275             const b2_lon = get<0>(b2);
 276         bool is_a_reversed = a1_lon > a2_lon || a1_lon == a2_lon && get<1>(a1) > get<1>(a2);
 277         bool is_b_reversed = b1_lon > b2_lon || b1_lon == b2_lon && get<1>(b1) > get<1>(b2);
 278         */
 279
 280         point1_type const& p0 = range_p.at(0);
 281         point1_type const& p1 = range_p.at(1);
 282         point2_type const& q0 = range_q.at(0);
 283         point2_type const& q1 = range_q.at(1);
 284
 285         bool const is_p_reversed = get<1>(p0) > get<1>(p1);
 286         bool const is_q_reversed = get<1>(q0) > get<1>(q1);
 287
 288         // Call apply with original segments and ordered points
 289         return apply<Policy>(segment_type1(p0, p1),
 290                              segment_type2(q0, q1),
 291                              (is_p_reversed ? p1 : p0),
 292                              (is_p_reversed ? p0 : p1),
 293                              (is_q_reversed ? q1 : q0),
 294                              (is_q_reversed ? q0 : q1),
 295                              is_p_reversed, is_q_reversed);
 296     }
 297
 298 private:
 299     // Relate segments a and b
 300     template
 301     <
 302         typename Policy,
 303         typename Segment1,
 304         typename Segment2,
 305         typename Point1,
 306         typename Point2
 307     >
 308     inline typename Policy::return_type apply(Segment1 const& a, Segment2 const& b,
 309                                               Point1 const& a1, Point1 const& a2,
 310                                               Point2 const& b1, Point2 const& b2,
 311                                               bool is_a_reversed, bool is_b_reversed) const
 312     {
 313         BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment1>) );
 314         BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment2>) );
 315
 316         typedef typename select_calculation_type
 317             <Segment1, Segment2, CalculationType>::type calc_t;
 318
 319         typedef srs::spheroid<calc_t> spheroid_type;
 320
 321         static const calc_t c0 = 0;
 322
 323         // normalized spheroid
 324         spheroid_type spheroid = formula::unit_spheroid<spheroid_type>(m_spheroid);
 325
 326         // TODO: check only 2 first coordinates here?
 327         bool a_is_point = equals_point_point(a1, a2);
 328         bool b_is_point = equals_point_point(b1, b2);
 329
 330         if(a_is_point && b_is_point)
 331         {
 332             return equals_point_point(a1, b2)
 333                 ? Policy::degenerate(a, true)
 334                 : Policy::disjoint()
 335                 ;
 336         }
 337
 338         calc_t const a1_lon = get_as_radian<0>(a1);
 339         calc_t const a1_lat = get_as_radian<1>(a1);
 340         calc_t const a2_lon = get_as_radian<0>(a2);
 341         calc_t const a2_lat = get_as_radian<1>(a2);
 342         calc_t const b1_lon = get_as_radian<0>(b1);
 343         calc_t const b1_lat = get_as_radian<1>(b1);
 344         calc_t const b2_lon = get_as_radian<0>(b2);
 345         calc_t const b2_lat = get_as_radian<1>(b2);
 346
 347         side_info sides;
 348
 349         // NOTE: potential optimization, don't calculate distance at this point
 350         // this would require to reimplement inverse strategy to allow
 351         // calculation of distance if needed, probably also storing intermediate
 352         // results somehow inside an object.
 353         typedef typename FormulaPolicy::template inverse<calc_t, true, true, false, false, false> inverse_dist_azi;
 354         typedef typename inverse_dist_azi::result_type inverse_result;
 355
 356         // TODO: no need to call inverse formula if we know that the points are equal
 357         // distance can be set to 0 in this case and azimuth may be not calculated
 358         bool is_equal_a1_b1 = equals_point_point(a1, b1);
 359         bool is_equal_a2_b1 = equals_point_point(a2, b1);
 360         bool degen_neq_coords = false;
 361
 362         inverse_result res_b1_b2, res_b1_a1, res_b1_a2;
 363         if (! b_is_point)
 364         {
 365             res_b1_b2 = inverse_dist_azi::apply(b1_lon, b1_lat, b2_lon, b2_lat, spheroid);
 366             if (math::equals(res_b1_b2.distance, c0))
 367             {
 368                 b_is_point = true;
 369                 degen_neq_coords = true;
 370             }
 371             else
 372             {
 373                 res_b1_a1 = inverse_dist_azi::apply(b1_lon, b1_lat, a1_lon, a1_lat, spheroid);
 374                 if (math::equals(res_b1_a1.distance, c0))
 375                 {
 376                     is_equal_a1_b1 = true;
 377                 }
 378                 res_b1_a2 = inverse_dist_azi::apply(b1_lon, b1_lat, a2_lon, a2_lat, spheroid);
 379                 if (math::equals(res_b1_a2.distance, c0))
 380                 {
 381                     is_equal_a2_b1 = true;
 382                 }
 383                 sides.set<0>(is_equal_a1_b1 ? 0 : formula::azimuth_side_value(res_b1_a1.azimuth, res_b1_b2.azimuth),
 384                              is_equal_a2_b1 ? 0 : formula::azimuth_side_value(res_b1_a2.azimuth, res_b1_b2.azimuth));
 385                 if (sides.same<0>())
 386                 {
 387                     // Both points are at the same side of other segment, we can leave
 388                     return Policy::disjoint();
 389                 }
 390             }
 391         }
 392
 393         bool is_equal_a1_b2 = equals_point_point(a1, b2);
 394
 395         inverse_result res_a1_a2, res_a1_b1, res_a1_b2;
 396         if (! a_is_point)
 397         {
 398             res_a1_a2 = inverse_dist_azi::apply(a1_lon, a1_lat, a2_lon, a2_lat, spheroid);
 399             if (math::equals(res_a1_a2.distance, c0))
 400             {
 401                 a_is_point = true;
 402                 degen_neq_coords = true;
 403             }
 404             else
 405             {
 406                 res_a1_b1 = inverse_dist_azi::apply(a1_lon, a1_lat, b1_lon, b1_lat, spheroid);
 407                 if (math::equals(res_a1_b1.distance, c0))
 408                 {
 409                     is_equal_a1_b1 = true;
 410                 }
 411                 res_a1_b2 = inverse_dist_azi::apply(a1_lon, a1_lat, b2_lon, b2_lat, spheroid);
 412                 if (math::equals(res_a1_b2.distance, c0))
 413                 {
 414                     is_equal_a1_b2 = true;
 415                 }
 416                 sides.set<1>(is_equal_a1_b1 ? 0 : formula::azimuth_side_value(res_a1_b1.azimuth, res_a1_a2.azimuth),
 417                              is_equal_a1_b2 ? 0 : formula::azimuth_side_value(res_a1_b2.azimuth, res_a1_a2.azimuth));
 418                 if (sides.same<1>())
 419                 {
 420                     // Both points are at the same side of other segment, we can leave
 421                     return Policy::disjoint();
 422                 }
 423             }
 424         }
 425
 426         if(a_is_point && b_is_point)
 427         {
 428             return is_equal_a1_b2
 429                 ? Policy::degenerate(a, true)
 430                 : Policy::disjoint()
 431                 ;
 432         }
 433
 434         // NOTE: at this point the segments may still be disjoint
 435         // NOTE: at this point one of the segments may be degenerated
 436
 437         bool collinear = sides.collinear();
 438
 439         if (! collinear)
 440         {
 441             // WARNING: the side strategy doesn't have the info about the other
 442             // segment so it may return results inconsistent with this intersection
 443             // strategy, as it checks both segments for consistency
 444
 445             if (sides.get<0, 0>() == 0 && sides.get<0, 1>() == 0)
 446             {
 447                 collinear = true;
 448                 sides.set<1>(0, 0);
 449             }
 450             else if (sides.get<1, 0>() == 0 && sides.get<1, 1>() == 0)
 451             {
 452                 collinear = true;
 453                 sides.set<0>(0, 0);
 454             }
 455         }
 456
 457         if (collinear)
 458         {
 459             if (a_is_point)
 460             {
 461                 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);
 462             }
 463             else if (b_is_point)
 464             {
 465                 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);
 466             }
 467             else
 468             {
 469                 calc_t dist_a1_a2, dist_a1_b1, dist_a1_b2;
 470                 calc_t dist_b1_b2, dist_b1_a1, dist_b1_a2;
 471                 // use shorter segment
 472                 if (res_a1_a2.distance <= res_b1_b2.distance)
 473                 {
 474                     calculate_collinear_data(a1, a2, b1, b2, res_a1_a2, res_a1_b1, res_a1_b2, dist_a1_a2, dist_a1_b1);
 475                     calculate_collinear_data(a1, a2, b2, b1, res_a1_a2, res_a1_b2, res_a1_b1, dist_a1_a2, dist_a1_b2);
 476                     dist_b1_b2 = dist_a1_b2 - dist_a1_b1;
 477                     dist_b1_a1 = -dist_a1_b1;
 478                     dist_b1_a2 = dist_a1_a2 - dist_a1_b1;
 479                 }
 480                 else
 481                 {
 482                     calculate_collinear_data(b1, b2, a1, a2, res_b1_b2, res_b1_a1, res_b1_a2, dist_b1_b2, dist_b1_a1);
 483                     calculate_collinear_data(b1, b2, a2, a1, res_b1_b2, res_b1_a2, res_b1_a1, dist_b1_b2, dist_b1_a2);
 484                     dist_a1_a2 = dist_b1_a2 - dist_b1_a1;
 485                     dist_a1_b1 = -dist_b1_a1;
 486                     dist_a1_b2 = dist_b1_b2 - dist_b1_a1;
 487                 }
 488
 489                 // NOTE: this is probably not needed
 490                 int a1_on_b = position_value(c0, dist_a1_b1, dist_a1_b2);
 491                 int a2_on_b = position_value(dist_a1_a2, dist_a1_b1, dist_a1_b2);
 492                 int b1_on_a = position_value(c0, dist_b1_a1, dist_b1_a2);
 493                 int b2_on_a = position_value(dist_b1_b2, dist_b1_a1, dist_b1_a2);
 494
 495                 if ((a1_on_b < 1 && a2_on_b < 1) || (a1_on_b > 3 && a2_on_b > 3))
 496                 {
 497                     return Policy::disjoint();
 498                 }
 499
 500                 if (a1_on_b == 1)
 501                 {
 502                     dist_b1_a1 = 0;
 503                     dist_a1_b1 = 0;
 504                 }
 505                 else if (a1_on_b == 3)
 506                 {
 507                     dist_b1_a1 = dist_b1_b2;
 508                     dist_a1_b2 = 0;
 509                 }
 510
 511                 if (a2_on_b == 1)
 512                 {
 513                     dist_b1_a2 = 0;
 514                     dist_a1_b1 = dist_a1_a2;
 515                 }
 516                 else if (a2_on_b == 3)
 517                 {
 518                     dist_b1_a2 = dist_b1_b2;
 519                     dist_a1_b2 = dist_a1_a2;
 520                 }
 521
 522                 bool opposite = ! same_direction(res_a1_a2.azimuth, res_b1_b2.azimuth);
 523
 524                 // NOTE: If segment was reversed opposite, positions and segment ratios has to be altered
 525                 if (is_a_reversed)
 526                 {
 527                     // opposite
 528                     opposite = ! opposite;
 529                     // positions
 530                     std::swap(a1_on_b, a2_on_b);
 531                     b1_on_a = 4 - b1_on_a;
 532                     b2_on_a = 4 - b2_on_a;
 533                     // distances for ratios
 534                     std::swap(dist_b1_a1, dist_b1_a2);
 535                     dist_a1_b1 = dist_a1_a2 - dist_a1_b1;
 536                     dist_a1_b2 = dist_a1_a2 - dist_a1_b2;
 537                 }
 538                 if (is_b_reversed)
 539                 {
 540                     // opposite
 541                     opposite = ! opposite;
 542                     // positions
 543                     a1_on_b = 4 - a1_on_b;
 544                     a2_on_b = 4 - a2_on_b;
 545                     std::swap(b1_on_a, b2_on_a);
 546                     // distances for ratios
 547                     dist_b1_a1 = dist_b1_b2 - dist_b1_a1;
 548                     dist_b1_a2 = dist_b1_b2 - dist_b1_a2;
 549                     std::swap(dist_a1_b1, dist_a1_b2);
 550                 }
 551
 552                 segment_ratio<calc_t> ra_from(dist_b1_a1, dist_b1_b2);
 553                 segment_ratio<calc_t> ra_to(dist_b1_a2, dist_b1_b2);
 554                 segment_ratio<calc_t> rb_from(dist_a1_b1, dist_a1_a2);
 555                 segment_ratio<calc_t> rb_to(dist_a1_b2, dist_a1_a2);
 556
 557                 return Policy::segments_collinear(a, b, opposite,
 558                     a1_on_b, a2_on_b, b1_on_a, b2_on_a,
 559                     ra_from, ra_to, rb_from, rb_to);
 560             }
 561         }
 562         else // crossing or touching
 563         {
 564             if (a_is_point || b_is_point)
 565             {
 566                 return Policy::disjoint();
 567             }
 568
 569             calc_t lon = 0, lat = 0;
 570             intersection_point_flag ip_flag;
 571             calc_t dist_a1_a2, dist_a1_i1, dist_b1_b2, dist_b1_i1;
 572             if (calculate_ip_data(a1, a2, b1, b2,
 573                                   a1_lon, a1_lat, a2_lon, a2_lat,
 574                                   b1_lon, b1_lat, b2_lon, b2_lat,
 575                                   res_a1_a2, res_a1_b1, res_a1_b2,
 576                                   res_b1_b2, res_b1_a1, res_b1_a2,
 577                                   sides, spheroid,
 578                                   lon, lat,
 579                                   dist_a1_a2, dist_a1_i1, dist_b1_b2, dist_b1_i1,
 580                                   ip_flag))
 581             {
 582                 // NOTE: If segment was reversed sides and segment ratios has to be altered
 583                 if (is_a_reversed)
 584                 {
 585                     // sides
 586                     sides_reverse_segment<0>(sides);
 587                     // distance for ratio
 588                     dist_a1_i1 = dist_a1_a2 - dist_a1_i1;
 589                     // ip flag
 590                     ip_flag_reverse_segment(ip_flag, ipi_at_a1, ipi_at_a2);
 591                 }
 592                 if (is_b_reversed)
 593                 {
 594                     // sides
 595                     sides_reverse_segment<1>(sides);
 596                     // distance for ratio
 597                     dist_b1_i1 = dist_b1_b2 - dist_b1_i1;
 598                     // ip flag
 599                     ip_flag_reverse_segment(ip_flag, ipi_at_b1, ipi_at_b2);
 600                 }
 601
 602                 // intersects
 603                 segment_intersection_info
 604                     <
 605                         calc_t,
 606                         segment_ratio<calc_t>
 607                     > sinfo;
 608
 609                 sinfo.lon = lon;
 610                 sinfo.lat = lat;
 611                 sinfo.robust_ra.assign(dist_a1_i1, dist_a1_a2);
 612                 sinfo.robust_rb.assign(dist_b1_i1, dist_b1_b2);
 613                 sinfo.ip_flag = ip_flag;
 614
 615                 return Policy::segments_crosses(sides, sinfo, a, b);
 616             }
 617             else
 618             {
 619                 return Policy::disjoint();
 620             }
 621         }
 622     }
 623
 624     template <typename Policy, typename CalcT, typename Segment, typename Point1, typename Point2, typename ResultInverse>
 625     static inline typename Policy::return_type
 626         collinear_one_degenerated(Segment const& segment, bool degenerated_a,
 627                                   Point1 const& a1, Point1 const& a2,
 628                                   Point2 const& b1, Point2 const& b2,
 629                                   ResultInverse const& res_a1_a2,
 630                                   ResultInverse const& res_a1_b1,
 631                                   ResultInverse const& res_a1_b2,
 632                                   bool is_other_reversed,
 633                                   bool degen_neq_coords)
 634     {
 635         CalcT dist_1_2, dist_1_o;
 636         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))
 637         {
 638             return Policy::disjoint();
 639         }
 640
 641         // NOTE: If segment was reversed segment ratio has to be altered
 642         if (is_other_reversed)
 643         {
 644             // distance for ratio
 645             dist_1_o = dist_1_2 - dist_1_o;
 646         }
 647
 648         return Policy::one_degenerate(segment, segment_ratio<CalcT>(dist_1_o, dist_1_2), degenerated_a);
 649     }
 650
 651     // TODO: instead of checks below test bi against a1 and a2 here?
 652     //       in order to make this independent from is_near()
 653     template <typename Point1, typename Point2, typename ResultInverse, typename CalcT>
 654     static inline bool calculate_collinear_data(Point1 const& a1, Point1 const& a2, // in
 655                                                 Point2 const& b1, Point2 const& /*b2*/, // in
 656                                                 ResultInverse const& res_a1_a2,     // in
 657                                                 ResultInverse const& res_a1_b1,     // in
 658                                                 ResultInverse const& res_a1_b2,     // in
 659                                                 CalcT& dist_a1_a2,                  // out
 660                                                 CalcT& dist_a1_b1,                  // out
 661                                                 bool degen_neq_coords = false)      // in
 662     {
 663         dist_a1_a2 = res_a1_a2.distance;
 664
 665         dist_a1_b1 = res_a1_b1.distance;
 666         if (! same_direction(res_a1_b1.azimuth, res_a1_a2.azimuth))
 667         {
 668             dist_a1_b1 = -dist_a1_b1;
 669         }
 670
 671         // if b1 is close a1
 672         if (is_endpoint_equal(dist_a1_b1, a1, b1))
 673         {
 674             dist_a1_b1 = 0;
 675             return true;
 676         }
 677         // if b1 is close a2
 678         else if (is_endpoint_equal(dist_a1_a2 - dist_a1_b1, a2, b1))
 679         {
 680             dist_a1_b1 = dist_a1_a2;
 681             return true;
 682         }
 683
 684         // check the other endpoint of degenerated segment near a pole
 685         if (degen_neq_coords)
 686         {
 687             static CalcT const c0 = 0;
 688             if (math::equals(res_a1_b2.distance, c0))
 689             {
 690                 dist_a1_b1 = 0;
 691                 return true;
 692             }
 693             else if (math::equals(dist_a1_a2 - res_a1_b2.distance, c0))
 694             {
 695                 dist_a1_b1 = dist_a1_a2;
 696                 return true;
 697             }
 698         }
 699
 700         // or i1 is on b
 701         return segment_ratio<CalcT>(dist_a1_b1, dist_a1_a2).on_segment();
 702     }
 703
 704     template <typename Point1, typename Point2, typename CalcT, typename ResultInverse, typename Spheroid_>
 705     static inline bool calculate_ip_data(Point1 const& a1, Point1 const& a2,       // in
 706                                          Point2 const& b1, Point2 const& b2,       // in
 707                                          CalcT const& a1_lon, CalcT const& a1_lat, // in
 708                                          CalcT const& a2_lon, CalcT const& a2_lat, // in
 709                                          CalcT const& b1_lon, CalcT const& b1_lat, // in
 710                                          CalcT const& b2_lon, CalcT const& b2_lat, // in
 711                                          ResultInverse const& res_a1_a2,           // in
 712                                          ResultInverse const& res_a1_b1,           // in
 713                                          ResultInverse const& res_a1_b2,           // in
 714                                          ResultInverse const& res_b1_b2,           // in
 715                                          ResultInverse const& res_b1_a1,           // in
 716                                          ResultInverse const& res_b1_a2,           // in
 717                                          side_info const& sides,                   // in
 718                                          Spheroid_ const& spheroid,                // in
 719                                          CalcT & lon, CalcT & lat,             // out
 720                                          CalcT& dist_a1_a2, CalcT& dist_a1_ip, // out
 721                                          CalcT& dist_b1_b2, CalcT& dist_b1_ip, // out
 722                                          intersection_point_flag& ip_flag)     // out
 723     {
 724         dist_a1_a2 = res_a1_a2.distance;
 725         dist_b1_b2 = res_b1_b2.distance;
 726
 727         // assign the IP if some endpoints overlap
 728         if (equals_point_point(a1, b1))
 729         {
 730             lon = a1_lon;
 731             lat = a1_lat;
 732             dist_a1_ip = 0;
 733             dist_b1_ip = 0;
 734             ip_flag = ipi_at_a1;
 735             return true;
 736         }
 737         else if (equals_point_point(a1, b2))
 738         {
 739             lon = a1_lon;
 740             lat = a1_lat;
 741             dist_a1_ip = 0;
 742             dist_b1_ip = dist_b1_b2;
 743             ip_flag = ipi_at_a1;
 744             return true;
 745         }
 746         else if (equals_point_point(a2, b1))
 747         {
 748             lon = a2_lon;
 749             lat = a2_lat;
 750             dist_a1_ip = dist_a1_a2;
 751             dist_b1_ip = 0;
 752             ip_flag = ipi_at_a2;
 753             return true;
 754         }
 755         else if (equals_point_point(a2, b2))
 756         {
 757             lon = a2_lon;
 758             lat = a2_lat;
 759             dist_a1_ip = dist_a1_a2;
 760             dist_b1_ip = dist_b1_b2;
 761             ip_flag = ipi_at_a2;
 762             return true;
 763         }
 764
 765         // at this point we know that the endpoints doesn't overlap
 766         // check cases when an endpoint lies on the other geodesic
 767         if (sides.template get<0, 0>() == 0) // a1 wrt b
 768         {
 769             if (res_b1_a1.distance <= res_b1_b2.distance
 770                 && same_direction(res_b1_a1.azimuth, res_b1_b2.azimuth))
 771             {
 772                 lon = a1_lon;
 773                 lat = a1_lat;
 774                 dist_a1_ip = 0;
 775                 dist_b1_ip = res_b1_a1.distance;
 776                 ip_flag = ipi_at_a1;
 777                 return true;
 778             }
 779             else
 780             {
 781                 return false;
 782             }
 783         }
 784         else if (sides.template get<0, 1>() == 0) // a2 wrt b
 785         {
 786             if (res_b1_a2.distance <= res_b1_b2.distance
 787                 && same_direction(res_b1_a2.azimuth, res_b1_b2.azimuth))
 788             {
 789                 lon = a2_lon;
 790                 lat = a2_lat;
 791                 dist_a1_ip = res_a1_a2.distance;
 792                 dist_b1_ip = res_b1_a2.distance;
 793                 ip_flag = ipi_at_a2;
 794                 return true;
 795             }
 796             else
 797             {
 798                 return false;
 799             }
 800         }
 801         else if (sides.template get<1, 0>() == 0) // b1 wrt a
 802         {
 803             if (res_a1_b1.distance <= res_a1_a2.distance
 804                 && same_direction(res_a1_b1.azimuth, res_a1_a2.azimuth))
 805             {
 806                 lon = b1_lon;
 807                 lat = b1_lat;
 808                 dist_a1_ip = res_a1_b1.distance;
 809                 dist_b1_ip = 0;
 810                 ip_flag = ipi_at_b1;
 811                 return true;
 812             }
 813             else
 814             {
 815                 return false;
 816             }
 817         }
 818         else if (sides.template get<1, 1>() == 0) // b2 wrt a
 819         {
 820             if (res_a1_b2.distance <= res_a1_a2.distance
 821                 && same_direction(res_a1_b2.azimuth, res_a1_a2.azimuth))
 822             {
 823                 lon = b2_lon;
 824                 lat = b2_lat;
 825                 dist_a1_ip = res_a1_b2.distance;
 826                 dist_b1_ip = res_b1_b2.distance;
 827                 ip_flag = ipi_at_b2;
 828                 return true;
 829             }
 830             else
 831             {
 832                 return false;
 833             }
 834         }
 835
 836         // At this point neither the endpoints overlaps
 837         // nor any andpoint lies on the other geodesic
 838         // So the endpoints should lie on the opposite sides of both geodesics
 839
 840         bool const ok = formula::sjoberg_intersection<CalcT, FormulaPolicy::template inverse, Order>
 841                         ::apply(a1_lon, a1_lat, a2_lon, a2_lat, res_a1_a2.azimuth,
 842                                 b1_lon, b1_lat, b2_lon, b2_lat, res_b1_b2.azimuth,
 843                                 lon, lat, spheroid);
 844
 845         if (! ok)
 846         {
 847             return false;
 848         }
 849
 850         typedef typename FormulaPolicy::template inverse<CalcT, true, true, false, false, false> inverse_dist_azi;
 851         typedef typename inverse_dist_azi::result_type inverse_result;
 852
 853         inverse_result const res_a1_ip = inverse_dist_azi::apply(a1_lon, a1_lat, lon, lat, spheroid);
 854         dist_a1_ip = res_a1_ip.distance;
 855         if (! same_direction(res_a1_ip.azimuth, res_a1_a2.azimuth))
 856         {
 857             dist_a1_ip = -dist_a1_ip;
 858         }
 859
 860         bool is_on_a = segment_ratio<CalcT>(dist_a1_ip, dist_a1_a2).on_segment();
 861         // NOTE: not fully consistent with equals_point_point() since radians are always used.
 862         bool is_on_a1 = math::equals(lon, a1_lon) && math::equals(lat, a1_lat);
 863         bool is_on_a2 = math::equals(lon, a2_lon) && math::equals(lat, a2_lat);
 864
 865         if (! (is_on_a || is_on_a1 || is_on_a2))
 866         {
 867             return false;
 868         }
 869
 870         inverse_result const res_b1_ip = inverse_dist_azi::apply(b1_lon, b1_lat, lon, lat, spheroid);
 871         dist_b1_ip = res_b1_ip.distance;
 872         if (! same_direction(res_b1_ip.azimuth, res_b1_b2.azimuth))
 873         {
 874             dist_b1_ip = -dist_b1_ip;
 875         }
 876
 877         bool is_on_b = segment_ratio<CalcT>(dist_b1_ip, dist_b1_b2).on_segment();
 878         // NOTE: not fully consistent with equals_point_point() since radians are always used.
 879         bool is_on_b1 = math::equals(lon, b1_lon) && math::equals(lat, b1_lat);
 880         bool is_on_b2 = math::equals(lon, b2_lon) && math::equals(lat, b2_lat);
 881
 882         if (! (is_on_b || is_on_b1 || is_on_b2))
 883         {
 884             return false;
 885         }
 886
 887         typedef typename FormulaPolicy::template inverse<CalcT, true, false, false, false, false> inverse_dist;
 888
 889         ip_flag = ipi_inters;
 890
 891         if (is_on_b1)
 892         {
 893             lon = b1_lon;
 894             lat = b1_lat;
 895             dist_a1_ip = inverse_dist::apply(a1_lon, a1_lat, lon, lat, spheroid).distance; // for consistency
 896             dist_b1_ip = 0;
 897             ip_flag = ipi_at_b1;
 898         }
 899         else if (is_on_b2)
 900         {
 901             lon = b2_lon;
 902             lat = b2_lat;
 903             dist_a1_ip = inverse_dist::apply(a1_lon, a1_lat, lon, lat, spheroid).distance; // for consistency
 904             dist_b1_ip = res_b1_b2.distance;
 905             ip_flag = ipi_at_b2;
 906         }
 907
 908         if (is_on_a1)
 909         {
 910             lon = a1_lon;
 911             lat = a1_lat;
 912             dist_a1_ip = 0;
 913             dist_b1_ip = inverse_dist::apply(b1_lon, b1_lat, lon, lat, spheroid).distance; // for consistency
 914             ip_flag = ipi_at_a1;
 915         }
 916         else if (is_on_a2)
 917         {
 918             lon = a2_lon;
 919             lat = a2_lat;
 920             dist_a1_ip = res_a1_a2.distance;
 921             dist_b1_ip = inverse_dist::apply(b1_lon, b1_lat, lon, lat, spheroid).distance; // for consistency
 922             ip_flag = ipi_at_a2;
 923         }
 924
 925         return true;
 926     }
 927
 928     template <typename CalcT, typename P1, typename P2>
 929     static inline bool is_endpoint_equal(CalcT const& dist,
 930                                          P1 const& ai, P2 const& b1)
 931     {
 932         static CalcT const c0 = 0;
 933         return is_near(dist) && (math::equals(dist, c0) || equals_point_point(ai, b1));
 934     }
 935
 936     template <typename CalcT>
 937     static inline bool is_near(CalcT const& dist)
 938     {
 939         // NOTE: This strongly depends on the Inverse method
 940         CalcT const small_number = CalcT(std::is_same<CalcT, float>::value ? 0.0001 : 0.00000001);
 941         return math::abs(dist) <= small_number;
 942     }
 943
 944     template <typename ProjCoord1, typename ProjCoord2>
 945     static inline int position_value(ProjCoord1 const& ca1,
 946                                      ProjCoord2 const& cb1,
 947                                      ProjCoord2 const& cb2)
 948     {
 949         // S1x  0   1    2     3   4
 950         // S2       |---------->
 951         return math::equals(ca1, cb1) ? 1
 952              : math::equals(ca1, cb2) ? 3
 953              : cb1 < cb2 ?
 954                 ( ca1 < cb1 ? 0
 955                 : ca1 > cb2 ? 4
 956                 : 2 )
 957               : ( ca1 > cb1 ? 0
 958                 : ca1 < cb2 ? 4
 959                 : 2 );
 960     }
 961
 962     template <typename CalcT>
 963     static inline bool same_direction(CalcT const& azimuth1, CalcT const& azimuth2)
 964     {
 965         // distance between two angles normalized to (-180, 180]
 966         CalcT const angle_diff = math::longitude_distance_signed<radian>(azimuth1, azimuth2);
 967         return math::abs(angle_diff) <= math::half_pi<CalcT>();
 968     }
 969
 970     template <int Which>
 971     static inline void sides_reverse_segment(side_info & sides)
 972     {
 973         // names assuming segment A is reversed (Which == 0)
 974         int a1_wrt_b = sides.template get<Which, 0>();
 975         int a2_wrt_b = sides.template get<Which, 1>();
 976         std::swap(a1_wrt_b, a2_wrt_b);
 977         sides.template set<Which>(a1_wrt_b, a2_wrt_b);
 978         int b1_wrt_a = sides.template get<1 - Which, 0>();
 979         int b2_wrt_a = sides.template get<1 - Which, 1>();
 980         sides.template set<1 - Which>(-b1_wrt_a, -b2_wrt_a);
 981     }
 982
 983     static inline void ip_flag_reverse_segment(intersection_point_flag & ip_flag,
 984                                                intersection_point_flag const& ipi_at_p1,
 985                                                intersection_point_flag const& ipi_at_p2)
 986     {
 987         ip_flag = ip_flag == ipi_at_p1 ? ipi_at_p2 :
 988                   ip_flag == ipi_at_p2 ? ipi_at_p1 :
 989                   ip_flag;
 990     }
 991
 992     template <typename Point1, typename Point2>
 993     static inline bool equals_point_point(Point1 const& point1, Point2 const& point2)
 994     {
 995         return detail::equals::equals_point_point(point1, point2,
 996                                                   point_in_point_strategy_type());
 997     }
 998
 999 private:
1000     Spheroid m_spheroid;
1001 };
1002
1003
1004 }} // namespace strategy::intersection
1005
1006 }} // namespace boost::geometry
1007
1008
1009 #endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_INTERSECTION_HPP