summaryrefslogtreecommitdiff
path: root/boost/geometry/strategies/spherical/intersection.hpp
diff options
context:
space:
mode:
Diffstat (limited to 'boost/geometry/strategies/spherical/intersection.hpp')
-rw-r--r--boost/geometry/strategies/spherical/intersection.hpp701
1 files changed, 701 insertions, 0 deletions
diff --git a/boost/geometry/strategies/spherical/intersection.hpp b/boost/geometry/strategies/spherical/intersection.hpp
new file mode 100644
index 0000000000..4ffc853aad
--- /dev/null
+++ b/boost/geometry/strategies/spherical/intersection.hpp
@@ -0,0 +1,701 @@
+// Boost.Geometry
+
+// Copyright (c) 2016, Oracle and/or its affiliates.
+// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
+
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_INTERSECTION_HPP
+#define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_INTERSECTION_HPP
+
+#include <algorithm>
+
+#include <boost/geometry/core/cs.hpp>
+#include <boost/geometry/core/access.hpp>
+#include <boost/geometry/core/radian_access.hpp>
+#include <boost/geometry/core/tags.hpp>
+
+#include <boost/geometry/algorithms/detail/assign_values.hpp>
+#include <boost/geometry/algorithms/detail/assign_indexed_point.hpp>
+#include <boost/geometry/algorithms/detail/equals/point_point.hpp>
+#include <boost/geometry/algorithms/detail/recalculate.hpp>
+
+#include <boost/geometry/arithmetic/arithmetic.hpp>
+#include <boost/geometry/arithmetic/cross_product.hpp>
+#include <boost/geometry/arithmetic/dot_product.hpp>
+#include <boost/geometry/formulas/spherical.hpp>
+
+#include <boost/geometry/geometries/concepts/point_concept.hpp>
+#include <boost/geometry/geometries/concepts/segment_concept.hpp>
+
+#include <boost/geometry/policies/robustness/segment_ratio.hpp>
+
+#include <boost/geometry/strategies/side_info.hpp>
+#include <boost/geometry/strategies/intersection.hpp>
+#include <boost/geometry/strategies/intersection_result.hpp>
+
+#include <boost/geometry/util/math.hpp>
+#include <boost/geometry/util/select_calculation_type.hpp>
+
+
+namespace boost { namespace geometry
+{
+
+namespace strategy { namespace intersection
+{
+
+// NOTE:
+// The coordinates of crossing IP may be calculated with small precision in some cases.
+// For double, near the equator noticed error ~1e-9 so far greater than
+// machine epsilon which is ~1e-16. This error is ~0.04m.
+// E.g. consider two cases, one near the origin and the second one rotated by 90 deg around Z or SN axis.
+// After the conversion from spherical degrees to cartesian 3d the following coordinates
+// are calculated:
+// for sph (-1 -1, 1 1) deg cart3d ys are -0.017449748351250485 and 0.017449748351250485
+// for sph (89 -1, 91 1) deg cart3d xs are 0.017449748351250571 and -0.017449748351250450
+// During the conversion degrees must first be converted to radians and then radians
+// are passed into trigonometric functions. The error may have several causes:
+// 1. Radians cannot represent exactly the same angles as degrees.
+// 2. Different longitudes are passed into sin() for x, corresponding to cos() for y,
+// and for different angle the error of the result may be different.
+// 3. These non-corresponding cartesian coordinates are used in calculation,
+// e.g. multiplied several times in cross and dot products.
+// If it was a problem this strategy could e.g. "normalize" longitudes before the conversion using the source units
+// by rotating the globe around Z axis, so moving longitudes always the same way towards the origin,
+// assuming this could help which is not clear.
+// For now, intersection points near the endpoints are checked explicitly if needed (if the IP is near the endpoint)
+// to generate precise result for them. Only the crossing (i) case may suffer from lower precision.
+
+template <typename Policy, typename CalculationType = void>
+struct relate_spherical_segments
+{
+ typedef typename Policy::return_type return_type;
+
+ enum intersection_point_flag { ipi_inters = 0, ipi_at_a1, ipi_at_a2, ipi_at_b1, ipi_at_b2 };
+
+ template <typename CoordinateType, typename SegmentRatio, typename Vector3d>
+ struct segment_intersection_info
+ {
+ typedef typename select_most_precise
+ <
+ CoordinateType, double
+ >::type promoted_type;
+
+ promoted_type comparable_length_a() const
+ {
+ return robust_ra.denominator();
+ }
+
+ promoted_type comparable_length_b() const
+ {
+ return robust_rb.denominator();
+ }
+
+ template <typename Point, typename Segment1, typename Segment2>
+ void assign_a(Point& point, Segment1 const& a, Segment2 const& b) const
+ {
+ assign(point, a, b);
+ }
+ template <typename Point, typename Segment1, typename Segment2>
+ void assign_b(Point& point, Segment1 const& a, Segment2 const& b) const
+ {
+ assign(point, a, b);
+ }
+
+ template <typename Point, typename Segment1, typename Segment2>
+ void assign(Point& point, Segment1 const& a, Segment2 const& b) const
+ {
+ if (ip_flag == ipi_inters)
+ {
+ // TODO: assign the rest of coordinates
+ point = formula::cart3d_to_sph<Point>(intersection_point);
+ }
+ else if (ip_flag == ipi_at_a1)
+ {
+ detail::assign_point_from_index<0>(a, point);
+ }
+ else if (ip_flag == ipi_at_a2)
+ {
+ detail::assign_point_from_index<1>(a, point);
+ }
+ else if (ip_flag == ipi_at_b1)
+ {
+ detail::assign_point_from_index<0>(b, point);
+ }
+ else // ip_flag == ipi_at_b2
+ {
+ detail::assign_point_from_index<1>(b, point);
+ }
+ }
+
+ Vector3d intersection_point;
+ SegmentRatio robust_ra;
+ SegmentRatio robust_rb;
+ intersection_point_flag ip_flag;
+ };
+
+ // Relate segments a and b
+ template <typename Segment1, typename Segment2, typename RobustPolicy>
+ static inline return_type apply(Segment1 const& a, Segment2 const& b,
+ RobustPolicy const& robust_policy)
+ {
+ typedef typename point_type<Segment1>::type point1_t;
+ typedef typename point_type<Segment2>::type point2_t;
+ point1_t a1, a2;
+ point2_t b1, b2;
+
+ // TODO: use indexed_point_view if possible?
+ detail::assign_point_from_index<0>(a, a1);
+ detail::assign_point_from_index<1>(a, a2);
+ detail::assign_point_from_index<0>(b, b1);
+ detail::assign_point_from_index<1>(b, b2);
+
+ return apply(a, b, robust_policy, a1, a2, b1, b2);
+ }
+
+ // Relate segments a and b
+ template <typename Segment1, typename Segment2, typename RobustPolicy, typename Point1, typename Point2>
+ static inline return_type apply(Segment1 const& a, Segment2 const& b,
+ RobustPolicy const&,
+ Point1 const& a1, Point1 const& a2, Point2 const& b1, Point2 const& b2)
+ {
+ BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment1>) );
+ BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment2>) );
+
+ // TODO: check only 2 first coordinates here?
+ using geometry::detail::equals::equals_point_point;
+ bool a_is_point = equals_point_point(a1, a2);
+ bool b_is_point = equals_point_point(b1, b2);
+
+ if(a_is_point && b_is_point)
+ {
+ return equals_point_point(a1, b2)
+ ? Policy::degenerate(a, true)
+ : Policy::disjoint()
+ ;
+ }
+
+ typedef typename select_calculation_type
+ <Segment1, Segment2, CalculationType>::type calc_t;
+
+ calc_t const c0 = 0;
+ calc_t const c1 = 1;
+
+ typedef model::point<calc_t, 3, cs::cartesian> vec3d_t;
+
+ using namespace formula;
+ vec3d_t const a1v = sph_to_cart3d<vec3d_t>(a1);
+ vec3d_t const a2v = sph_to_cart3d<vec3d_t>(a2);
+ vec3d_t const b1v = sph_to_cart3d<vec3d_t>(b1);
+ vec3d_t const b2v = sph_to_cart3d<vec3d_t>(b2);
+
+ vec3d_t norm1 = cross_product(a1v, a2v);
+ vec3d_t norm2 = cross_product(b1v, b2v);
+
+ side_info sides;
+ // not normalized normals, the same as in SSF
+ sides.set<0>(sph_side_value(norm2, a1v), sph_side_value(norm2, a2v));
+ if (sides.same<0>())
+ {
+ // Both points are at same side of other segment, we can leave
+ return Policy::disjoint();
+ }
+ // not normalized normals, the same as in SSF
+ sides.set<1>(sph_side_value(norm1, b1v), sph_side_value(norm1, b2v));
+ if (sides.same<1>())
+ {
+ // Both points are at same side of other segment, we can leave
+ return Policy::disjoint();
+ }
+
+ // NOTE: at this point the segments may still be disjoint
+
+ bool collinear = sides.collinear();
+
+ calc_t const len1 = math::sqrt(dot_product(norm1, norm1));
+ calc_t const len2 = math::sqrt(dot_product(norm2, norm2));
+
+ // point or opposite sides of a sphere, assume point
+ if (math::equals(len1, c0))
+ {
+ a_is_point = true;
+ if (sides.get<0, 0>() == 0 || sides.get<0, 1>() == 0)
+ {
+ sides.set<0>(0, 0);
+ }
+ }
+ else
+ {
+ // normalize
+ divide_value(norm1, len1);
+ }
+
+ if (math::equals(len2, c0))
+ {
+ b_is_point = true;
+ if (sides.get<1, 0>() == 0 || sides.get<1, 1>() == 0)
+ {
+ sides.set<1>(0, 0);
+ }
+ }
+ else
+ {
+ // normalize
+ divide_value(norm2, len2);
+ }
+
+ // check both degenerated once more
+ if (a_is_point && b_is_point)
+ {
+ return equals_point_point(a1, b2)
+ ? Policy::degenerate(a, true)
+ : Policy::disjoint()
+ ;
+ }
+
+ // NOTE: at this point one of the segments may be degenerated
+ // and the segments may still be disjoint
+
+ calc_t dot_n1n2 = dot_product(norm1, norm2);
+
+ // NOTE: this is technically not needed since theoretically above sides
+ // are calculated, but just in case check the normals.
+ // Have in mind that SSF side strategy doesn't check this.
+ // collinear if normals are equal or opposite: cos(a) in {-1, 1}
+ if (!collinear && math::equals(math::abs(dot_n1n2), c1))
+ {
+ collinear = true;
+ sides.set<0>(0, 0);
+ sides.set<1>(0, 0);
+ }
+
+ if (collinear)
+ {
+ if (a_is_point)
+ {
+ return collinear_one_degenerted<calc_t>(a, true, b1, b2, a1, a2, b1v, b2v, norm2, a1v);
+ }
+ else if (b_is_point)
+ {
+ // b2 used to be consistent with (degenerated) checks above (is it needed?)
+ return collinear_one_degenerted<calc_t>(b, false, a1, a2, b1, b2, a1v, a2v, norm1, b1v);
+ }
+ else
+ {
+ calc_t dist_a1_a2, dist_a1_b1, dist_a1_b2;
+ calc_t dist_b1_b2, dist_b1_a1, dist_b1_a2;
+ // use shorter segment
+ if (len1 <= len2)
+ {
+ calculate_collinear_data(a1, a2, b1, b2, a1v, a2v, norm1, b1v, dist_a1_a2, dist_a1_b1);
+ calculate_collinear_data(a1, a2, b1, b2, a1v, a2v, norm1, b2v, dist_a1_a2, dist_a1_b2);
+ dist_b1_b2 = dist_a1_b2 - dist_a1_b1;
+ dist_b1_a1 = -dist_a1_b1;
+ dist_b1_a2 = dist_a1_a2 - dist_a1_b1;
+ }
+ else
+ {
+ calculate_collinear_data(b1, b2, a1, a2, b1v, b2v, norm2, a1v, dist_b1_b2, dist_b1_a1);
+ calculate_collinear_data(b1, b2, a1, a2, b1v, b2v, norm2, a2v, dist_b1_b2, dist_b1_a2);
+ dist_a1_a2 = dist_b1_a2 - dist_b1_a1;
+ dist_a1_b1 = -dist_b1_a1;
+ dist_a1_b2 = dist_b1_b2 - dist_b1_a1;
+ }
+
+ segment_ratio<calc_t> ra_from(dist_b1_a1, dist_b1_b2);
+ segment_ratio<calc_t> ra_to(dist_b1_a2, dist_b1_b2);
+ segment_ratio<calc_t> rb_from(dist_a1_b1, dist_a1_a2);
+ segment_ratio<calc_t> rb_to(dist_a1_b2, dist_a1_a2);
+
+ // NOTE: this is probably not needed
+ int const a1_wrt_b = position_value(c0, dist_a1_b1, dist_a1_b2);
+ int const a2_wrt_b = position_value(dist_a1_a2, dist_a1_b1, dist_a1_b2);
+ int const b1_wrt_a = position_value(c0, dist_b1_a1, dist_b1_a2);
+ int const b2_wrt_a = position_value(dist_b1_b2, dist_b1_a1, dist_b1_a2);
+
+ if (a1_wrt_b == 1)
+ {
+ ra_from.assign(0, dist_b1_b2);
+ rb_from.assign(0, dist_a1_a2);
+ }
+ else if (a1_wrt_b == 3)
+ {
+ ra_from.assign(dist_b1_b2, dist_b1_b2);
+ rb_to.assign(0, dist_a1_a2);
+ }
+
+ if (a2_wrt_b == 1)
+ {
+ ra_to.assign(0, dist_b1_b2);
+ rb_from.assign(dist_a1_a2, dist_a1_a2);
+ }
+ else if (a2_wrt_b == 3)
+ {
+ ra_to.assign(dist_b1_b2, dist_b1_b2);
+ rb_to.assign(dist_a1_a2, dist_a1_a2);
+ }
+
+ if ((a1_wrt_b < 1 && a2_wrt_b < 1) || (a1_wrt_b > 3 && a2_wrt_b > 3))
+ {
+ return Policy::disjoint();
+ }
+
+ bool const opposite = dot_n1n2 < c0;
+
+ return Policy::segments_collinear(a, b, opposite,
+ a1_wrt_b, a2_wrt_b, b1_wrt_a, b2_wrt_a,
+ ra_from, ra_to, rb_from, rb_to);
+ }
+ }
+ else // crossing
+ {
+ if (a_is_point || b_is_point)
+ {
+ return Policy::disjoint();
+ }
+
+ vec3d_t i1;
+ intersection_point_flag ip_flag;
+ calc_t dist_a1_a2, dist_a1_i1, dist_b1_b2, dist_b1_i1;
+ if (calculate_ip_data(a1, a2, b1, b2, a1v, a2v, b1v, b2v, norm1, norm2, sides,
+ i1, dist_a1_a2, dist_a1_i1, dist_b1_b2, dist_b1_i1, ip_flag))
+ {
+ // intersects
+ segment_intersection_info
+ <
+ calc_t,
+ segment_ratio<calc_t>,
+ vec3d_t
+ > sinfo;
+
+ sinfo.robust_ra.assign(dist_a1_i1, dist_a1_a2);
+ sinfo.robust_rb.assign(dist_b1_i1, dist_b1_b2);
+ sinfo.intersection_point = i1;
+ sinfo.ip_flag = ip_flag;
+
+ return Policy::segments_crosses(sides, sinfo, a, b);
+ }
+ else
+ {
+ return Policy::disjoint();
+ }
+ }
+ }
+
+private:
+ template <typename CalcT, typename Segment, typename Point1, typename Point2, typename Vec3d>
+ static inline return_type collinear_one_degenerted(Segment const& segment, bool degenerated_a,
+ Point1 const& a1, Point1 const& a2,
+ Point2 const& b1, Point2 const& b2,
+ Vec3d const& v1, Vec3d const& v2, Vec3d const& norm,
+ Vec3d const& vother)
+ {
+ CalcT dist_1_2, dist_1_o;
+ return ! calculate_collinear_data(a1, a2, b1, b2, v1, v2, norm, vother, dist_1_2, dist_1_o)
+ ? Policy::disjoint()
+ : Policy::one_degenerate(segment, segment_ratio<CalcT>(dist_1_o, dist_1_2), degenerated_a);
+ }
+
+ template <typename Point1, typename Point2, typename Vec3d, typename CalcT>
+ static inline bool calculate_collinear_data(Point1 const& a1, Point1 const& a2,
+ Point2 const& b1, Point2 const& b2,
+ Vec3d const& a1v, // in
+ Vec3d const& a2v, // in
+ Vec3d const& norm1, // in
+ Vec3d const& b1v_or_b2v, // in
+ CalcT& dist_a1_a2, CalcT& dist_a1_i1) // out
+ {
+ // calculate dist_a1_a2 and dist_a1_i1
+ calculate_dists(a1v, a2v, norm1, b1v_or_b2v, dist_a1_a2, dist_a1_i1);
+
+ // if i1 is close to a1 and b1 or b2 is equal to a1
+ if (is_endpoint_equal(dist_a1_i1, a1, b1, b2))
+ {
+ dist_a1_i1 = 0;
+ return true;
+ }
+ // or i1 is close to a2 and b1 or b2 is equal to a2
+ else if (is_endpoint_equal(dist_a1_a2 - dist_a1_i1, a2, b1, b2))
+ {
+ dist_a1_i1 = dist_a1_a2;
+ return true;
+ }
+
+ // or i1 is on b
+ return segment_ratio<CalcT>(dist_a1_i1, dist_a1_a2).on_segment();
+ }
+
+ template <typename Point1, typename Point2, typename Vec3d, typename CalcT>
+ static inline bool calculate_ip_data(Point1 const& a1, Point1 const& a2, // in
+ Point2 const& b1, Point2 const& b2, // in
+ Vec3d const& a1v, Vec3d const& a2v, // in
+ Vec3d const& b1v, Vec3d const& b2v, // in
+ Vec3d const& norm1, Vec3d const& norm2, // in
+ side_info const& sides, // in
+ Vec3d & i1, // in-out
+ CalcT& dist_a1_a2, CalcT& dist_a1_i1, // out
+ CalcT& dist_b1_b2, CalcT& dist_b1_i1, // out
+ intersection_point_flag& ip_flag) // out
+ {
+ // great circles intersections
+ i1 = cross_product(norm1, norm2);
+ // NOTE: the length should be greater than 0 at this point
+ // if the normals were not normalized and their dot product
+ // not checked before this function is called the length
+ // should be checked here (math::equals(len, c0))
+ CalcT const len = math::sqrt(dot_product(i1, i1));
+ divide_value(i1, len); // normalize i1
+
+ calculate_dists(a1v, a2v, norm1, i1, dist_a1_a2, dist_a1_i1);
+
+ // choose the opposite side of the globe if the distance is shorter
+ {
+ CalcT const d = abs_distance(dist_a1_a2, dist_a1_i1);
+ if (d > CalcT(0))
+ {
+ CalcT const dist_a1_i2 = dist_of_i2(dist_a1_i1);
+ CalcT const d2 = abs_distance(dist_a1_a2, dist_a1_i2);
+ if (d2 < d)
+ {
+ dist_a1_i1 = dist_a1_i2;
+ multiply_value(i1, CalcT(-1)); // the opposite intersection
+ }
+ }
+ }
+
+ bool is_on_a = false, is_near_a1 = false, is_near_a2 = false;
+ if (! is_potentially_crossing(dist_a1_a2, dist_a1_i1, is_on_a, is_near_a1, is_near_a2))
+ {
+ return false;
+ }
+
+ calculate_dists(b1v, b2v, norm2, i1, dist_b1_b2, dist_b1_i1);
+
+ bool is_on_b = false, is_near_b1 = false, is_near_b2 = false;
+ if (! is_potentially_crossing(dist_b1_b2, dist_b1_i1, is_on_b, is_near_b1, is_near_b2))
+ {
+ return false;
+ }
+
+ // reassign the IP if some endpoints overlap
+ using geometry::detail::equals::equals_point_point;
+ if (is_near_a1)
+ {
+ if (is_near_b1 && equals_point_point(a1, b1))
+ {
+ dist_a1_i1 = 0;
+ dist_b1_i1 = 0;
+ //i1 = a1v;
+ ip_flag = ipi_at_a1;
+ return true;
+ }
+
+ if (is_near_b2 && equals_point_point(a1, b2))
+ {
+ dist_a1_i1 = 0;
+ dist_b1_i1 = dist_b1_b2;
+ //i1 = a1v;
+ ip_flag = ipi_at_a1;
+ return true;
+ }
+ }
+
+ if (is_near_a2)
+ {
+ if (is_near_b1 && equals_point_point(a2, b1))
+ {
+ dist_a1_i1 = dist_a1_a2;
+ dist_b1_i1 = 0;
+ //i1 = a2v;
+ ip_flag = ipi_at_a2;
+ return true;
+ }
+
+ if (is_near_b2 && equals_point_point(a2, b2))
+ {
+ dist_a1_i1 = dist_a1_a2;
+ dist_b1_i1 = dist_b1_b2;
+ //i1 = a2v;
+ ip_flag = ipi_at_a2;
+ return true;
+ }
+ }
+
+ // at this point we know that the endpoints doesn't overlap
+ // reassign IP and distance if the IP is on a segment and one of
+ // the endpoints of the other segment lies on the former segment
+ if (is_on_a)
+ {
+ if (is_near_b1 && sides.template get<1, 0>() == 0) // b1 wrt a
+ {
+ dist_b1_i1 = 0;
+ //i1 = b1v;
+ ip_flag = ipi_at_b1;
+ return true;
+ }
+
+ if (is_near_b2 && sides.template get<1, 1>() == 0) // b2 wrt a
+ {
+ dist_b1_i1 = dist_b1_b2;
+ //i1 = b2v;
+ ip_flag = ipi_at_b2;
+ return true;
+ }
+ }
+
+ if (is_on_b)
+ {
+ if (is_near_a1 && sides.template get<0, 0>() == 0) // a1 wrt b
+ {
+ dist_a1_i1 = 0;
+ //i1 = a1v;
+ ip_flag = ipi_at_a1;
+ return true;
+ }
+
+ if (is_near_a2 && sides.template get<0, 1>() == 0) // a2 wrt b
+ {
+ dist_a1_i1 = dist_a1_a2;
+ //i1 = a2v;
+ ip_flag = ipi_at_a2;
+ return true;
+ }
+ }
+
+ ip_flag = ipi_inters;
+
+ return is_on_a && is_on_b;
+ }
+
+ template <typename Vec3d, typename CalcT>
+ static inline void calculate_dists(Vec3d const& a1v, // in
+ Vec3d const& a2v, // in
+ Vec3d const& norm1, // in
+ Vec3d const& i1, // in
+ CalcT& dist_a1_a2, CalcT& dist_a1_i1) // out
+ {
+ CalcT const c0 = 0;
+ CalcT const c1 = 1;
+ CalcT const c2 = 2;
+ CalcT const c4 = 4;
+
+ CalcT cos_a1_a2 = dot_product(a1v, a2v);
+ dist_a1_a2 = -cos_a1_a2 + c1; // [1, -1] -> [0, 2] representing [0, pi]
+
+ CalcT cos_a1_i1 = dot_product(a1v, i1);
+ dist_a1_i1 = -cos_a1_i1 + c1; // [0, 2] representing [0, pi]
+ if (dot_product(norm1, cross_product(a1v, i1)) < c0) // left or right of a1 on a
+ {
+ dist_a1_i1 = -dist_a1_i1; // [0, 2] -> [0, -2] representing [0, -pi]
+ }
+ if (dist_a1_i1 <= -c2) // <= -pi
+ {
+ dist_a1_i1 += c4; // += 2pi
+ }
+ }
+
+ // the dist of the ip on the other side of the sphere
+ template <typename CalcT>
+ static inline CalcT dist_of_i2(CalcT const& dist_a1_i1)
+ {
+ CalcT const c2 = 2;
+ CalcT const c4 = 4;
+
+ CalcT dist_a1_i2 = dist_a1_i1 - c2; // dist_a1_i2 = dist_a1_i1 - pi;
+ if (dist_a1_i2 <= -c2) // <= -pi
+ {
+ dist_a1_i2 += c4; // += 2pi;
+ }
+ return dist_a1_i2;
+ }
+
+ template <typename CalcT>
+ static inline CalcT abs_distance(CalcT const& dist_a1_a2, CalcT const& dist_a1_i1)
+ {
+ if (dist_a1_i1 < CalcT(0))
+ return -dist_a1_i1;
+ else if (dist_a1_i1 > dist_a1_a2)
+ return dist_a1_i1 - dist_a1_a2;
+ else
+ return CalcT(0);
+ }
+
+ template <typename CalcT>
+ static inline bool is_potentially_crossing(CalcT const& dist_a1_a2, CalcT const& dist_a1_i1, // in
+ bool& is_on_a, bool& is_near_a1, bool& is_near_a2) // out
+ {
+ is_on_a = segment_ratio<CalcT>(dist_a1_i1, dist_a1_a2).on_segment();
+ is_near_a1 = is_near(dist_a1_i1);
+ is_near_a2 = is_near(dist_a1_a2 - dist_a1_i1);
+ return is_on_a || is_near_a1 || is_near_a2;
+ }
+
+ template <typename CalcT, typename P1, typename P2>
+ static inline bool is_endpoint_equal(CalcT const& dist,
+ P1 const& ai, P2 const& b1, P2 const& b2)
+ {
+ using geometry::detail::equals::equals_point_point;
+ return is_near(dist) && (equals_point_point(ai, b1) || equals_point_point(ai, b2));
+ }
+
+ template <typename CalcT>
+ static inline bool is_near(CalcT const& dist)
+ {
+ CalcT const small_number = CalcT(boost::is_same<CalcT, float>::value ? 0.0001 : 0.00000001);
+ return math::abs(dist) <= small_number;
+ }
+
+ template <typename ProjCoord1, typename ProjCoord2>
+ static inline int position_value(ProjCoord1 const& ca1,
+ ProjCoord2 const& cb1,
+ ProjCoord2 const& cb2)
+ {
+ // S1x 0 1 2 3 4
+ // S2 |---------->
+ return math::equals(ca1, cb1) ? 1
+ : math::equals(ca1, cb2) ? 3
+ : cb1 < cb2 ?
+ ( ca1 < cb1 ? 0
+ : ca1 > cb2 ? 4
+ : 2 )
+ : ( ca1 > cb1 ? 0
+ : ca1 < cb2 ? 4
+ : 2 );
+ }
+};
+
+
+#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
+namespace services
+{
+
+/*template <typename Policy, typename CalculationType>
+struct default_strategy<spherical_polar_tag, Policy, CalculationType>
+{
+ typedef relate_spherical_segments<Policy, CalculationType> type;
+};*/
+
+template <typename Policy, typename CalculationType>
+struct default_strategy<spherical_equatorial_tag, Policy, CalculationType>
+{
+ typedef relate_spherical_segments<Policy, CalculationType> type;
+};
+
+template <typename Policy, typename CalculationType>
+struct default_strategy<geographic_tag, Policy, CalculationType>
+{
+ typedef relate_spherical_segments<Policy, CalculationType> type;
+};
+
+} // namespace services
+#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
+
+
+}} // namespace strategy::intersection
+
+}} // namespace boost::geometry
+
+
+#endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_INTERSECTION_HPP