diff options
Diffstat (limited to 'boost/geometry/strategies/spherical/intersection.hpp')
| -rw-r--r-- | boost/geometry/strategies/spherical/intersection.hpp | 525 |
1 files changed, 402 insertions, 123 deletions
diff --git a/boost/geometry/strategies/spherical/intersection.hpp b/boost/geometry/strategies/spherical/intersection.hpp index 4ffc853aad..5d37583333 100644 --- a/boost/geometry/strategies/spherical/intersection.hpp +++ b/boost/geometry/strategies/spherical/intersection.hpp @@ -1,6 +1,6 @@ // Boost.Geometry -// Copyright (c) 2016, Oracle and/or its affiliates. +// Copyright (c) 2016-2017, 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, @@ -25,6 +25,7 @@ #include <boost/geometry/arithmetic/arithmetic.hpp> #include <boost/geometry/arithmetic/cross_product.hpp> #include <boost/geometry/arithmetic/dot_product.hpp> +#include <boost/geometry/arithmetic/normalize.hpp> #include <boost/geometry/formulas/spherical.hpp> #include <boost/geometry/geometries/concepts/point_concept.hpp> @@ -32,9 +33,16 @@ #include <boost/geometry/policies/robustness/segment_ratio.hpp> -#include <boost/geometry/strategies/side_info.hpp> +#include <boost/geometry/strategies/agnostic/point_in_poly_winding.hpp> +#include <boost/geometry/strategies/covered_by.hpp> #include <boost/geometry/strategies/intersection.hpp> #include <boost/geometry/strategies/intersection_result.hpp> +#include <boost/geometry/strategies/side.hpp> +#include <boost/geometry/strategies/side_info.hpp> +#include <boost/geometry/strategies/spherical/area.hpp> +#include <boost/geometry/strategies/spherical/distance_haversine.hpp> +#include <boost/geometry/strategies/spherical/ssf.hpp> +#include <boost/geometry/strategies/within.hpp> #include <boost/geometry/util/math.hpp> #include <boost/geometry/util/select_calculation_type.hpp> @@ -68,13 +76,80 @@ namespace strategy { namespace intersection // 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 +template +< + typename CalcPolicy, + typename CalculationType = void +> +struct ecef_segments { - typedef typename Policy::return_type return_type; + typedef side::spherical_side_formula<CalculationType> side_strategy_type; + + static inline side_strategy_type get_side_strategy() + { + return side_strategy_type(); + } + + template <typename Geometry1, typename Geometry2> + struct point_in_geometry_strategy + { + typedef strategy::within::winding + < + typename point_type<Geometry1>::type, + typename point_type<Geometry2>::type, + side_strategy_type, + CalculationType + > type; + }; + + template <typename Geometry1, typename Geometry2> + static inline typename point_in_geometry_strategy<Geometry1, Geometry2>::type + get_point_in_geometry_strategy() + { + typedef typename point_in_geometry_strategy + < + Geometry1, Geometry2 + >::type strategy_type; + return strategy_type(); + } + + template <typename Geometry> + struct area_strategy + { + typedef area::spherical + < + typename point_type<Geometry>::type, + CalculationType + > type; + }; + + template <typename Geometry> + static inline typename area_strategy<Geometry>::type get_area_strategy() + { + typedef typename area_strategy<Geometry>::type strategy_type; + return strategy_type(); + } + + template <typename Geometry> + struct distance_strategy + { + typedef distance::haversine + < + typename coordinate_type<Geometry>::type, + CalculationType + > type; + }; + + template <typename Geometry> + static inline typename distance_strategy<Geometry>::type get_distance_strategy() + { + typedef typename distance_strategy<Geometry>::type strategy_type; + return strategy_type(); + } enum intersection_point_flag { ipi_inters = 0, ipi_at_a1, ipi_at_a2, ipi_at_b1, ipi_at_b2 }; + // segment_intersection_info cannot outlive relate_ecef_segments template <typename CoordinateType, typename SegmentRatio, typename Vector3d> struct segment_intersection_info { @@ -83,6 +158,10 @@ struct relate_spherical_segments CoordinateType, double >::type promoted_type; + segment_intersection_info(CalcPolicy const& calc) + : calc_policy(calc) + {} + promoted_type comparable_length_a() const { return robust_ra.denominator(); @@ -110,7 +189,7 @@ struct relate_spherical_segments if (ip_flag == ipi_inters) { // TODO: assign the rest of coordinates - point = formula::cart3d_to_sph<Point>(intersection_point); + point = calc_policy.template from_cart3d<Point>(intersection_point); } else if (ip_flag == ipi_at_a1) { @@ -134,12 +213,21 @@ struct relate_spherical_segments SegmentRatio robust_ra; SegmentRatio robust_rb; intersection_point_flag ip_flag; + + CalcPolicy const& calc_policy; }; // 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) + template + < + typename Segment1, + typename Segment2, + typename Policy, + typename RobustPolicy + > + static inline typename Policy::return_type + apply(Segment1 const& a, Segment2 const& b, + Policy const& policy, RobustPolicy const& robust_policy) { typedef typename point_type<Segment1>::type point1_t; typedef typename point_type<Segment2>::type point2_t; @@ -152,15 +240,30 @@ struct relate_spherical_segments 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); + return apply(a, b, policy, 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) + template + < + typename Segment1, + typename Segment2, + typename Policy, + typename RobustPolicy, + typename Point1, + typename Point2 + > + static inline typename Policy::return_type + apply(Segment1 const& a, Segment2 const& b, + Policy const&, RobustPolicy const&, + Point1 const& a1, Point1 const& a2, Point2 const& b1, Point2 const& b2) { + // For now create it using default constructor. In the future it could + // be stored in strategy. However then apply() wouldn't be static and + // all relops and setops would have to take the strategy or model. + // Initialize explicitly to prevent compiler errors in case of PoD type + CalcPolicy const calc_policy = CalcPolicy(); + BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment1>) ); BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment2>) ); @@ -185,25 +288,29 @@ struct relate_spherical_segments 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 const a1v = calc_policy.template to_cart3d<vec3d_t>(a1); + vec3d_t const a2v = calc_policy.template to_cart3d<vec3d_t>(a2); + vec3d_t const b1v = calc_policy.template to_cart3d<vec3d_t>(b1); + vec3d_t const b2v = calc_policy.template 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)); + + typename CalcPolicy::template plane<vec3d_t> + plane2 = calc_policy.get_plane(b1v, b2v); + + // not normalized normals, the same as in side strategy + sides.set<0>(plane2.side_value(a1v), plane2.side_value(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)); + + typename CalcPolicy::template plane<vec3d_t> + plane1 = calc_policy.get_plane(a1v, a2v); + + // not normalized normals, the same as in side strategy + sides.set<1>(plane1.side_value(b1v), plane1.side_value(b2v)); if (sides.same<1>()) { // Both points are at same side of other segment, we can leave @@ -212,13 +319,10 @@ struct relate_spherical_segments // 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)); + calc_t len1, len2; - // point or opposite sides of a sphere, assume point - if (math::equals(len1, c0)) + // point or opposite sides of a sphere/spheroid, assume point + if (! detail::vec_normalize(plane1.normal, len1)) { a_is_point = true; if (sides.get<0, 0>() == 0 || sides.get<0, 1>() == 0) @@ -226,13 +330,8 @@ struct relate_spherical_segments sides.set<0>(0, 0); } } - else - { - // normalize - divide_value(norm1, len1); - } - if (math::equals(len2, c0)) + if (! detail::vec_normalize(plane2.normal, len2)) { b_is_point = true; if (sides.get<1, 0>() == 0 || sides.get<1, 1>() == 0) @@ -240,11 +339,6 @@ struct relate_spherical_segments sides.set<1>(0, 0); } } - else - { - // normalize - divide_value(norm2, len2); - } // check both degenerated once more if (a_is_point && b_is_point) @@ -255,16 +349,42 @@ struct relate_spherical_segments ; } + // NOTE: at this point the segments may still be 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); + bool collinear = sides.collinear(); + + if (! collinear) + { + // NOTE: for some approximations it's possible that both points may lie + // on the same geodesic but still some of the sides may be != 0. + // This is e.g. true for long segments represented as elliptic arcs + // with origin different than the center of the coordinate system. + // So make the sides consistent + + // WARNING: the side strategy doesn't have the info about the other + // segment so it may return results inconsistent with this intersection + // strategy, as it checks both segments for consistency + + if (sides.get<0, 0>() == 0 && sides.get<0, 1>() == 0) + { + collinear = true; + sides.set<1>(0, 0); + } + else if (sides.get<1, 0>() == 0 && sides.get<1, 1>() == 0) + { + collinear = true; + sides.set<0>(0, 0); + } + } + + calc_t dot_n1n2 = dot_product(plane1.normal, plane2.normal); // 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)) + if (! collinear && math::equals(math::abs(dot_n1n2), c1)) { collinear = true; sides.set<0>(0, 0); @@ -275,12 +395,12 @@ struct relate_spherical_segments { if (a_is_point) { - return collinear_one_degenerted<calc_t>(a, true, b1, b2, a1, a2, b1v, b2v, norm2, a1v); + return collinear_one_degenerated<Policy, calc_t>(a, true, b1, b2, a1, a2, b1v, b2v, plane2, 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); + return collinear_one_degenerated<Policy, calc_t>(b, false, a1, a2, b1, b2, a1v, a2v, plane1, b1v); } else { @@ -289,16 +409,16 @@ struct relate_spherical_segments // 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); + calculate_collinear_data(a1, a2, b1, b2, a1v, a2v, plane1, b1v, dist_a1_a2, dist_a1_b1); + calculate_collinear_data(a1, a2, b1, b2, a1v, a2v, plane1, 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); + calculate_collinear_data(b1, b2, a1, a2, b1v, b2v, plane2, a1v, dist_b1_b2, dist_b1_a1); + calculate_collinear_data(b1, b2, a1, a2, b1v, b2v, plane2, 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; @@ -359,7 +479,8 @@ struct relate_spherical_segments 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, + if (calculate_ip_data(a1, a2, b1, b2, a1v, a2v, b1v, b2v, + plane1, plane2, calc_policy, sides, i1, dist_a1_a2, dist_a1_i1, dist_b1_b2, dist_b1_i1, ip_flag)) { // intersects @@ -368,7 +489,7 @@ struct relate_spherical_segments calc_t, segment_ratio<calc_t>, vec3d_t - > sinfo; + > sinfo(calc_policy); sinfo.robust_ra.assign(dist_a1_i1, dist_a1_a2); sinfo.robust_rb.assign(dist_b1_i1, dist_b1_b2); @@ -385,30 +506,32 @@ struct relate_spherical_segments } 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) + template <typename Policy, typename CalcT, typename Segment, typename Point1, typename Point2, typename Vec3d, typename Plane> + static inline typename Policy::return_type + collinear_one_degenerated(Segment const& segment, bool degenerated_a, + Point1 const& a1, Point1 const& a2, + Point2 const& b1, Point2 const& b2, + Vec3d const& v1, Vec3d const& v2, + Plane const& plane, + 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) + return ! calculate_collinear_data(a1, a2, b1, b2, v1, v2, plane, 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 + template <typename Point1, typename Point2, typename Vec3d, typename Plane, typename CalcT> + static inline bool calculate_collinear_data(Point1 const& a1, Point1 const& a2, // in + Point2 const& b1, Point2 const& b2, // in + Vec3d const& a1v, // in + Vec3d const& a2v, // in + Plane const& plane1, // 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); + calculate_dists(a1v, a2v, plane1, 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)) @@ -427,54 +550,53 @@ private: return segment_ratio<CalcT>(dist_a1_i1, dist_a1_a2).on_segment(); } - template <typename Point1, typename Point2, typename Vec3d, typename CalcT> + template <typename Point1, typename Point2, typename Vec3d, typename Plane, 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 + Plane const& plane1, // in + Plane const& plane2, // in + CalcPolicy const& calc_policy, // in + side_info const& sides, // in + Vec3d & ip, // out + CalcT& dist_a1_a2, CalcT& dist_a1_ip, // out + CalcT& dist_b1_b2, CalcT& dist_b1_ip, // 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); + Vec3d ip1, ip2; + calc_policy.intersection_points(plane1, plane2, ip1, ip2); + calculate_dists(a1v, a2v, plane1, ip1, dist_a1_a2, dist_a1_ip); + ip = ip1; + // choose the opposite side of the globe if the distance is shorter { - CalcT const d = abs_distance(dist_a1_a2, dist_a1_i1); + CalcT const d = abs_distance(dist_a1_a2, dist_a1_ip); if (d > CalcT(0)) { - CalcT const dist_a1_i2 = dist_of_i2(dist_a1_i1); + // TODO: this should be ok not only for sphere + // but requires more investigation + CalcT const dist_a1_i2 = dist_of_i2(dist_a1_ip); 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 + dist_a1_ip = dist_a1_i2; + ip = ip2; } } } 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)) + if (! is_potentially_crossing(dist_a1_a2, dist_a1_ip, is_on_a, is_near_a1, is_near_a2)) { return false; } - calculate_dists(b1v, b2v, norm2, i1, dist_b1_b2, dist_b1_i1); + calculate_dists(b1v, b2v, plane2, ip, dist_b1_b2, dist_b1_ip); 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)) + if (! is_potentially_crossing(dist_b1_b2, dist_b1_ip, is_on_b, is_near_b1, is_near_b2)) { return false; } @@ -485,8 +607,8 @@ private: { if (is_near_b1 && equals_point_point(a1, b1)) { - dist_a1_i1 = 0; - dist_b1_i1 = 0; + dist_a1_ip = 0; + dist_b1_ip = 0; //i1 = a1v; ip_flag = ipi_at_a1; return true; @@ -494,8 +616,8 @@ private: if (is_near_b2 && equals_point_point(a1, b2)) { - dist_a1_i1 = 0; - dist_b1_i1 = dist_b1_b2; + dist_a1_ip = 0; + dist_b1_ip = dist_b1_b2; //i1 = a1v; ip_flag = ipi_at_a1; return true; @@ -506,8 +628,8 @@ private: { if (is_near_b1 && equals_point_point(a2, b1)) { - dist_a1_i1 = dist_a1_a2; - dist_b1_i1 = 0; + dist_a1_ip = dist_a1_a2; + dist_b1_ip = 0; //i1 = a2v; ip_flag = ipi_at_a2; return true; @@ -515,8 +637,8 @@ private: if (is_near_b2 && equals_point_point(a2, b2)) { - dist_a1_i1 = dist_a1_a2; - dist_b1_i1 = dist_b1_b2; + dist_a1_ip = dist_a1_a2; + dist_b1_ip = dist_b1_b2; //i1 = a2v; ip_flag = ipi_at_a2; return true; @@ -530,7 +652,7 @@ private: { if (is_near_b1 && sides.template get<1, 0>() == 0) // b1 wrt a { - dist_b1_i1 = 0; + dist_b1_ip = 0; //i1 = b1v; ip_flag = ipi_at_b1; return true; @@ -538,7 +660,7 @@ private: if (is_near_b2 && sides.template get<1, 1>() == 0) // b2 wrt a { - dist_b1_i1 = dist_b1_b2; + dist_b1_ip = dist_b1_b2; //i1 = b2v; ip_flag = ipi_at_b2; return true; @@ -549,7 +671,7 @@ private: { if (is_near_a1 && sides.template get<0, 0>() == 0) // a1 wrt b { - dist_a1_i1 = 0; + dist_a1_ip = 0; //i1 = a1v; ip_flag = ipi_at_a1; return true; @@ -557,7 +679,7 @@ private: if (is_near_a2 && sides.template get<0, 1>() == 0) // a2 wrt b { - dist_a1_i1 = dist_a1_a2; + dist_a1_ip = dist_a1_a2; //i1 = a2v; ip_flag = ipi_at_a2; return true; @@ -569,24 +691,26 @@ private: 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 + template <typename Vec3d, typename Plane, typename CalcT> + static inline void calculate_dists(Vec3d const& a1v, // in + Vec3d const& a2v, // in + Plane const& plane1, // in + Vec3d const& i1, // in + CalcT& dist_a1_a2, // out + CalcT& dist_a1_i1) // out { - CalcT const c0 = 0; + //CalcT const c0 = 0; CalcT const c1 = 1; CalcT const c2 = 2; CalcT const c4 = 4; - CalcT cos_a1_a2 = dot_product(a1v, a2v); + CalcT cos_a1_a2 = plane1.cos_angle_between(a1v, a2v); dist_a1_a2 = -cos_a1_a2 + c1; // [1, -1] -> [0, 2] representing [0, pi] - CalcT cos_a1_i1 = dot_product(a1v, i1); + bool is_forward = true; + CalcT cos_a1_i1 = plane1.cos_angle_between(a1v, i1, is_forward); 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 + if (! is_forward) // left or right of a1 on a { dist_a1_i1 = -dist_a1_i1; // [0, 2] -> [0, -2] representing [0, -pi] } @@ -666,27 +790,117 @@ private: } }; +struct spherical_segments_calc_policy +{ + template <typename Point, typename Point3d> + static Point from_cart3d(Point3d const& point_3d) + { + return formula::cart3d_to_sph<Point>(point_3d); + } + + template <typename Point3d, typename Point> + static Point3d to_cart3d(Point const& point) + { + return formula::sph_to_cart3d<Point3d>(point); + } + + template <typename Point3d> + struct plane + { + typedef typename coordinate_type<Point3d>::type coord_t; + + // not normalized + plane(Point3d const& p1, Point3d const& p2) + : normal(cross_product(p1, p2)) + {} + + int side_value(Point3d const& pt) const + { + return formula::sph_side_value(normal, pt); + } + + static coord_t cos_angle_between(Point3d const& p1, Point3d const& p2) + { + return dot_product(p1, p2); + } + + coord_t cos_angle_between(Point3d const& p1, Point3d const& p2, bool & is_forward) const + { + coord_t const c0 = 0; + is_forward = dot_product(normal, cross_product(p1, p2)) >= c0; + return dot_product(p1, p2); + } + + Point3d normal; + }; + + template <typename Point3d> + static plane<Point3d> get_plane(Point3d const& p1, Point3d const& p2) + { + return plane<Point3d>(p1, p2); + } + + template <typename Point3d> + static bool intersection_points(plane<Point3d> const& plane1, + plane<Point3d> const& plane2, + Point3d & ip1, Point3d & ip2) + { + typedef typename coordinate_type<Point3d>::type coord_t; + + ip1 = cross_product(plane1.normal, plane2.normal); + // 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)) + coord_t const len = math::sqrt(dot_product(ip1, ip1)); + divide_value(ip1, len); // normalize i1 + + ip2 = ip1; + multiply_value(ip2, coord_t(-1)); + + return true; + } +}; + + +template +< + typename CalculationType = void +> +struct spherical_segments + : ecef_segments + < + spherical_segments_calc_policy, + CalculationType + > +{}; + #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS namespace services { -/*template <typename Policy, typename CalculationType> -struct default_strategy<spherical_polar_tag, Policy, CalculationType> +/*template <typename CalculationType> +struct default_strategy<spherical_polar_tag, CalculationType> { - typedef relate_spherical_segments<Policy, CalculationType> type; + typedef spherical_segments<CalculationType> type; };*/ -template <typename Policy, typename CalculationType> -struct default_strategy<spherical_equatorial_tag, Policy, CalculationType> +template <typename CalculationType> +struct default_strategy<spherical_equatorial_tag, CalculationType> { - typedef relate_spherical_segments<Policy, CalculationType> type; + typedef spherical_segments<CalculationType> type; }; -template <typename Policy, typename CalculationType> -struct default_strategy<geographic_tag, Policy, CalculationType> +template <typename CalculationType> +struct default_strategy<geographic_tag, CalculationType> { - typedef relate_spherical_segments<Policy, CalculationType> type; + // NOTE: Spherical strategy returns the same result as the geographic one + // representing segments as great elliptic arcs. If the elliptic arcs are + // not great elliptic arcs (the origin not in the center of the coordinate + // system) then there may be problems with consistency of the side and + // intersection strategies. + typedef spherical_segments<CalculationType> type; }; } // namespace services @@ -695,6 +909,71 @@ struct default_strategy<geographic_tag, Policy, CalculationType> }} // namespace strategy::intersection + +namespace strategy +{ + +namespace within { namespace services +{ + +template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> +struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, linear_tag, spherical_tag, spherical_tag> +{ + typedef strategy::intersection::spherical_segments<> type; +}; + +template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> +struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, polygonal_tag, spherical_tag, spherical_tag> +{ + typedef strategy::intersection::spherical_segments<> type; +}; + +template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> +struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, linear_tag, spherical_tag, spherical_tag> +{ + typedef strategy::intersection::spherical_segments<> type; +}; + +template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> +struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, polygonal_tag, spherical_tag, spherical_tag> +{ + typedef strategy::intersection::spherical_segments<> type; +}; + +}} // within::services + +namespace covered_by { namespace services +{ + +template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> +struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, linear_tag, spherical_tag, spherical_tag> +{ + typedef strategy::intersection::spherical_segments<> type; +}; + +template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> +struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, polygonal_tag, spherical_tag, spherical_tag> +{ + typedef strategy::intersection::spherical_segments<> type; +}; + +template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> +struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, linear_tag, spherical_tag, spherical_tag> +{ + typedef strategy::intersection::spherical_segments<> type; +}; + +template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> +struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, polygonal_tag, spherical_tag, spherical_tag> +{ + typedef strategy::intersection::spherical_segments<> type; +}; + +}} // within::services + +} // strategy + + }} // namespace boost::geometry |