summaryrefslogtreecommitdiff
path: root/boost/geometry/strategies/geographic/intersection.hpp
diff options
context:
space:
mode:
Diffstat (limited to 'boost/geometry/strategies/geographic/intersection.hpp')
-rw-r--r--boost/geometry/strategies/geographic/intersection.hpp148
1 files changed, 115 insertions, 33 deletions
diff --git a/boost/geometry/strategies/geographic/intersection.hpp b/boost/geometry/strategies/geographic/intersection.hpp
index 1708c274c0..59a40f281d 100644
--- a/boost/geometry/strategies/geographic/intersection.hpp
+++ b/boost/geometry/strategies/geographic/intersection.hpp
@@ -34,6 +34,7 @@
#include <boost/geometry/strategies/geographic/area.hpp>
#include <boost/geometry/strategies/geographic/distance.hpp>
+#include <boost/geometry/strategies/geographic/envelope_segment.hpp>
#include <boost/geometry/strategies/geographic/parameters.hpp>
#include <boost/geometry/strategies/geographic/side.hpp>
#include <boost/geometry/strategies/intersection.hpp>
@@ -135,6 +136,14 @@ struct geographic_segments
return strategy_type(m_spheroid);
}
+ typedef envelope::geographic_segment<FormulaPolicy, Spheroid, CalculationType>
+ envelope_strategy_type;
+
+ inline envelope_strategy_type get_envelope_strategy() const
+ {
+ return envelope_strategy_type(m_spheroid);
+ }
+
enum intersection_point_flag { ipi_inters = 0, ipi_at_a1, ipi_at_a2, ipi_at_b1, ipi_at_b2 };
template <typename CoordinateType, typename SegmentRatio>
@@ -280,6 +289,8 @@ private:
typedef typename select_calculation_type
<Segment1, Segment2, CalculationType>::type calc_t;
+ static const calc_t c0 = 0;
+
// normalized spheroid
srs::spheroid<calc_t> spheroid = normalized_spheroid<calc_t>(m_spheroid);
@@ -316,31 +327,80 @@ private:
// TODO: no need to call inverse formula if we know that the points are equal
// distance can be set to 0 in this case and azimuth may be not calculated
- bool const is_equal_a1_b1 = equals_point_point(a1, b1);
- bool const is_equal_a2_b1 = equals_point_point(a2, b1);
+ bool is_equal_a1_b1 = equals_point_point(a1, b1);
+ bool is_equal_a2_b1 = equals_point_point(a2, b1);
+ bool degen_neq_coords = false;
- inverse_result res_b1_b2 = inverse_dist_azi::apply(b1_lon, b1_lat, b2_lon, b2_lat, spheroid);
- inverse_result res_b1_a1 = inverse_dist_azi::apply(b1_lon, b1_lat, a1_lon, a1_lat, spheroid);
- inverse_result res_b1_a2 = inverse_dist_azi::apply(b1_lon, b1_lat, a2_lon, a2_lat, spheroid);
- sides.set<0>(is_equal_a1_b1 ? 0 : formula::azimuth_side_value(res_b1_a1.azimuth, res_b1_b2.azimuth),
- is_equal_a2_b1 ? 0 : formula::azimuth_side_value(res_b1_a2.azimuth, res_b1_b2.azimuth));
- if (sides.same<0>())
+ inverse_result res_b1_b2, res_b1_a1, res_b1_a2;
+ if (! b_is_point)
{
- // Both points are at the same side of other segment, we can leave
- return Policy::disjoint();
+ res_b1_b2 = inverse_dist_azi::apply(b1_lon, b1_lat, b2_lon, b2_lat, spheroid);
+ if (math::equals(res_b1_b2.distance, c0))
+ {
+ b_is_point = true;
+ degen_neq_coords = true;
+ }
+ else
+ {
+ res_b1_a1 = inverse_dist_azi::apply(b1_lon, b1_lat, a1_lon, a1_lat, spheroid);
+ if (math::equals(res_b1_a1.distance, c0))
+ {
+ is_equal_a1_b1 = true;
+ }
+ res_b1_a2 = inverse_dist_azi::apply(b1_lon, b1_lat, a2_lon, a2_lat, spheroid);
+ if (math::equals(res_b1_a2.distance, c0))
+ {
+ is_equal_a2_b1 = true;
+ }
+ sides.set<0>(is_equal_a1_b1 ? 0 : formula::azimuth_side_value(res_b1_a1.azimuth, res_b1_b2.azimuth),
+ is_equal_a2_b1 ? 0 : formula::azimuth_side_value(res_b1_a2.azimuth, res_b1_b2.azimuth));
+ if (sides.same<0>())
+ {
+ // Both points are at the same side of other segment, we can leave
+ return Policy::disjoint();
+ }
+ }
}
- bool const is_equal_a1_b2 = equals_point_point(a1, b2);
+ bool is_equal_a1_b2 = equals_point_point(a1, b2);
- inverse_result res_a1_a2 = inverse_dist_azi::apply(a1_lon, a1_lat, a2_lon, a2_lat, spheroid);
- inverse_result res_a1_b1 = inverse_dist_azi::apply(a1_lon, a1_lat, b1_lon, b1_lat, spheroid);
- inverse_result res_a1_b2 = inverse_dist_azi::apply(a1_lon, a1_lat, b2_lon, b2_lat, spheroid);
- sides.set<1>(is_equal_a1_b1 ? 0 : formula::azimuth_side_value(res_a1_b1.azimuth, res_a1_a2.azimuth),
- is_equal_a1_b2 ? 0 : formula::azimuth_side_value(res_a1_b2.azimuth, res_a1_a2.azimuth));
- if (sides.same<1>())
+ inverse_result res_a1_a2, res_a1_b1, res_a1_b2;
+ if (! a_is_point)
{
- // Both points are at the same side of other segment, we can leave
- return Policy::disjoint();
+ res_a1_a2 = inverse_dist_azi::apply(a1_lon, a1_lat, a2_lon, a2_lat, spheroid);
+ if (math::equals(res_a1_a2.distance, c0))
+ {
+ a_is_point = true;
+ degen_neq_coords = true;
+ }
+ else
+ {
+ res_a1_b1 = inverse_dist_azi::apply(a1_lon, a1_lat, b1_lon, b1_lat, spheroid);
+ if (math::equals(res_a1_b1.distance, c0))
+ {
+ is_equal_a1_b1 = true;
+ }
+ res_a1_b2 = inverse_dist_azi::apply(a1_lon, a1_lat, b2_lon, b2_lat, spheroid);
+ if (math::equals(res_a1_b2.distance, c0))
+ {
+ is_equal_a1_b2 = true;
+ }
+ sides.set<1>(is_equal_a1_b1 ? 0 : formula::azimuth_side_value(res_a1_b1.azimuth, res_a1_a2.azimuth),
+ is_equal_a1_b2 ? 0 : formula::azimuth_side_value(res_a1_b2.azimuth, res_a1_a2.azimuth));
+ if (sides.same<1>())
+ {
+ // Both points are at the same side of other segment, we can leave
+ return Policy::disjoint();
+ }
+ }
+ }
+
+ if(a_is_point && b_is_point)
+ {
+ return is_equal_a1_b2
+ ? Policy::degenerate(a, true)
+ : Policy::disjoint()
+ ;
}
// NOTE: at this point the segments may still be disjoint
@@ -370,11 +430,11 @@ private:
{
if (a_is_point)
{
- return collinear_one_degenerated<Policy, calc_t>(a, true, b1, b2, a1, a2, res_b1_b2, res_b1_a1, is_b_reversed);
+ return collinear_one_degenerated<Policy, calc_t>(a, true, b1, b2, a1, a2, res_b1_b2, res_b1_a1, res_b1_a2, is_b_reversed, degen_neq_coords);
}
else if (b_is_point)
{
- return collinear_one_degenerated<Policy, calc_t>(b, false, a1, a2, b1, b2, res_a1_a2, res_a1_b1, is_a_reversed);
+ return collinear_one_degenerated<Policy, calc_t>(b, false, a1, a2, b1, b2, res_a1_a2, res_a1_b1, res_a1_b2, is_a_reversed, degen_neq_coords);
}
else
{
@@ -383,16 +443,16 @@ private:
// use shorter segment
if (res_a1_a2.distance <= res_b1_b2.distance)
{
- calculate_collinear_data(a1, a2, b1, b2, res_a1_a2, res_a1_b1, dist_a1_a2, dist_a1_b1);
- calculate_collinear_data(a1, a2, b1, b2, res_a1_a2, res_a1_b2, dist_a1_a2, dist_a1_b2);
+ calculate_collinear_data(a1, a2, b1, b2, res_a1_a2, res_a1_b1, res_a1_b2, dist_a1_a2, dist_a1_b1);
+ calculate_collinear_data(a1, a2, b1, b2, res_a1_a2, res_a1_b2, res_a1_b1, 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, res_b1_b2, res_b1_a1, dist_b1_b2, dist_b1_a1);
- calculate_collinear_data(b1, b2, a1, a2, res_b1_b2, res_b1_a2, dist_b1_b2, dist_b1_a2);
+ calculate_collinear_data(b1, b2, a1, a2, res_b1_b2, res_b1_a1, res_b1_a2, dist_b1_b2, dist_b1_a1);
+ calculate_collinear_data(b1, b2, a1, a2, res_b1_b2, res_b1_a2, res_b1_a1, 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;
@@ -540,11 +600,13 @@ private:
Point1 const& a1, Point1 const& a2,
Point2 const& b1, Point2 const& b2,
ResultInverse const& res_a1_a2,
- ResultInverse const& res_a1_bi,
- bool is_other_reversed)
+ ResultInverse const& res_a1_b1,
+ ResultInverse const& res_a1_b2,
+ bool is_other_reversed,
+ bool degen_neq_coords)
{
CalcT dist_1_2, dist_1_o;
- if (! calculate_collinear_data(a1, a2, b1, b2, res_a1_a2, res_a1_bi, dist_1_2, dist_1_o))
+ if (! calculate_collinear_data(a1, a2, b1, b2, res_a1_a2, res_a1_b1, res_a1_b2, dist_1_2, dist_1_o, degen_neq_coords))
{
return Policy::disjoint();
}
@@ -565,13 +627,16 @@ private:
static inline bool calculate_collinear_data(Point1 const& a1, Point1 const& a2, // in
Point2 const& b1, Point2 const& b2, // in
ResultInverse const& res_a1_a2, // in
- ResultInverse const& res_a1_bi, // in
- CalcT& dist_a1_a2, CalcT& dist_a1_bi) // out
+ ResultInverse const& res_a1_b1, // in
+ ResultInverse const& res_a1_b2, // in
+ CalcT& dist_a1_a2, // out
+ CalcT& dist_a1_bi, // out
+ bool degen_neq_coords = false) // in
{
dist_a1_a2 = res_a1_a2.distance;
- dist_a1_bi = res_a1_bi.distance;
- if (! same_direction(res_a1_bi.azimuth, res_a1_a2.azimuth))
+ dist_a1_bi = res_a1_b1.distance;
+ if (! same_direction(res_a1_b1.azimuth, res_a1_a2.azimuth))
{
dist_a1_bi = -dist_a1_bi;
}
@@ -589,6 +654,22 @@ private:
return true;
}
+ // check the other endpoint of a very short segment near the pole
+ if (degen_neq_coords)
+ {
+ static CalcT const c0 = 0;
+ if (math::equals(res_a1_b2.distance, c0))
+ {
+ dist_a1_bi = 0;
+ return true;
+ }
+ else if (math::equals(dist_a1_a2 - res_a1_b2.distance, c0))
+ {
+ dist_a1_bi = dist_a1_a2;
+ return true;
+ }
+ }
+
// or i1 is on b
return segment_ratio<CalcT>(dist_a1_bi, dist_a1_a2).on_segment();
}
@@ -816,8 +897,9 @@ private:
static inline bool is_endpoint_equal(CalcT const& dist,
P1 const& ai, P2 const& b1, P2 const& b2)
{
+ static CalcT const c0 = 0;
using geometry::detail::equals::equals_point_point;
- return is_near(dist) && (equals_point_point(ai, b1) || equals_point_point(ai, b2));
+ return is_near(dist) && (equals_point_point(ai, b1) || equals_point_point(ai, b2) || math::equals(dist, c0));
}
template <typename CalcT>