strategies/geographic/intersection.hpp

0001 // Boost.Geometry
0002
0003 // Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland.
0004
0005 // Copyright (c) 2016-2021, Oracle and/or its affiliates.
0006 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0007
0008 // Use, modification and distribution is subject to the Boost Software License,
0009 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0010 // http://www.boost.org/LICENSE_1_0.txt)
0011
0012 #ifndef BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_INTERSECTION_HPP
0013 #define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_INTERSECTION_HPP
0014
0015 #include <algorithm>
0016 #include <type_traits>
0017
0018 #include <boost/geometry/core/cs.hpp>
0019 #include <boost/geometry/core/access.hpp>
0020 #include <boost/geometry/core/radian_access.hpp>
0021 #include <boost/geometry/core/tags.hpp>
0022
0023 #include <boost/geometry/algorithms/detail/assign_values.hpp>
0024 #include <boost/geometry/algorithms/detail/assign_indexed_point.hpp>
0025 #include <boost/geometry/algorithms/detail/equals/point_point.hpp>
0026 #include <boost/geometry/algorithms/detail/recalculate.hpp>
0027
0028 #include <boost/geometry/formulas/andoyer_inverse.hpp>
0029 #include <boost/geometry/formulas/sjoberg_intersection.hpp>
0030 #include <boost/geometry/formulas/spherical.hpp>
0031 #include <boost/geometry/formulas/unit_spheroid.hpp>
0032
0033 #include <boost/geometry/geometries/concepts/point_concept.hpp>
0034 #include <boost/geometry/geometries/concepts/segment_concept.hpp>
0035 #include <boost/geometry/geometries/segment.hpp>
0036
0037 #include <boost/geometry/policies/robustness/segment_ratio.hpp>
0038
0039 #include <boost/geometry/srs/spheroid.hpp>
0040
0041 #include <boost/geometry/strategy/geographic/area.hpp>
0042 #include <boost/geometry/strategy/geographic/envelope.hpp>
0043 #include <boost/geometry/strategy/geographic/expand_segment.hpp>
0044 #include <boost/geometry/strategy/spherical/expand_box.hpp>
0045
0046 #include <boost/geometry/strategies/geographic/disjoint_segment_box.hpp>
0047 #include <boost/geometry/strategies/geographic/distance.hpp>
0048 #include <boost/geometry/strategies/geographic/parameters.hpp>
0049 #include <boost/geometry/strategies/geographic/point_in_poly_winding.hpp>
0050 #include <boost/geometry/strategies/geographic/side.hpp>
0051 #include <boost/geometry/strategies/spherical/disjoint_box_box.hpp>
0052 #include <boost/geometry/strategies/spherical/point_in_point.hpp>
0053 #include <boost/geometry/strategies/intersection.hpp>
0054 #include <boost/geometry/strategies/intersection_result.hpp>
0055 #include <boost/geometry/strategies/side_info.hpp>
0056
0057 #include <boost/geometry/util/math.hpp>
0058 #include <boost/geometry/util/select_calculation_type.hpp>
0059
0060
0061 namespace boost { namespace geometry
0062 {
0063
0064 namespace strategy { namespace intersection
0065 {
0066
0067 // CONSIDER: Improvement of the robustness/accuracy/repeatability by
0068 // moving all segments to 0 longitude
0069 // picking latitudes closer to 0
0070 // etc.
0071
0072 template
0073 <
0074     typename FormulaPolicy = strategy::andoyer,
0075     std::size_t Order = strategy::default_order<FormulaPolicy>::value,
0076     typename Spheroid = srs::spheroid<double>,
0077     typename CalculationType = void
0078 >
0079 struct geographic_segments
0080 {
0081     typedef geographic_tag cs_tag;
0082
0083     enum intersection_point_flag { ipi_inters = 0, ipi_at_a1, ipi_at_a2, ipi_at_b1, ipi_at_b2 };
0084
0085     template <typename CoordinateType, typename SegmentRatio>
0086     struct segment_intersection_info
0087     {
0088         template <typename Point, typename Segment1, typename Segment2>
0089         void calculate(Point& point, Segment1 const& a, Segment2 const& b) const
0090         {
0091             if (ip_flag == ipi_inters)
0092             {
0093                 // TODO: assign the rest of coordinates
0094                 set_from_radian<0>(point, lon);
0095                 set_from_radian<1>(point, lat);
0096             }
0097             else if (ip_flag == ipi_at_a1)
0098             {
0099                 detail::assign_point_from_index<0>(a, point);
0100             }
0101             else if (ip_flag == ipi_at_a2)
0102             {
0103                 detail::assign_point_from_index<1>(a, point);
0104             }
0105             else if (ip_flag == ipi_at_b1)
0106             {
0107                 detail::assign_point_from_index<0>(b, point);
0108             }
0109             else // ip_flag == ipi_at_b2
0110             {
0111                 detail::assign_point_from_index<1>(b, point);
0112             }
0113         }
0114
0115         CoordinateType lon;
0116         CoordinateType lat;
0117         SegmentRatio robust_ra;
0118         SegmentRatio robust_rb;
0119         intersection_point_flag ip_flag;
0120     };
0121
0122     explicit geographic_segments(Spheroid const& spheroid = Spheroid())
0123         : m_spheroid(spheroid)
0124     {}
0125
0126     Spheroid model() const
0127     {
0128         return m_spheroid;
0129     }
0130
0131     // Relate segments a and b
0132     template
0133     <
0134         typename UniqueSubRange1,
0135         typename UniqueSubRange2,
0136         typename Policy
0137     >
0138     inline typename Policy::return_type apply(UniqueSubRange1 const& range_p,
0139                                               UniqueSubRange2 const& range_q,
0140                                               Policy const&) const
0141     {
0142         typedef typename UniqueSubRange1::point_type point1_type;
0143         typedef typename UniqueSubRange2::point_type point2_type;
0144         typedef model::referring_segment<point1_type const> segment_type1;
0145         typedef model::referring_segment<point2_type const> segment_type2;
0146
0147         BOOST_CONCEPT_ASSERT( (concepts::ConstPoint<point1_type>) );
0148         BOOST_CONCEPT_ASSERT( (concepts::ConstPoint<point2_type>) );
0149
0150         /*
0151         typename coordinate_type<Point1>::type
0152             const a1_lon = get<0>(a1),
0153             const a2_lon = get<0>(a2);
0154         typename coordinate_type<Point2>::type
0155             const b1_lon = get<0>(b1),
0156             const b2_lon = get<0>(b2);
0157         bool is_a_reversed = a1_lon > a2_lon || a1_lon == a2_lon && get<1>(a1) > get<1>(a2);
0158         bool is_b_reversed = b1_lon > b2_lon || b1_lon == b2_lon && get<1>(b1) > get<1>(b2);
0159         */
0160
0161         point1_type const& p0 = range_p.at(0);
0162         point1_type const& p1 = range_p.at(1);
0163         point2_type const& q0 = range_q.at(0);
0164         point2_type const& q1 = range_q.at(1);
0165
0166         bool const is_p_reversed = get<1>(p0) > get<1>(p1);
0167         bool const is_q_reversed = get<1>(q0) > get<1>(q1);
0168
0169         // Call apply with original segments and ordered points
0170         return apply<Policy>(segment_type1(p0, p1),
0171                              segment_type2(q0, q1),
0172                              (is_p_reversed ? p1 : p0),
0173                              (is_p_reversed ? p0 : p1),
0174                              (is_q_reversed ? q1 : q0),
0175                              (is_q_reversed ? q0 : q1),
0176                              is_p_reversed, is_q_reversed);
0177     }
0178
0179 private:
0180     // Relate segments a and b
0181     template
0182     <
0183         typename Policy,
0184         typename Segment1,
0185         typename Segment2,
0186         typename Point1,
0187         typename Point2
0188     >
0189     inline typename Policy::return_type apply(Segment1 const& a, Segment2 const& b,
0190                                               Point1 const& a1, Point1 const& a2,
0191                                               Point2 const& b1, Point2 const& b2,
0192                                               bool is_a_reversed, bool is_b_reversed) const
0193     {
0194         BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment1>) );
0195         BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment2>) );
0196
0197         typedef typename select_calculation_type
0198             <Segment1, Segment2, CalculationType>::type calc_t;
0199
0200         typedef srs::spheroid<calc_t> spheroid_type;
0201
0202         static const calc_t c0 = 0;
0203
0204         // normalized spheroid
0205         spheroid_type spheroid = formula::unit_spheroid<spheroid_type>(m_spheroid);
0206
0207         // TODO: check only 2 first coordinates here?
0208         bool a_is_point = equals_point_point(a1, a2);
0209         bool b_is_point = equals_point_point(b1, b2);
0210
0211         if(a_is_point && b_is_point)
0212         {
0213             return equals_point_point(a1, b2)
0214                 ? Policy::degenerate(a, true)
0215                 : Policy::disjoint()
0216                 ;
0217         }
0218
0219         calc_t const a1_lon = get_as_radian<0>(a1);
0220         calc_t const a1_lat = get_as_radian<1>(a1);
0221         calc_t const a2_lon = get_as_radian<0>(a2);
0222         calc_t const a2_lat = get_as_radian<1>(a2);
0223         calc_t const b1_lon = get_as_radian<0>(b1);
0224         calc_t const b1_lat = get_as_radian<1>(b1);
0225         calc_t const b2_lon = get_as_radian<0>(b2);
0226         calc_t const b2_lat = get_as_radian<1>(b2);
0227
0228         side_info sides;
0229
0230         // NOTE: potential optimization, don't calculate distance at this point
0231         // this would require to reimplement inverse strategy to allow
0232         // calculation of distance if needed, probably also storing intermediate
0233         // results somehow inside an object.
0234         typedef typename FormulaPolicy::template inverse<calc_t, true, true, false, false, false> inverse_dist_azi;
0235         typedef typename inverse_dist_azi::result_type inverse_result;
0236
0237         // TODO: no need to call inverse formula if we know that the points are equal
0238         // distance can be set to 0 in this case and azimuth may be not calculated
0239         bool is_equal_a1_b1 = equals_point_point(a1, b1);
0240         bool is_equal_a2_b1 = equals_point_point(a2, b1);
0241         bool degen_neq_coords = false;
0242
0243         inverse_result res_b1_b2, res_b1_a1, res_b1_a2;
0244         if (! b_is_point)
0245         {
0246             res_b1_b2 = inverse_dist_azi::apply(b1_lon, b1_lat, b2_lon, b2_lat, spheroid);
0247             if (math::equals(res_b1_b2.distance, c0))
0248             {
0249                 b_is_point = true;
0250                 degen_neq_coords = true;
0251             }
0252             else
0253             {
0254                 res_b1_a1 = inverse_dist_azi::apply(b1_lon, b1_lat, a1_lon, a1_lat, spheroid);
0255                 if (math::equals(res_b1_a1.distance, c0))
0256                 {
0257                     is_equal_a1_b1 = true;
0258                 }
0259                 res_b1_a2 = inverse_dist_azi::apply(b1_lon, b1_lat, a2_lon, a2_lat, spheroid);
0260                 if (math::equals(res_b1_a2.distance, c0))
0261                 {
0262                     is_equal_a2_b1 = true;
0263                 }
0264                 sides.set<0>(is_equal_a1_b1 ? 0 : formula::azimuth_side_value(res_b1_a1.azimuth, res_b1_b2.azimuth),
0265                              is_equal_a2_b1 ? 0 : formula::azimuth_side_value(res_b1_a2.azimuth, res_b1_b2.azimuth));
0266                 if (sides.same<0>())
0267                 {
0268                     // Both points are at the same side of other segment, we can leave
0269                     return Policy::disjoint();
0270                 }
0271             }
0272         }
0273
0274         bool is_equal_a1_b2 = equals_point_point(a1, b2);
0275
0276         inverse_result res_a1_a2, res_a1_b1, res_a1_b2;
0277         if (! a_is_point)
0278         {
0279             res_a1_a2 = inverse_dist_azi::apply(a1_lon, a1_lat, a2_lon, a2_lat, spheroid);
0280             if (math::equals(res_a1_a2.distance, c0))
0281             {
0282                 a_is_point = true;
0283                 degen_neq_coords = true;
0284             }
0285             else
0286             {
0287                 res_a1_b1 = inverse_dist_azi::apply(a1_lon, a1_lat, b1_lon, b1_lat, spheroid);
0288                 if (math::equals(res_a1_b1.distance, c0))
0289                 {
0290                     is_equal_a1_b1 = true;
0291                 }
0292                 res_a1_b2 = inverse_dist_azi::apply(a1_lon, a1_lat, b2_lon, b2_lat, spheroid);
0293                 if (math::equals(res_a1_b2.distance, c0))
0294                 {
0295                     is_equal_a1_b2 = true;
0296                 }
0297                 sides.set<1>(is_equal_a1_b1 ? 0 : formula::azimuth_side_value(res_a1_b1.azimuth, res_a1_a2.azimuth),
0298                              is_equal_a1_b2 ? 0 : formula::azimuth_side_value(res_a1_b2.azimuth, res_a1_a2.azimuth));
0299                 if (sides.same<1>())
0300                 {
0301                     // Both points are at the same side of other segment, we can leave
0302                     return Policy::disjoint();
0303                 }
0304             }
0305         }
0306
0307         if(a_is_point && b_is_point)
0308         {
0309             return is_equal_a1_b2
0310                 ? Policy::degenerate(a, true)
0311                 : Policy::disjoint()
0312                 ;
0313         }
0314
0315         // NOTE: at this point the segments may still be disjoint
0316         // NOTE: at this point one of the segments may be degenerated
0317
0318         bool collinear = sides.collinear();
0319
0320         if (! collinear)
0321         {
0322             // WARNING: the side strategy doesn't have the info about the other
0323             // segment so it may return results inconsistent with this intersection
0324             // strategy, as it checks both segments for consistency
0325
0326             if (sides.get<0, 0>() == 0 && sides.get<0, 1>() == 0)
0327             {
0328                 collinear = true;
0329                 sides.set<1>(0, 0);
0330             }
0331             else if (sides.get<1, 0>() == 0 && sides.get<1, 1>() == 0)
0332             {
0333                 collinear = true;
0334                 sides.set<0>(0, 0);
0335             }
0336         }
0337
0338         if (collinear)
0339         {
0340             if (a_is_point)
0341             {
0342                 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);
0343             }
0344             else if (b_is_point)
0345             {
0346                 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);
0347             }
0348             else
0349             {
0350                 calc_t dist_a1_a2, dist_a1_b1, dist_a1_b2;
0351                 calc_t dist_b1_b2, dist_b1_a1, dist_b1_a2;
0352                 // use shorter segment
0353                 if (res_a1_a2.distance <= res_b1_b2.distance)
0354                 {
0355                     calculate_collinear_data(a1, a2, b1, b2, res_a1_a2, res_a1_b1, res_a1_b2, dist_a1_a2, dist_a1_b1);
0356                     calculate_collinear_data(a1, a2, b2, b1, res_a1_a2, res_a1_b2, res_a1_b1, dist_a1_a2, dist_a1_b2);
0357                     dist_b1_b2 = dist_a1_b2 - dist_a1_b1;
0358                     dist_b1_a1 = -dist_a1_b1;
0359                     dist_b1_a2 = dist_a1_a2 - dist_a1_b1;
0360                 }
0361                 else
0362                 {
0363                     calculate_collinear_data(b1, b2, a1, a2, res_b1_b2, res_b1_a1, res_b1_a2, dist_b1_b2, dist_b1_a1);
0364                     calculate_collinear_data(b1, b2, a2, a1, res_b1_b2, res_b1_a2, res_b1_a1, dist_b1_b2, dist_b1_a2);
0365                     dist_a1_a2 = dist_b1_a2 - dist_b1_a1;
0366                     dist_a1_b1 = -dist_b1_a1;
0367                     dist_a1_b2 = dist_b1_b2 - dist_b1_a1;
0368                 }
0369
0370                 // NOTE: this is probably not needed
0371                 int a1_on_b = position_value(c0, dist_a1_b1, dist_a1_b2);
0372                 int a2_on_b = position_value(dist_a1_a2, dist_a1_b1, dist_a1_b2);
0373                 int b1_on_a = position_value(c0, dist_b1_a1, dist_b1_a2);
0374                 int b2_on_a = position_value(dist_b1_b2, dist_b1_a1, dist_b1_a2);
0375
0376                 if ((a1_on_b < 1 && a2_on_b < 1) || (a1_on_b > 3 && a2_on_b > 3))
0377                 {
0378                     return Policy::disjoint();
0379                 }
0380
0381                 if (a1_on_b == 1)
0382                 {
0383                     dist_b1_a1 = 0;
0384                     dist_a1_b1 = 0;
0385                 }
0386                 else if (a1_on_b == 3)
0387                 {
0388                     dist_b1_a1 = dist_b1_b2;
0389                     dist_a1_b2 = 0;
0390                 }
0391
0392                 if (a2_on_b == 1)
0393                 {
0394                     dist_b1_a2 = 0;
0395                     dist_a1_b1 = dist_a1_a2;
0396                 }
0397                 else if (a2_on_b == 3)
0398                 {
0399                     dist_b1_a2 = dist_b1_b2;
0400                     dist_a1_b2 = dist_a1_a2;
0401                 }
0402
0403                 bool opposite = ! same_direction(res_a1_a2.azimuth, res_b1_b2.azimuth);
0404
0405                 // NOTE: If segment was reversed opposite, positions and segment ratios has to be altered
0406                 if (is_a_reversed)
0407                 {
0408                     // opposite
0409                     opposite = ! opposite;
0410                     // positions
0411                     std::swap(a1_on_b, a2_on_b);
0412                     b1_on_a = 4 - b1_on_a;
0413                     b2_on_a = 4 - b2_on_a;
0414                     // distances for ratios
0415                     std::swap(dist_b1_a1, dist_b1_a2);
0416                     dist_a1_b1 = dist_a1_a2 - dist_a1_b1;
0417                     dist_a1_b2 = dist_a1_a2 - dist_a1_b2;
0418                 }
0419                 if (is_b_reversed)
0420                 {
0421                     // opposite
0422                     opposite = ! opposite;
0423                     // positions
0424                     a1_on_b = 4 - a1_on_b;
0425                     a2_on_b = 4 - a2_on_b;
0426                     std::swap(b1_on_a, b2_on_a);
0427                     // distances for ratios
0428                     dist_b1_a1 = dist_b1_b2 - dist_b1_a1;
0429                     dist_b1_a2 = dist_b1_b2 - dist_b1_a2;
0430                     std::swap(dist_a1_b1, dist_a1_b2);
0431                 }
0432
0433                 segment_ratio<calc_t> ra_from(dist_b1_a1, dist_b1_b2);
0434                 segment_ratio<calc_t> ra_to(dist_b1_a2, dist_b1_b2);
0435                 segment_ratio<calc_t> rb_from(dist_a1_b1, dist_a1_a2);
0436                 segment_ratio<calc_t> rb_to(dist_a1_b2, dist_a1_a2);
0437
0438                 return Policy::segments_collinear(a, b, opposite,
0439                     a1_on_b, a2_on_b, b1_on_a, b2_on_a,
0440                     ra_from, ra_to, rb_from, rb_to);
0441             }
0442         }
0443         else // crossing or touching
0444         {
0445             if (a_is_point || b_is_point)
0446             {
0447                 return Policy::disjoint();
0448             }
0449
0450             calc_t lon = 0, lat = 0;
0451             intersection_point_flag ip_flag;
0452             calc_t dist_a1_a2, dist_a1_i1, dist_b1_b2, dist_b1_i1;
0453             if (calculate_ip_data(a1, a2, b1, b2,
0454                                   a1_lon, a1_lat, a2_lon, a2_lat,
0455                                   b1_lon, b1_lat, b2_lon, b2_lat,
0456                                   res_a1_a2, res_a1_b1, res_a1_b2,
0457                                   res_b1_b2, res_b1_a1, res_b1_a2,
0458                                   sides, spheroid,
0459                                   lon, lat,
0460                                   dist_a1_a2, dist_a1_i1, dist_b1_b2, dist_b1_i1,
0461                                   ip_flag))
0462             {
0463                 // NOTE: If segment was reversed sides and segment ratios has to be altered
0464                 if (is_a_reversed)
0465                 {
0466                     // sides
0467                     sides_reverse_segment<0>(sides);
0468                     // distance for ratio
0469                     dist_a1_i1 = dist_a1_a2 - dist_a1_i1;
0470                     // ip flag
0471                     ip_flag_reverse_segment(ip_flag, ipi_at_a1, ipi_at_a2);
0472                 }
0473                 if (is_b_reversed)
0474                 {
0475                     // sides
0476                     sides_reverse_segment<1>(sides);
0477                     // distance for ratio
0478                     dist_b1_i1 = dist_b1_b2 - dist_b1_i1;
0479                     // ip flag
0480                     ip_flag_reverse_segment(ip_flag, ipi_at_b1, ipi_at_b2);
0481                 }
0482
0483                 // intersects
0484                 segment_intersection_info
0485                     <
0486                         calc_t,
0487                         segment_ratio<calc_t>
0488                     > sinfo;
0489
0490                 sinfo.lon = lon;
0491                 sinfo.lat = lat;
0492                 sinfo.robust_ra.assign(dist_a1_i1, dist_a1_a2);
0493                 sinfo.robust_rb.assign(dist_b1_i1, dist_b1_b2);
0494                 sinfo.ip_flag = ip_flag;
0495
0496                 return Policy::segments_crosses(sides, sinfo, a, b);
0497             }
0498             else
0499             {
0500                 return Policy::disjoint();
0501             }
0502         }
0503     }
0504
0505     template <typename Policy, typename CalcT, typename Segment, typename Point1, typename Point2, typename ResultInverse>
0506     static inline typename Policy::return_type
0507         collinear_one_degenerated(Segment const& segment, bool degenerated_a,
0508                                   Point1 const& a1, Point1 const& a2,
0509                                   Point2 const& b1, Point2 const& b2,
0510                                   ResultInverse const& res_a1_a2,
0511                                   ResultInverse const& res_a1_b1,
0512                                   ResultInverse const& res_a1_b2,
0513                                   bool is_other_reversed,
0514                                   bool degen_neq_coords)
0515     {
0516         CalcT dist_1_2, dist_1_o;
0517         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))
0518         {
0519             return Policy::disjoint();
0520         }
0521
0522         // NOTE: If segment was reversed segment ratio has to be altered
0523         if (is_other_reversed)
0524         {
0525             // distance for ratio
0526             dist_1_o = dist_1_2 - dist_1_o;
0527         }
0528
0529         return Policy::one_degenerate(segment, segment_ratio<CalcT>(dist_1_o, dist_1_2), degenerated_a);
0530     }
0531
0532     // TODO: instead of checks below test bi against a1 and a2 here?
0533     //       in order to make this independent from is_near()
0534     template <typename Point1, typename Point2, typename ResultInverse, typename CalcT>
0535     static inline bool calculate_collinear_data(Point1 const& a1, Point1 const& a2, // in
0536                                                 Point2 const& b1, Point2 const& /*b2*/, // in
0537                                                 ResultInverse const& res_a1_a2,     // in
0538                                                 ResultInverse const& res_a1_b1,     // in
0539                                                 ResultInverse const& res_a1_b2,     // in
0540                                                 CalcT& dist_a1_a2,                  // out
0541                                                 CalcT& dist_a1_b1,                  // out
0542                                                 bool degen_neq_coords = false)      // in
0543     {
0544         dist_a1_a2 = res_a1_a2.distance;
0545
0546         dist_a1_b1 = res_a1_b1.distance;
0547         if (! same_direction(res_a1_b1.azimuth, res_a1_a2.azimuth))
0548         {
0549             dist_a1_b1 = -dist_a1_b1;
0550         }
0551
0552         // if b1 is close a1
0553         if (is_endpoint_equal(dist_a1_b1, a1, b1))
0554         {
0555             dist_a1_b1 = 0;
0556             return true;
0557         }
0558         // if b1 is close a2
0559         else if (is_endpoint_equal(dist_a1_a2 - dist_a1_b1, a2, b1))
0560         {
0561             dist_a1_b1 = dist_a1_a2;
0562             return true;
0563         }
0564
0565         // check the other endpoint of degenerated segment near a pole
0566         if (degen_neq_coords)
0567         {
0568             static CalcT const c0 = 0;
0569             if (math::equals(res_a1_b2.distance, c0))
0570             {
0571                 dist_a1_b1 = 0;
0572                 return true;
0573             }
0574             else if (math::equals(dist_a1_a2 - res_a1_b2.distance, c0))
0575             {
0576                 dist_a1_b1 = dist_a1_a2;
0577                 return true;
0578             }
0579         }
0580
0581         // or i1 is on b
0582         return segment_ratio<CalcT>(dist_a1_b1, dist_a1_a2).on_segment();
0583     }
0584
0585     template <typename Point1, typename Point2, typename CalcT, typename ResultInverse, typename Spheroid_>
0586     static inline bool calculate_ip_data(Point1 const& a1, Point1 const& a2,       // in
0587                                          Point2 const& b1, Point2 const& b2,       // in
0588                                          CalcT const& a1_lon, CalcT const& a1_lat, // in
0589                                          CalcT const& a2_lon, CalcT const& a2_lat, // in
0590                                          CalcT const& b1_lon, CalcT const& b1_lat, // in
0591                                          CalcT const& b2_lon, CalcT const& b2_lat, // in
0592                                          ResultInverse const& res_a1_a2,           // in
0593                                          ResultInverse const& res_a1_b1,           // in
0594                                          ResultInverse const& res_a1_b2,           // in
0595                                          ResultInverse const& res_b1_b2,           // in
0596                                          ResultInverse const& res_b1_a1,           // in
0597                                          ResultInverse const& res_b1_a2,           // in
0598                                          side_info const& sides,                   // in
0599                                          Spheroid_ const& spheroid,                // in
0600                                          CalcT & lon, CalcT & lat,             // out
0601                                          CalcT& dist_a1_a2, CalcT& dist_a1_ip, // out
0602                                          CalcT& dist_b1_b2, CalcT& dist_b1_ip, // out
0603                                          intersection_point_flag& ip_flag)     // out
0604     {
0605         dist_a1_a2 = res_a1_a2.distance;
0606         dist_b1_b2 = res_b1_b2.distance;
0607
0608         // assign the IP if some endpoints overlap
0609         if (equals_point_point(a1, b1))
0610         {
0611             lon = a1_lon;
0612             lat = a1_lat;
0613             dist_a1_ip = 0;
0614             dist_b1_ip = 0;
0615             ip_flag = ipi_at_a1;
0616             return true;
0617         }
0618         else if (equals_point_point(a1, b2))
0619         {
0620             lon = a1_lon;
0621             lat = a1_lat;
0622             dist_a1_ip = 0;
0623             dist_b1_ip = dist_b1_b2;
0624             ip_flag = ipi_at_a1;
0625             return true;
0626         }
0627         else if (equals_point_point(a2, b1))
0628         {
0629             lon = a2_lon;
0630             lat = a2_lat;
0631             dist_a1_ip = dist_a1_a2;
0632             dist_b1_ip = 0;
0633             ip_flag = ipi_at_a2;
0634             return true;
0635         }
0636         else if (equals_point_point(a2, b2))
0637         {
0638             lon = a2_lon;
0639             lat = a2_lat;
0640             dist_a1_ip = dist_a1_a2;
0641             dist_b1_ip = dist_b1_b2;
0642             ip_flag = ipi_at_a2;
0643             return true;
0644         }
0645
0646         // at this point we know that the endpoints doesn't overlap
0647         // check cases when an endpoint lies on the other geodesic
0648         if (sides.template get<0, 0>() == 0) // a1 wrt b
0649         {
0650             if (res_b1_a1.distance <= res_b1_b2.distance
0651                 && same_direction(res_b1_a1.azimuth, res_b1_b2.azimuth))
0652             {
0653                 lon = a1_lon;
0654                 lat = a1_lat;
0655                 dist_a1_ip = 0;
0656                 dist_b1_ip = res_b1_a1.distance;
0657                 ip_flag = ipi_at_a1;
0658                 return true;
0659             }
0660             else
0661             {
0662                 return false;
0663             }
0664         }
0665         else if (sides.template get<0, 1>() == 0) // a2 wrt b
0666         {
0667             if (res_b1_a2.distance <= res_b1_b2.distance
0668                 && same_direction(res_b1_a2.azimuth, res_b1_b2.azimuth))
0669             {
0670                 lon = a2_lon;
0671                 lat = a2_lat;
0672                 dist_a1_ip = res_a1_a2.distance;
0673                 dist_b1_ip = res_b1_a2.distance;
0674                 ip_flag = ipi_at_a2;
0675                 return true;
0676             }
0677             else
0678             {
0679                 return false;
0680             }
0681         }
0682         else if (sides.template get<1, 0>() == 0) // b1 wrt a
0683         {
0684             if (res_a1_b1.distance <= res_a1_a2.distance
0685                 && same_direction(res_a1_b1.azimuth, res_a1_a2.azimuth))
0686             {
0687                 lon = b1_lon;
0688                 lat = b1_lat;
0689                 dist_a1_ip = res_a1_b1.distance;
0690                 dist_b1_ip = 0;
0691                 ip_flag = ipi_at_b1;
0692                 return true;
0693             }
0694             else
0695             {
0696                 return false;
0697             }
0698         }
0699         else if (sides.template get<1, 1>() == 0) // b2 wrt a
0700         {
0701             if (res_a1_b2.distance <= res_a1_a2.distance
0702                 && same_direction(res_a1_b2.azimuth, res_a1_a2.azimuth))
0703             {
0704                 lon = b2_lon;
0705                 lat = b2_lat;
0706                 dist_a1_ip = res_a1_b2.distance;
0707                 dist_b1_ip = res_b1_b2.distance;
0708                 ip_flag = ipi_at_b2;
0709                 return true;
0710             }
0711             else
0712             {
0713                 return false;
0714             }
0715         }
0716
0717         // At this point neither the endpoints overlaps
0718         // nor any andpoint lies on the other geodesic
0719         // So the endpoints should lie on the opposite sides of both geodesics
0720
0721         bool const ok = formula::sjoberg_intersection<CalcT, FormulaPolicy::template inverse, Order>
0722                         ::apply(a1_lon, a1_lat, a2_lon, a2_lat, res_a1_a2.azimuth,
0723                                 b1_lon, b1_lat, b2_lon, b2_lat, res_b1_b2.azimuth,
0724                                 lon, lat, spheroid);
0725
0726         if (! ok)
0727         {
0728             return false;
0729         }
0730
0731         typedef typename FormulaPolicy::template inverse<CalcT, true, true, false, false, false> inverse_dist_azi;
0732         typedef typename inverse_dist_azi::result_type inverse_result;
0733
0734         inverse_result const res_a1_ip = inverse_dist_azi::apply(a1_lon, a1_lat, lon, lat, spheroid);
0735         dist_a1_ip = res_a1_ip.distance;
0736         if (! same_direction(res_a1_ip.azimuth, res_a1_a2.azimuth))
0737         {
0738             dist_a1_ip = -dist_a1_ip;
0739         }
0740
0741         bool is_on_a = segment_ratio<CalcT>(dist_a1_ip, dist_a1_a2).on_segment();
0742         // NOTE: not fully consistent with equals_point_point() since radians are always used.
0743         bool is_on_a1 = math::equals(lon, a1_lon) && math::equals(lat, a1_lat);
0744         bool is_on_a2 = math::equals(lon, a2_lon) && math::equals(lat, a2_lat);
0745
0746         if (! (is_on_a || is_on_a1 || is_on_a2))
0747         {
0748             return false;
0749         }
0750
0751         inverse_result const res_b1_ip = inverse_dist_azi::apply(b1_lon, b1_lat, lon, lat, spheroid);
0752         dist_b1_ip = res_b1_ip.distance;
0753         if (! same_direction(res_b1_ip.azimuth, res_b1_b2.azimuth))
0754         {
0755             dist_b1_ip = -dist_b1_ip;
0756         }
0757
0758         bool is_on_b = segment_ratio<CalcT>(dist_b1_ip, dist_b1_b2).on_segment();
0759         // NOTE: not fully consistent with equals_point_point() since radians are always used.
0760         bool is_on_b1 = math::equals(lon, b1_lon) && math::equals(lat, b1_lat);
0761         bool is_on_b2 = math::equals(lon, b2_lon) && math::equals(lat, b2_lat);
0762
0763         if (! (is_on_b || is_on_b1 || is_on_b2))
0764         {
0765             return false;
0766         }
0767
0768         typedef typename FormulaPolicy::template inverse<CalcT, true, false, false, false, false> inverse_dist;
0769
0770         ip_flag = ipi_inters;
0771
0772         if (is_on_b1)
0773         {
0774             lon = b1_lon;
0775             lat = b1_lat;
0776             dist_a1_ip = inverse_dist::apply(a1_lon, a1_lat, lon, lat, spheroid).distance; // for consistency
0777             dist_b1_ip = 0;
0778             ip_flag = ipi_at_b1;
0779         }
0780         else if (is_on_b2)
0781         {
0782             lon = b2_lon;
0783             lat = b2_lat;
0784             dist_a1_ip = inverse_dist::apply(a1_lon, a1_lat, lon, lat, spheroid).distance; // for consistency
0785             dist_b1_ip = res_b1_b2.distance;
0786             ip_flag = ipi_at_b2;
0787         }
0788
0789         if (is_on_a1)
0790         {
0791             lon = a1_lon;
0792             lat = a1_lat;
0793             dist_a1_ip = 0;
0794             dist_b1_ip = inverse_dist::apply(b1_lon, b1_lat, lon, lat, spheroid).distance; // for consistency
0795             ip_flag = ipi_at_a1;
0796         }
0797         else if (is_on_a2)
0798         {
0799             lon = a2_lon;
0800             lat = a2_lat;
0801             dist_a1_ip = res_a1_a2.distance;
0802             dist_b1_ip = inverse_dist::apply(b1_lon, b1_lat, lon, lat, spheroid).distance; // for consistency
0803             ip_flag = ipi_at_a2;
0804         }
0805
0806         return true;
0807     }
0808
0809     template <typename CalcT, typename P1, typename P2>
0810     static inline bool is_endpoint_equal(CalcT const& dist,
0811                                          P1 const& ai, P2 const& b1)
0812     {
0813         static CalcT const c0 = 0;
0814         return is_near(dist) && (math::equals(dist, c0) || equals_point_point(ai, b1));
0815     }
0816
0817     template <typename CalcT>
0818     static inline bool is_near(CalcT const& dist)
0819     {
0820         // NOTE: This strongly depends on the Inverse method
0821         CalcT const small_number = CalcT(std::is_same<CalcT, float>::value ? 0.0001 : 0.00000001);
0822         return math::abs(dist) <= small_number;
0823     }
0824
0825     template <typename ProjCoord1, typename ProjCoord2>
0826     static inline int position_value(ProjCoord1 const& ca1,
0827                                      ProjCoord2 const& cb1,
0828                                      ProjCoord2 const& cb2)
0829     {
0830         // S1x  0   1    2     3   4
0831         // S2       |---------->
0832         return math::equals(ca1, cb1) ? 1
0833              : math::equals(ca1, cb2) ? 3
0834              : cb1 < cb2 ?
0835                 ( ca1 < cb1 ? 0
0836                 : ca1 > cb2 ? 4
0837                 : 2 )
0838               : ( ca1 > cb1 ? 0
0839                 : ca1 < cb2 ? 4
0840                 : 2 );
0841     }
0842
0843     template <typename CalcT>
0844     static inline bool same_direction(CalcT const& azimuth1, CalcT const& azimuth2)
0845     {
0846         // distance between two angles normalized to (-180, 180]
0847         CalcT const angle_diff = math::longitude_distance_signed<radian>(azimuth1, azimuth2);
0848         return math::abs(angle_diff) <= math::half_pi<CalcT>();
0849     }
0850
0851     template <int Which>
0852     static inline void sides_reverse_segment(side_info & sides)
0853     {
0854         // names assuming segment A is reversed (Which == 0)
0855         int a1_wrt_b = sides.template get<Which, 0>();
0856         int a2_wrt_b = sides.template get<Which, 1>();
0857         std::swap(a1_wrt_b, a2_wrt_b);
0858         sides.template set<Which>(a1_wrt_b, a2_wrt_b);
0859         int b1_wrt_a = sides.template get<1 - Which, 0>();
0860         int b2_wrt_a = sides.template get<1 - Which, 1>();
0861         sides.template set<1 - Which>(-b1_wrt_a, -b2_wrt_a);
0862     }
0863
0864     static inline void ip_flag_reverse_segment(intersection_point_flag & ip_flag,
0865                                                intersection_point_flag const& ipi_at_p1,
0866                                                intersection_point_flag const& ipi_at_p2)
0867     {
0868         ip_flag = ip_flag == ipi_at_p1 ? ipi_at_p2 :
0869                   ip_flag == ipi_at_p2 ? ipi_at_p1 :
0870                   ip_flag;
0871     }
0872
0873     template <typename Point1, typename Point2>
0874     static inline bool equals_point_point(Point1 const& point1, Point2 const& point2)
0875     {
0876         return strategy::within::spherical_point_point::apply(point1, point2);
0877     }
0878
0879 private:
0880     Spheroid m_spheroid;
0881 };
0882
0883
0884 }} // namespace strategy::intersection
0885
0886 }} // namespace boost::geometry
0887
0888
0889 #endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_INTERSECTION_HPP