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

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