strategies/spherical/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_SPHERICAL_INTERSECTION_HPP
0013 #define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_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/arithmetic/arithmetic.hpp>
0029 #include <boost/geometry/arithmetic/cross_product.hpp>
0030 #include <boost/geometry/arithmetic/dot_product.hpp>
0031 #include <boost/geometry/arithmetic/normalize.hpp>
0032 #include <boost/geometry/formulas/spherical.hpp>
0033
0034 #include <boost/geometry/geometries/concepts/point_concept.hpp>
0035 #include <boost/geometry/geometries/concepts/segment_concept.hpp>
0036 #include <boost/geometry/geometries/segment.hpp>
0037
0038 #include <boost/geometry/policies/robustness/segment_ratio.hpp>
0039
0040 #include <boost/geometry/strategy/spherical/area.hpp>
0041 #include <boost/geometry/strategy/spherical/envelope.hpp>
0042 #include <boost/geometry/strategy/spherical/expand_box.hpp>
0043 #include <boost/geometry/strategy/spherical/expand_segment.hpp>
0044
0045 #include <boost/geometry/strategies/covered_by.hpp>
0046 #include <boost/geometry/strategies/intersection.hpp>
0047 #include <boost/geometry/strategies/intersection_result.hpp>
0048 #include <boost/geometry/strategies/side.hpp>
0049 #include <boost/geometry/strategies/side_info.hpp>
0050 #include <boost/geometry/strategies/spherical/disjoint_box_box.hpp>
0051 #include <boost/geometry/strategies/spherical/disjoint_segment_box.hpp>
0052 #include <boost/geometry/strategies/spherical/distance_haversine.hpp>
0053 #include <boost/geometry/strategies/spherical/point_in_point.hpp>
0054 #include <boost/geometry/strategies/spherical/point_in_poly_winding.hpp>
0055 #include <boost/geometry/strategies/spherical/ssf.hpp>
0056 #include <boost/geometry/strategies/within.hpp>
0057
0058 #include <boost/geometry/util/math.hpp>
0059 #include <boost/geometry/util/select_calculation_type.hpp>
0060
0061
0062 namespace boost { namespace geometry
0063 {
0064
0065 namespace strategy { namespace intersection
0066 {
0067
0068 // NOTE:
0069 // The coordinates of crossing IP may be calculated with small precision in some cases.
0070 // For double, near the equator noticed error ~1e-9 so far greater than
0071 // machine epsilon which is ~1e-16. This error is ~0.04m.
0072 // E.g. consider two cases, one near the origin and the second one rotated by 90 deg around Z or SN axis.
0073 // After the conversion from spherical degrees to cartesian 3d the following coordinates
0074 // are calculated:
0075 // for sph (-1 -1,  1 1) deg cart3d ys are -0.017449748351250485 and  0.017449748351250485
0076 // for sph (89 -1, 91 1) deg cart3d xs are  0.017449748351250571 and -0.017449748351250450
0077 // During the conversion degrees must first be converted to radians and then radians
0078 // are passed into trigonometric functions. The error may have several causes:
0079 // 1. Radians cannot represent exactly the same angles as degrees.
0080 // 2. Different longitudes are passed into sin() for x, corresponding to cos() for y,
0081 //    and for different angle the error of the result may be different.
0082 // 3. These non-corresponding cartesian coordinates are used in calculation,
0083 //    e.g. multiplied several times in cross and dot products.
0084 // If it was a problem this strategy could e.g. "normalize" longitudes before the conversion using the source units
0085 // by rotating the globe around Z axis, so moving longitudes always the same way towards the origin,
0086 // assuming this could help which is not clear.
0087 // For now, intersection points near the endpoints are checked explicitly if needed (if the IP is near the endpoint)
0088 // to generate precise result for them. Only the crossing (i) case may suffer from lower precision.
0089
0090 template
0091 <
0092     typename CalcPolicy,
0093     typename CalculationType = void
0094 >
0095 struct ecef_segments
0096 {
0097     typedef spherical_tag cs_tag;
0098
0099     enum intersection_point_flag { ipi_inters = 0, ipi_at_a1, ipi_at_a2, ipi_at_b1, ipi_at_b2 };
0100
0101     // segment_intersection_info cannot outlive relate_ecef_segments
0102     template <typename CoordinateType, typename SegmentRatio, typename Vector3d>
0103     struct segment_intersection_info
0104     {
0105         segment_intersection_info(CalcPolicy const& calc)
0106             : calc_policy(calc)
0107         {}
0108
0109         template <typename Point, typename Segment1, typename Segment2>
0110         void calculate(Point& point, Segment1 const& a, Segment2 const& b) const
0111         {
0112             if (ip_flag == ipi_inters)
0113             {
0114                 // TODO: assign the rest of coordinates
0115                 point = calc_policy.template from_cart3d<Point>(intersection_point);
0116             }
0117             else if (ip_flag == ipi_at_a1)
0118             {
0119                 detail::assign_point_from_index<0>(a, point);
0120             }
0121             else if (ip_flag == ipi_at_a2)
0122             {
0123                 detail::assign_point_from_index<1>(a, point);
0124             }
0125             else if (ip_flag == ipi_at_b1)
0126             {
0127                 detail::assign_point_from_index<0>(b, point);
0128             }
0129             else // ip_flag == ipi_at_b2
0130             {
0131                 detail::assign_point_from_index<1>(b, point);
0132             }
0133         }
0134
0135         Vector3d intersection_point;
0136         SegmentRatio robust_ra;
0137         SegmentRatio robust_rb;
0138         intersection_point_flag ip_flag;
0139
0140         CalcPolicy const& calc_policy;
0141     };
0142
0143     // Relate segments a and b
0144     template
0145     <
0146         typename UniqueSubRange1,
0147         typename UniqueSubRange2,
0148         typename Policy
0149     >
0150     static inline typename Policy::return_type
0151         apply(UniqueSubRange1 const& range_p, UniqueSubRange2 const& range_q,
0152               Policy const&)
0153     {
0154         // For now create it using default constructor. In the future it could
0155         //  be stored in strategy. However then apply() wouldn't be static and
0156         //  all relops and setops would have to take the strategy or model.
0157         // Initialize explicitly to prevent compiler errors in case of PoD type
0158         CalcPolicy const calc_policy = CalcPolicy();
0159
0160         typedef typename UniqueSubRange1::point_type point1_type;
0161         typedef typename UniqueSubRange2::point_type point2_type;
0162
0163         BOOST_CONCEPT_ASSERT( (concepts::ConstPoint<point1_type>) );
0164         BOOST_CONCEPT_ASSERT( (concepts::ConstPoint<point2_type>) );
0165
0166         point1_type const& a1 = range_p.at(0);
0167         point1_type const& a2 = range_p.at(1);
0168         point2_type const& b1 = range_q.at(0);
0169         point2_type const& b2 = range_q.at(1);
0170
0171         typedef model::referring_segment<point1_type const> segment1_type;
0172         typedef model::referring_segment<point2_type const> segment2_type;
0173         segment1_type const a(a1, a2);
0174         segment2_type const b(b1, b2);
0175
0176         // TODO: check only 2 first coordinates here?
0177         bool a_is_point = equals_point_point(a1, a2);
0178         bool b_is_point = equals_point_point(b1, b2);
0179
0180         if(a_is_point && b_is_point)
0181         {
0182             return equals_point_point(a1, b2)
0183                 ? Policy::degenerate(a, true)
0184                 : Policy::disjoint()
0185                 ;
0186         }
0187
0188         typedef typename select_calculation_type
0189             <segment1_type, segment2_type, CalculationType>::type calc_t;
0190
0191         calc_t const c0 = 0;
0192         calc_t const c1 = 1;
0193
0194         typedef model::point<calc_t, 3, cs::cartesian> vec3d_t;
0195
0196         vec3d_t const a1v = calc_policy.template to_cart3d<vec3d_t>(a1);
0197         vec3d_t const a2v = calc_policy.template to_cart3d<vec3d_t>(a2);
0198         vec3d_t const b1v = calc_policy.template to_cart3d<vec3d_t>(b1);
0199         vec3d_t const b2v = calc_policy.template to_cart3d<vec3d_t>(b2);
0200
0201         bool degen_neq_coords = false;
0202         side_info sides;
0203
0204         typename CalcPolicy::template plane<vec3d_t>
0205             plane2 = calc_policy.get_plane(b1v, b2v);
0206
0207         calc_t dist_b1_b2 = 0;
0208         if (! b_is_point)
0209         {
0210             calculate_dist(b1v, b2v, plane2, dist_b1_b2);
0211             if (math::equals(dist_b1_b2, c0))
0212             {
0213                 degen_neq_coords = true;
0214                 b_is_point = true;
0215                 dist_b1_b2 = 0;
0216             }
0217             else
0218             {
0219                 // not normalized normals, the same as in side strategy
0220                 sides.set<0>(plane2.side_value(a1v), plane2.side_value(a2v));
0221                 if (sides.same<0>())
0222                 {
0223                     // Both points are at same side of other segment, we can leave
0224                     return Policy::disjoint();
0225                 }
0226             }
0227         }
0228
0229         typename CalcPolicy::template plane<vec3d_t>
0230             plane1 = calc_policy.get_plane(a1v, a2v);
0231
0232         calc_t dist_a1_a2 = 0;
0233         if (! a_is_point)
0234         {
0235             calculate_dist(a1v, a2v, plane1, dist_a1_a2);
0236             if (math::equals(dist_a1_a2, c0))
0237             {
0238                 degen_neq_coords = true;
0239                 a_is_point = true;
0240                 dist_a1_a2 = 0;
0241             }
0242             else
0243             {
0244                 // not normalized normals, the same as in side strategy
0245                 sides.set<1>(plane1.side_value(b1v), plane1.side_value(b2v));
0246                 if (sides.same<1>())
0247                 {
0248                     // Both points are at same side of other segment, we can leave
0249                     return Policy::disjoint();
0250                 }
0251             }
0252         }
0253
0254         // NOTE: at this point the segments may still be disjoint
0255
0256         calc_t len1 = 0;
0257         // point or opposite sides of a sphere/spheroid, assume point
0258         if (! a_is_point && ! detail::vec_normalize(plane1.normal, len1))
0259         {
0260             a_is_point = true;
0261             if (sides.get<0, 0>() == 0 || sides.get<0, 1>() == 0)
0262             {
0263                 sides.set<0>(0, 0);
0264             }
0265         }
0266
0267         calc_t len2 = 0;
0268         if (! b_is_point && ! detail::vec_normalize(plane2.normal, len2))
0269         {
0270             b_is_point = true;
0271             if (sides.get<1, 0>() == 0 || sides.get<1, 1>() == 0)
0272             {
0273                 sides.set<1>(0, 0);
0274             }
0275         }
0276
0277         // check both degenerated once more
0278         if (a_is_point && b_is_point)
0279         {
0280             return equals_point_point(a1, b2)
0281                 ? Policy::degenerate(a, true)
0282                 : Policy::disjoint()
0283                 ;
0284         }
0285
0286         // NOTE: at this point the segments may still be disjoint
0287         // NOTE: at this point one of the segments may be degenerated
0288
0289         bool collinear = sides.collinear();
0290
0291         if (! collinear)
0292         {
0293             // NOTE: for some approximations it's possible that both points may lie
0294             // on the same geodesic but still some of the sides may be != 0.
0295             // This is e.g. true for long segments represented as elliptic arcs
0296             // with origin different than the center of the coordinate system.
0297             // So make the sides consistent
0298
0299             // WARNING: the side strategy doesn't have the info about the other
0300             // segment so it may return results inconsistent with this intersection
0301             // strategy, as it checks both segments for consistency
0302
0303             if (sides.get<0, 0>() == 0 && sides.get<0, 1>() == 0)
0304             {
0305                 collinear = true;
0306                 sides.set<1>(0, 0);
0307             }
0308             else if (sides.get<1, 0>() == 0 && sides.get<1, 1>() == 0)
0309             {
0310                 collinear = true;
0311                 sides.set<0>(0, 0);
0312             }
0313         }
0314
0315         calc_t dot_n1n2 = dot_product(plane1.normal, plane2.normal);
0316
0317         // NOTE: this is technically not needed since theoretically above sides
0318         //       are calculated, but just in case check the normals.
0319         //       Have in mind that SSF side strategy doesn't check this.
0320         // collinear if normals are equal or opposite: cos(a) in {-1, 1}
0321         if (! collinear && math::equals(math::abs(dot_n1n2), c1))
0322         {
0323             collinear = true;
0324             sides.set<0>(0, 0);
0325             sides.set<1>(0, 0);
0326         }
0327
0328         if (collinear)
0329         {
0330             if (a_is_point)
0331             {
0332                 return collinear_one_degenerated<Policy, calc_t>(a, true, b1, b2, a1, a2, b1v, b2v,
0333                                                                  plane2, a1v, a2v, dist_b1_b2, degen_neq_coords);
0334             }
0335             else if (b_is_point)
0336             {
0337                 // b2 used to be consistent with (degenerated) checks above (is it needed?)
0338                 return collinear_one_degenerated<Policy, calc_t>(b, false, a1, a2, b1, b2, a1v, a2v,
0339                                                                  plane1, b1v, b2v, dist_a1_a2, degen_neq_coords);
0340             }
0341             else
0342             {
0343                 calc_t dist_a1_b1, dist_a1_b2;
0344                 calc_t dist_b1_a1, dist_b1_a2;
0345                 calculate_collinear_data(a1, a2, b1, b2, a1v, a2v, plane1, b1v, b2v, dist_a1_a2, dist_a1_b1);
0346                 calculate_collinear_data(a1, a2, b2, b1, a1v, a2v, plane1, b2v, b1v, dist_a1_a2, dist_a1_b2);
0347                 calculate_collinear_data(b1, b2, a1, a2, b1v, b2v, plane2, a1v, a2v, dist_b1_b2, dist_b1_a1);
0348                 calculate_collinear_data(b1, b2, a2, a1, b1v, b2v, plane2, a2v, a1v, dist_b1_b2, dist_b1_a2);
0349                 // NOTE: The following optimization causes problems with consitency
0350                 // It may either be caused by numerical issues or the way how distance is coded:
0351                 //   as cosine of angle scaled and translated, see: calculate_dist()
0352                 /*dist_b1_b2 = dist_a1_b2 - dist_a1_b1;
0353                 dist_b1_a1 = -dist_a1_b1;
0354                 dist_b1_a2 = dist_a1_a2 - dist_a1_b1;
0355                 dist_a1_a2 = dist_b1_a2 - dist_b1_a1;
0356                 dist_a1_b1 = -dist_b1_a1;
0357                 dist_a1_b2 = dist_b1_b2 - dist_b1_a1;*/
0358
0359                 segment_ratio<calc_t> ra_from(dist_b1_a1, dist_b1_b2);
0360                 segment_ratio<calc_t> ra_to(dist_b1_a2, dist_b1_b2);
0361                 segment_ratio<calc_t> rb_from(dist_a1_b1, dist_a1_a2);
0362                 segment_ratio<calc_t> rb_to(dist_a1_b2, dist_a1_a2);
0363
0364                 // NOTE: this is probably not needed
0365                 int const a1_wrt_b = position_value(c0, dist_a1_b1, dist_a1_b2);
0366                 int const a2_wrt_b = position_value(dist_a1_a2, dist_a1_b1, dist_a1_b2);
0367                 int const b1_wrt_a = position_value(c0, dist_b1_a1, dist_b1_a2);
0368                 int const b2_wrt_a = position_value(dist_b1_b2, dist_b1_a1, dist_b1_a2);
0369
0370                 if (a1_wrt_b == 1)
0371                 {
0372                     ra_from.assign(0, dist_b1_b2);
0373                     rb_from.assign(0, dist_a1_a2);
0374                 }
0375                 else if (a1_wrt_b == 3)
0376                 {
0377                     ra_from.assign(dist_b1_b2, dist_b1_b2);
0378                     rb_to.assign(0, dist_a1_a2);
0379                 }
0380
0381                 if (a2_wrt_b == 1)
0382                 {
0383                     ra_to.assign(0, dist_b1_b2);
0384                     rb_from.assign(dist_a1_a2, dist_a1_a2);
0385                 }
0386                 else if (a2_wrt_b == 3)
0387                 {
0388                     ra_to.assign(dist_b1_b2, dist_b1_b2);
0389                     rb_to.assign(dist_a1_a2, dist_a1_a2);
0390                 }
0391
0392                 if ((a1_wrt_b < 1 && a2_wrt_b < 1) || (a1_wrt_b > 3 && a2_wrt_b > 3))
0393                 {
0394                     return Policy::disjoint();
0395                 }
0396
0397                 bool const opposite = dot_n1n2 < c0;
0398
0399                 return Policy::segments_collinear(a, b, opposite,
0400                     a1_wrt_b, a2_wrt_b, b1_wrt_a, b2_wrt_a,
0401                     ra_from, ra_to, rb_from, rb_to);
0402             }
0403         }
0404         else // crossing
0405         {
0406             if (a_is_point || b_is_point)
0407             {
0408                 return Policy::disjoint();
0409             }
0410
0411             vec3d_t i1;
0412             intersection_point_flag ip_flag;
0413             calc_t dist_a1_i1, dist_b1_i1;
0414             if (calculate_ip_data(a1, a2, b1, b2, a1v, a2v, b1v, b2v,
0415                                   plane1, plane2, calc_policy,
0416                                   sides, dist_a1_a2, dist_b1_b2,
0417                                   i1, dist_a1_i1, dist_b1_i1, ip_flag))
0418             {
0419                 // intersects
0420                 segment_intersection_info
0421                     <
0422                         calc_t,
0423                         segment_ratio<calc_t>,
0424                         vec3d_t
0425                     > sinfo(calc_policy);
0426
0427                 sinfo.robust_ra.assign(dist_a1_i1, dist_a1_a2);
0428                 sinfo.robust_rb.assign(dist_b1_i1, dist_b1_b2);
0429                 sinfo.intersection_point = i1;
0430                 sinfo.ip_flag = ip_flag;
0431
0432                 return Policy::segments_crosses(sides, sinfo, a, b);
0433             }
0434             else
0435             {
0436                 return Policy::disjoint();
0437             }
0438         }
0439     }
0440
0441 private:
0442     template <typename Policy, typename CalcT, typename Segment, typename Point1, typename Point2, typename Vec3d, typename Plane>
0443     static inline typename Policy::return_type
0444         collinear_one_degenerated(Segment const& segment, bool degenerated_a,
0445                                   Point1 const& a1, Point1 const& a2,
0446                                   Point2 const& b1, Point2 const& b2,
0447                                   Vec3d const& a1v, Vec3d const& a2v,
0448                                   Plane const& plane,
0449                                   Vec3d const& b1v, Vec3d const& b2v,
0450                                   CalcT const& dist_1_2,
0451                                   bool degen_neq_coords)
0452     {
0453         CalcT dist_1_o;
0454         return ! calculate_collinear_data(a1, a2, b1, b2, a1v, a2v, plane, b1v, b2v, dist_1_2, dist_1_o, degen_neq_coords)
0455                 ? Policy::disjoint()
0456                 : Policy::one_degenerate(segment, segment_ratio<CalcT>(dist_1_o, dist_1_2), degenerated_a);
0457     }
0458
0459     template <typename Point1, typename Point2, typename Vec3d, typename Plane, typename CalcT>
0460     static inline bool calculate_collinear_data(Point1 const& a1, Point1 const& a2, // in
0461                                                 Point2 const& b1, Point2 const& /*b2*/, // in
0462                                                 Vec3d const& a1v,                   // in
0463                                                 Vec3d const& a2v,                   // in
0464                                                 Plane const& plane1,                // in
0465                                                 Vec3d const& b1v,                   // in
0466                                                 Vec3d const& b2v,                   // in
0467                                                 CalcT const& dist_a1_a2,            // in
0468                                                 CalcT& dist_a1_b1,                  // out
0469                                                 bool degen_neq_coords = false)      // in
0470     {
0471         // calculate dist_a1_b1
0472         calculate_dist(a1v, a2v, plane1, b1v, dist_a1_b1);
0473
0474         // if b1 is equal to a1
0475         if (is_endpoint_equal(dist_a1_b1, a1, b1))
0476         {
0477             dist_a1_b1 = 0;
0478             return true;
0479         }
0480         // or b1 is equal to a2
0481         else if (is_endpoint_equal(dist_a1_a2 - dist_a1_b1, a2, b1))
0482         {
0483             dist_a1_b1 = dist_a1_a2;
0484             return true;
0485         }
0486
0487         // check the other endpoint of degenerated segment near a pole
0488         if (degen_neq_coords)
0489         {
0490             static CalcT const c0 = 0;
0491
0492             CalcT dist_a1_b2 = 0;
0493             calculate_dist(a1v, a2v, plane1, b2v, dist_a1_b2);
0494
0495             if (math::equals(dist_a1_b2, c0))
0496             {
0497                 dist_a1_b1 = 0;
0498                 return true;
0499             }
0500             else if (math::equals(dist_a1_a2 - dist_a1_b2, c0))
0501             {
0502                 dist_a1_b1 = dist_a1_a2;
0503                 return true;
0504             }
0505         }
0506
0507         // or i1 is on b
0508         return segment_ratio<CalcT>(dist_a1_b1, dist_a1_a2).on_segment();
0509     }
0510
0511     template <typename Point1, typename Point2, typename Vec3d, typename Plane, typename CalcT>
0512     static inline bool calculate_ip_data(Point1 const& a1, Point1 const& a2, // in
0513                                          Point2 const& b1, Point2 const& b2, // in
0514                                          Vec3d const& a1v, Vec3d const& a2v, // in
0515                                          Vec3d const& b1v, Vec3d const& b2v, // in
0516                                          Plane const& plane1,                // in
0517                                          Plane const& plane2,                // in
0518                                          CalcPolicy const& calc_policy,      // in
0519                                          side_info const& sides,             // in
0520                                          CalcT const& dist_a1_a2,            // in
0521                                          CalcT const& dist_b1_b2,            // in
0522                                          Vec3d & ip,                         // out
0523                                          CalcT& dist_a1_ip,                  // out
0524                                          CalcT& dist_b1_ip,                  // out
0525                                          intersection_point_flag& ip_flag)   // out
0526     {
0527         Vec3d ip1, ip2;
0528         calc_policy.intersection_points(plane1, plane2, ip1, ip2);
0529
0530         calculate_dist(a1v, a2v, plane1, ip1, dist_a1_ip);
0531         ip = ip1;
0532
0533         // choose the opposite side of the globe if the distance is shorter
0534         {
0535             CalcT const d = abs_distance(dist_a1_a2, dist_a1_ip);
0536             if (d > CalcT(0))
0537             {
0538                 // TODO: this should be ok not only for sphere
0539                 //       but requires more investigation
0540                 CalcT const dist_a1_i2 = dist_of_i2(dist_a1_ip);
0541                 CalcT const d2 = abs_distance(dist_a1_a2, dist_a1_i2);
0542                 if (d2 < d)
0543                 {
0544                     dist_a1_ip = dist_a1_i2;
0545                     ip = ip2;
0546                 }
0547             }
0548         }
0549
0550         bool is_on_a = false, is_near_a1 = false, is_near_a2 = false;
0551         if (! is_potentially_crossing(dist_a1_a2, dist_a1_ip, is_on_a, is_near_a1, is_near_a2))
0552         {
0553             return false;
0554         }
0555
0556         calculate_dist(b1v, b2v, plane2, ip, dist_b1_ip);
0557
0558         bool is_on_b = false, is_near_b1 = false, is_near_b2 = false;
0559         if (! is_potentially_crossing(dist_b1_b2, dist_b1_ip, is_on_b, is_near_b1, is_near_b2))
0560         {
0561             return false;
0562         }
0563
0564         // reassign the IP if some endpoints overlap
0565         if (is_near_a1)
0566         {
0567             if (is_near_b1 && equals_point_point(a1, b1))
0568             {
0569                 dist_a1_ip = 0;
0570                 dist_b1_ip = 0;
0571                 //i1 = a1v;
0572                 ip_flag = ipi_at_a1;
0573                 return true;
0574             }
0575
0576             if (is_near_b2 && equals_point_point(a1, b2))
0577             {
0578                 dist_a1_ip = 0;
0579                 dist_b1_ip = dist_b1_b2;
0580                 //i1 = a1v;
0581                 ip_flag = ipi_at_a1;
0582                 return true;
0583             }
0584         }
0585
0586         if (is_near_a2)
0587         {
0588             if (is_near_b1 && equals_point_point(a2, b1))
0589             {
0590                 dist_a1_ip = dist_a1_a2;
0591                 dist_b1_ip = 0;
0592                 //i1 = a2v;
0593                 ip_flag = ipi_at_a2;
0594                 return true;
0595             }
0596
0597             if (is_near_b2 && equals_point_point(a2, b2))
0598             {
0599                 dist_a1_ip = dist_a1_a2;
0600                 dist_b1_ip = dist_b1_b2;
0601                 //i1 = a2v;
0602                 ip_flag = ipi_at_a2;
0603                 return true;
0604             }
0605         }
0606
0607         // at this point we know that the endpoints doesn't overlap
0608         // reassign IP and distance if the IP is on a segment and one of
0609         //   the endpoints of the other segment lies on the former segment
0610         if (is_on_a)
0611         {
0612             if (is_near_b1 && sides.template get<1, 0>() == 0) // b1 wrt a
0613             {
0614                 calculate_dist(a1v, a2v, plane1, b1v, dist_a1_ip); // for consistency
0615                 dist_b1_ip = 0;
0616                 //i1 = b1v;
0617                 ip_flag = ipi_at_b1;
0618                 return true;
0619             }
0620
0621             if (is_near_b2 && sides.template get<1, 1>() == 0) // b2 wrt a
0622             {
0623                 calculate_dist(a1v, a2v, plane1, b2v, dist_a1_ip); // for consistency
0624                 dist_b1_ip = dist_b1_b2;
0625                 //i1 = b2v;
0626                 ip_flag = ipi_at_b2;
0627                 return true;
0628             }
0629         }
0630
0631         if (is_on_b)
0632         {
0633             if (is_near_a1 && sides.template get<0, 0>() == 0) // a1 wrt b
0634             {
0635                 dist_a1_ip = 0;
0636                 calculate_dist(b1v, b2v, plane2, a1v, dist_b1_ip); // for consistency
0637                 //i1 = a1v;
0638                 ip_flag = ipi_at_a1;
0639                 return true;
0640             }
0641
0642             if (is_near_a2 && sides.template get<0, 1>() == 0) // a2 wrt b
0643             {
0644                 dist_a1_ip = dist_a1_a2;
0645                 calculate_dist(b1v, b2v, plane2, a2v, dist_b1_ip); // for consistency
0646                 //i1 = a2v;
0647                 ip_flag = ipi_at_a2;
0648                 return true;
0649             }
0650         }
0651
0652         ip_flag = ipi_inters;
0653
0654         return is_on_a && is_on_b;
0655     }
0656
0657     template <typename Vec3d, typename Plane, typename CalcT>
0658     static inline void calculate_dist(Vec3d const& a1v,    // in
0659                                       Vec3d const& a2v,    // in
0660                                       Plane const& plane1, // in
0661                                       CalcT& dist_a1_a2)   // out
0662     {
0663         static CalcT const c1 = 1;
0664         CalcT const cos_a1_a2 = plane1.cos_angle_between(a1v, a2v);
0665         dist_a1_a2 = -cos_a1_a2 + c1; // [1, -1] -> [0, 2] representing [0, pi]
0666     }
0667
0668     template <typename Vec3d, typename Plane, typename CalcT>
0669     static inline void calculate_dist(Vec3d const& a1v,     // in
0670                                       Vec3d const& /*a2v*/, // in
0671                                       Plane const& plane1,  // in
0672                                       Vec3d const& i1,      // in
0673                                       CalcT& dist_a1_i1)    // out
0674     {
0675         static CalcT const c1 = 1;
0676         static CalcT const c2 = 2;
0677         static CalcT const c4 = 4;
0678
0679         bool is_forward = true;
0680         CalcT cos_a1_i1 = plane1.cos_angle_between(a1v, i1, is_forward);
0681         dist_a1_i1 = -cos_a1_i1 + c1; // [0, 2] representing [0, pi]
0682         if (! is_forward) // left or right of a1 on a
0683         {
0684             dist_a1_i1 = -dist_a1_i1; // [0, 2] -> [0, -2] representing [0, -pi]
0685         }
0686         if (dist_a1_i1 <= -c2) // <= -pi
0687         {
0688             dist_a1_i1 += c4; // += 2pi
0689         }
0690     }
0691     /*
0692     template <typename Vec3d, typename Plane, typename CalcT>
0693     static inline void calculate_dists(Vec3d const& a1v,    // in
0694                                        Vec3d const& a2v,    // in
0695                                        Plane const& plane1, // in
0696                                        Vec3d const& i1,     // in
0697                                        CalcT& dist_a1_a2, // out
0698                                        CalcT& dist_a1_i1) // out
0699     {
0700         calculate_dist(a1v, a2v, plane1, dist_a1_a2);
0701         calculate_dist(a1v, a2v, plane1, i1, dist_a1_i1);
0702     }
0703     */
0704     // the dist of the ip on the other side of the sphere
0705     template <typename CalcT>
0706     static inline CalcT dist_of_i2(CalcT const& dist_a1_i1)
0707     {
0708         CalcT const c2 = 2;
0709         CalcT const c4 = 4;
0710
0711         CalcT dist_a1_i2 = dist_a1_i1 - c2; // dist_a1_i2 = dist_a1_i1 - pi;
0712         if (dist_a1_i2 <= -c2)          // <= -pi
0713         {
0714             dist_a1_i2 += c4;           // += 2pi;
0715         }
0716         return dist_a1_i2;
0717     }
0718
0719     template <typename CalcT>
0720     static inline CalcT abs_distance(CalcT const& dist_a1_a2, CalcT const& dist_a1_i1)
0721     {
0722         if (dist_a1_i1 < CalcT(0))
0723             return -dist_a1_i1;
0724         else if (dist_a1_i1 > dist_a1_a2)
0725             return dist_a1_i1 - dist_a1_a2;
0726         else
0727             return CalcT(0);
0728     }
0729
0730     template <typename CalcT>
0731     static inline bool is_potentially_crossing(CalcT const& dist_a1_a2, CalcT const& dist_a1_i1, // in
0732                                                bool& is_on_a, bool& is_near_a1, bool& is_near_a2) // out
0733     {
0734         is_on_a = segment_ratio<CalcT>(dist_a1_i1, dist_a1_a2).on_segment();
0735         is_near_a1 = is_near(dist_a1_i1);
0736         is_near_a2 = is_near(dist_a1_a2 - dist_a1_i1);
0737         return is_on_a || is_near_a1 || is_near_a2;
0738     }
0739
0740     template <typename CalcT, typename P1, typename P2>
0741     static inline bool is_endpoint_equal(CalcT const& dist,
0742                                          P1 const& ai, P2 const& b1)
0743     {
0744         static CalcT const c0 = 0;
0745         return is_near(dist) && (math::equals(dist, c0) || equals_point_point(ai, b1));
0746     }
0747
0748     template <typename CalcT>
0749     static inline bool is_near(CalcT const& dist)
0750     {
0751         CalcT const small_number = CalcT(std::is_same<CalcT, float>::value ? 0.0001 : 0.00000001);
0752         return math::abs(dist) <= small_number;
0753     }
0754
0755     template <typename ProjCoord1, typename ProjCoord2>
0756     static inline int position_value(ProjCoord1 const& ca1,
0757                                      ProjCoord2 const& cb1,
0758                                      ProjCoord2 const& cb2)
0759     {
0760         // S1x  0   1    2     3   4
0761         // S2       |---------->
0762         return math::equals(ca1, cb1) ? 1
0763              : math::equals(ca1, cb2) ? 3
0764              : cb1 < cb2 ?
0765                 ( ca1 < cb1 ? 0
0766                 : ca1 > cb2 ? 4
0767                 : 2 )
0768               : ( ca1 > cb1 ? 0
0769                 : ca1 < cb2 ? 4
0770                 : 2 );
0771     }
0772
0773     template <typename Point1, typename Point2>
0774     static inline bool equals_point_point(Point1 const& point1, Point2 const& point2)
0775     {
0776         return strategy::within::spherical_point_point::apply(point1, point2);
0777     }
0778 };
0779
0780 struct spherical_segments_calc_policy
0781 {
0782     template <typename Point, typename Point3d>
0783     static Point from_cart3d(Point3d const& point_3d)
0784     {
0785         return formula::cart3d_to_sph<Point>(point_3d);
0786     }
0787
0788     template <typename Point3d, typename Point>
0789     static Point3d to_cart3d(Point const& point)
0790     {
0791         return formula::sph_to_cart3d<Point3d>(point);
0792     }
0793
0794     template <typename Point3d>
0795     struct plane
0796     {
0797         typedef typename coordinate_type<Point3d>::type coord_t;
0798
0799         // not normalized
0800         plane(Point3d const& p1, Point3d const& p2)
0801             : normal(cross_product(p1, p2))
0802         {}
0803
0804         int side_value(Point3d const& pt) const
0805         {
0806             return formula::sph_side_value(normal, pt);
0807         }
0808
0809         static coord_t cos_angle_between(Point3d const& p1, Point3d const& p2)
0810         {
0811             return dot_product(p1, p2);
0812         }
0813
0814         coord_t cos_angle_between(Point3d const& p1, Point3d const& p2, bool & is_forward) const
0815         {
0816             coord_t const c0 = 0;
0817             is_forward = dot_product(normal, cross_product(p1, p2)) >= c0;
0818             return dot_product(p1, p2);
0819         }
0820
0821         Point3d normal;
0822     };
0823
0824     template <typename Point3d>
0825     static plane<Point3d> get_plane(Point3d const& p1, Point3d const& p2)
0826     {
0827         return plane<Point3d>(p1, p2);
0828     }
0829
0830     template <typename Point3d>
0831     static bool intersection_points(plane<Point3d> const& plane1,
0832                                     plane<Point3d> const& plane2,
0833                                     Point3d & ip1, Point3d & ip2)
0834     {
0835         typedef typename coordinate_type<Point3d>::type coord_t;
0836
0837         ip1 = cross_product(plane1.normal, plane2.normal);
0838         // NOTE: the length should be greater than 0 at this point
0839         //       if the normals were not normalized and their dot product
0840         //       not checked before this function is called the length
0841         //       should be checked here (math::equals(len, c0))
0842         coord_t const len = math::sqrt(dot_product(ip1, ip1));
0843         geometry::detail::for_each_dimension<Point3d>([&](auto index)
0844         {
0845             coord_t const coord = get<index>(ip1) / len; // normalize
0846             set<index>(ip1, coord);
0847             set<index>(ip2, -coord);
0848         });
0849
0850         return true;
0851     }
0852 };
0853
0854
0855 template
0856 <
0857     typename CalculationType = void
0858 >
0859 struct spherical_segments
0860     : ecef_segments
0861         <
0862             spherical_segments_calc_policy,
0863             CalculationType
0864         >
0865 {};
0866
0867
0868 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
0869 namespace services
0870 {
0871
0872 /*template <typename CalculationType>
0873 struct default_strategy<spherical_polar_tag, CalculationType>
0874 {
0875     typedef spherical_segments<CalculationType> type;
0876 };*/
0877
0878 template <typename CalculationType>
0879 struct default_strategy<spherical_equatorial_tag, CalculationType>
0880 {
0881     typedef spherical_segments<CalculationType> type;
0882 };
0883
0884 template <typename CalculationType>
0885 struct default_strategy<geographic_tag, CalculationType>
0886 {
0887     // NOTE: Spherical strategy returns the same result as the geographic one
0888     // representing segments as great elliptic arcs. If the elliptic arcs are
0889     // not great elliptic arcs (the origin not in the center of the coordinate
0890     // system) then there may be problems with consistency of the side and
0891     // intersection strategies.
0892     typedef spherical_segments<CalculationType> type;
0893 };
0894
0895 } // namespace services
0896 #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
0897
0898
0899 }} // namespace strategy::intersection
0900
0901
0902 namespace strategy
0903 {
0904
0905 namespace within { namespace services
0906 {
0907
0908 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
0909 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, linear_tag, spherical_tag, spherical_tag>
0910 {
0911     typedef strategy::intersection::spherical_segments<> type;
0912 };
0913
0914 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
0915 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, polygonal_tag, spherical_tag, spherical_tag>
0916 {
0917     typedef strategy::intersection::spherical_segments<> type;
0918 };
0919
0920 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
0921 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, linear_tag, spherical_tag, spherical_tag>
0922 {
0923     typedef strategy::intersection::spherical_segments<> type;
0924 };
0925
0926 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
0927 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, polygonal_tag, spherical_tag, spherical_tag>
0928 {
0929     typedef strategy::intersection::spherical_segments<> type;
0930 };
0931
0932 }} // within::services
0933
0934 namespace covered_by { namespace services
0935 {
0936
0937 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
0938 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, linear_tag, spherical_tag, spherical_tag>
0939 {
0940     typedef strategy::intersection::spherical_segments<> type;
0941 };
0942
0943 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
0944 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, polygonal_tag, spherical_tag, spherical_tag>
0945 {
0946     typedef strategy::intersection::spherical_segments<> type;
0947 };
0948
0949 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
0950 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, linear_tag, spherical_tag, spherical_tag>
0951 {
0952     typedef strategy::intersection::spherical_segments<> type;
0953 };
0954
0955 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
0956 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, polygonal_tag, spherical_tag, spherical_tag>
0957 {
0958     typedef strategy::intersection::spherical_segments<> type;
0959 };
0960
0961 }} // within::services
0962
0963 } // strategy
0964
0965
0966 }} // namespace boost::geometry
0967
0968
0969 #endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_INTERSECTION_HPP