boost/geometry/strategies/geographic/intersection.hpp

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