diff options
Diffstat (limited to 'boost/geometry/strategies/spherical/intersection.hpp')
-rw-r--r-- | boost/geometry/strategies/spherical/intersection.hpp | 701 |
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 |