boost/geometry
2
0
Fork 0
mirror of https://github.com/boostorg/geometry.git synced 2026-02-10 23:42:12 +00:00
Code Issues Projects Releases Wiki Activity

[geometry] fix of several robustness issues in cart_intersect and get_turn_info found by testing buffer algorithm. Also restructured cart_intersect such that all robustness issues are handled in separate methods (could be policy later). Finally fixed ever circling iterator with range (for assignment)

Browse Source 
[SVN r77988]
This commit is contained in:
Barend Gehrels
2012-04-15 11:44:15 +00:00
parent f139ba3a6c
commit 2aed65855f
4 changed files with 432 additions and 182 deletions
Download Patch File Download Diff File Expand all files Collapse all files

73
include/boost/geometry/algorithms/detail/overlay/get_turn_info.hpp
View File

@@ -571,6 +571,9 @@ struct collinear : public base_turn_handler
: side_q
;
int const side_pk = SideStrategy::apply(qi, qj, pk);
int const side_qk = SideStrategy::apply(pi, pj, qk);
// See comments above,
// resulting in a strange sort of mathematic rule here:
// The arrival-info multiplied by the relevant side
@@ -578,53 +581,55 @@ struct collinear : public base_turn_handler
int const product = arrival * side_p_or_q;
if(product == 0)
// Robustness: side_p is supposed to be equal to side_pk (because p/q are collinear)
// and side_q to side_qk
bool const robustness_issue = side_pk != side_p || side_qk != side_q;
if (robustness_issue)
{
handle_robustness(ti, arrival, side_p, side_q, side_pk, side_qk);
}
else if(product == 0)
{
both(ti, operation_continue);
}
else
{
int const side_pk = SideStrategy::apply(qi, qj, pk);
int const side_qk = SideStrategy::apply(pi, pj, qk);
if (side_pk != side_p || side_qk != side_q)
{
//std::cout << "ROBUSTNESS -> Collinear "
// << " arr: " << arrival
// << " prod: " << product
// << " dir: " << side_p << " " << side_q
// << " rev: " << side_pk << " " << side_qk
// << std::endl;
handle_robustness(ti, arrival, product,
side_p, side_q, side_pk, side_qk);
}
else
{
// The normal case
ui_else_iu(product == 1, ti);
}
ui_else_iu(product == 1, ti);
}
}
static inline void handle_robustness(TurnInfo& ti,
int arrival, int product,
int side_p, int side_q,
int side_pk, int side_qk)
static inline void handle_robustness(TurnInfo& ti, int arrival,
int side_p, int side_q, int side_pk, int side_qk)
{
bool take_ui = product == 1;
if (product == arrival)
// We take the longer one, i.e. if q arrives in p (arrival == -1),
// then p exceeds q and we should take p for a union...
bool use_p_for_union = arrival == -1;
// ... unless one of the sides consistently directs to the other side
int const consistent_side_p = side_p == side_pk ? side_p : 0;
int const consistent_side_q = side_q == side_qk ? side_q : 0;
if (arrival == -1 && (consistent_side_p == -1 || consistent_side_q == 1))
{
if ((product == 1 && side_p == 1 && side_pk != 1)
|| (product == -1 && side_q == 1 && side_qk != 1))
{
//std::cout << "ROBUSTNESS: -> Reverse" << std::endl;
take_ui = ! take_ui;
}
use_p_for_union = false;
}
if (arrival == 1 && (consistent_side_p == 1 || consistent_side_q == -1))
{
use_p_for_union = true;
}
ui_else_iu(take_ui, ti);
//std::cout << "ROBUSTNESS -> Collinear "
// << " arr: " << arrival
// << " dir: " << side_p << " " << side_q
// << " rev: " << side_pk << " " << side_qk
// << " cst: " << cside_p << " " << cside_q
// << std::boolalpha << " " << use_p_for_union
// << std::endl;
ui_else_iu(use_p_for_union, ti);
}
};
template

134
include/boost/geometry/iterators/ever_circling_iterator.hpp
View File

@@ -95,67 +95,115 @@ private:
bool m_skip_first;
};
template <typename Range>
class ever_circling_range_iterator
: public boost::iterator_adaptor
<
ever_circling_range_iterator<Range>,
typename boost::range_iterator<Range>::type
>
struct ever_circling_range_iterator
: public boost::iterator_facade
<
ever_circling_range_iterator<Range>,
typename boost::range_value<Range>::type const,
boost::random_access_traversal_tag
>
{
public :
typedef typename boost::range_iterator<Range>::type iterator_type;
/// Constructor including the range it is based on
explicit inline ever_circling_range_iterator(Range& range)
: m_range(&range)
, m_iterator(boost::begin(range))
, m_size(boost::size(range))
, m_index(0)
{}
explicit inline ever_circling_range_iterator(Range& range,
bool skip_first = false)
: m_range(range)
, m_skip_first(skip_first)
/// Default constructor
explicit inline ever_circling_range_iterator()
: m_range(NULL)
, m_size(0)
, m_index(0)
{}
inline ever_circling_range_iterator<Range>& operator=(ever_circling_range_iterator<Range> const& source)
{
this->base_reference() = boost::begin(m_range);
m_range = source.m_range;
m_iterator = source.m_iterator;
m_size = source.m_size;
m_index = source.m_index;
return *this;
}
explicit inline ever_circling_range_iterator(Range& range, iterator_type start,
bool skip_first = false)
: m_range(range)
, m_skip_first(skip_first)
{
this->base_reference() = start;
}
/// Navigate to a certain position, should be in [start .. end], if at end
/// it will circle again.
inline void moveto(iterator_type it)
{
this->base_reference() = it;
check_end();
}
typedef std::ptrdiff_t difference_type;
private:
friend class boost::iterator_core_access;
inline void increment(bool possibly_skip = true)
inline typename boost::range_value<Range>::type const& dereference() const
{
(this->base_reference())++;
check_end(possibly_skip);
return *m_iterator;
}
inline void check_end(bool possibly_skip = true)
inline difference_type distance_to(ever_circling_range_iterator<Range> const& other) const
{
if (this->base_reference() == boost::end(m_range))
return other.m_index - this->m_index;
}
inline bool equal(ever_circling_range_iterator<Range> const& other) const
{
return this->m_range == other.m_range
&& this->m_index == other.m_index;
}
inline void increment()
{
++m_index;
if (m_index >= 0 && m_index < m_size)
{
this->base_reference() = boost::begin(m_range);
if (m_skip_first && possibly_skip)
{
increment(false);
}
++m_iterator;
}
else
{
update_iterator();
}
}
Range& m_range;
bool m_skip_first;
inline void decrement()
{
--m_index;
if (m_index >= 0 && m_index < m_size)
{
--m_iterator;
}
else
{
update_iterator();
}
}
inline void advance(difference_type n)
{
if (m_index >= 0 && m_index < m_size
&& m_index + n >= 0 && m_index + n < m_size)
{
m_index += n;
m_iterator += n;
}
else
{
m_index += n;
update_iterator();
}
}
inline void update_iterator()
{
while (m_index < 0)
{
m_index += m_size;
}
m_index = m_index % m_size;
this->m_iterator = boost::begin(*m_range) + m_index;
}
Range* m_range;
typename boost::range_iterator<Range>::type m_iterator;
difference_type m_size;
difference_type m_index;
};

375
include/boost/geometry/strategies/cartesian/cart_intersect.hpp
View File

@@ -120,6 +120,8 @@ struct relate_cartesian_segments
typedef side::side_by_triangle<coordinate_type> side;
side_info sides;
bool collinear_use_first = math::abs(dx_a) + math::abs(dx_b) >= math::abs(dy_a) + math::abs(dy_b);
sides.set<0>
(
side::apply(detail::get_from_index<0>(b)
@@ -141,48 +143,16 @@ struct relate_cartesian_segments
bool collinear = sides.collinear();
if ((sides.zero<0>() && ! sides.zero<1>()) || (sides.zero<1>() && ! sides.zero<0>()))
{
// If one of the segments is collinear, the other must be as well.
// So handle it as collinear.
// (In float/double epsilon margins it can easily occur that one or two of them are -1/1)
// sides.debug();
sides.set<0>(0,0);
sides.set<1>(0,0);
collinear = true;
}
if (sides.meeting())
{
// If two segments meet each other at their segment-points, two sides are zero,
// the other two are not (unless collinear but we don't mean those here).
// However, in near-epsilon ranges it can happen that two sides are zero
// but they do not meet at their segment-points.
// In that case they are nearly collinear and handled as such.
if (! point_equals
(
select(sides.zero_index<0>(), a),
select(sides.zero_index<1>(), b)
)
)
{
//std::cout << "ROBUSTNESS: Suspicious ";
//sides.debug();
//std::cout << std::endl;
//std::cout << " " << d << " " << da << std::endl;
//std::cout << " -> " << geometry::wkt(a) << " " << geometry::wkt(b) << std::endl;
//std::cout << " -> " << dx_a << " " << dy_a << " " << dx_b << " " << dy_b << std::endl;
sides.set<0>(0,0);
sides.set<1>(0,0);
collinear = true;
}
}
robustness_verify_collinear(a, b, sides, collinear);
robustness_verify_meeting(a, b, sides, collinear, collinear_use_first);
if (sides.same<0>() || sides.same<1>())
{
// Both points are at same side of other segment, we can leave
return Policy::disjoint();
if (robustness_verify_same_side(a, b, sides))
{
return Policy::disjoint();
}
}
// Degenerate cases: segments of single point, lying on other segment, non disjoint
@@ -225,62 +195,31 @@ struct relate_cartesian_segments
{
r = da / d;
// Ratio should lie between 0 and 1
// Also these three conditions might be of FP imprecision, the segments were actually (nearly) collinear
promoted_type const zero = 0;
promoted_type const one = 1;
promoted_type const epsilon = std::numeric_limits<double>::epsilon();
if (r < zero || r > one)
if (! robustness_verify_r(a, b, r))
{
if (double_check_disjoint<0>(a, b)
|| double_check_disjoint<1>(a, b))
{
// Can still be disjoint (even if not one is left or right from another)
// This is e.g. in case #snake4 of buffer test.
return Policy::disjoint();
}
//std::cout << "ROBUSTNESS: correction of r " << r << std::endl;
// sides.debug()
// ROBUSTNESS: the r value can in epsilon-cases much larger than 1, while (with perfect arithmetic)
// it should be one. It can be 1.14 or even 1.98049 or 2 (while still intersecting)
// If segments are crossing (we can see that with the sides)
// and one is inside the other, there must be an intersection point.
// We correct for that.
// This is (only) in case #ggl_list_20110820_christophe in unit tests
// If segments are touching (two sides zero), of course they should intersect
// This is (only) in case #buffer_rt_i in the unit tests)
// If one touches in the middle, they also should intersect (#buffer_rt_j)
// Note that even for ttmath r is occasionally > 1, e.g. 1.0000000000000000000000036191231203575
if (r > one)
{
r = one;
}
else if (r < zero)
{
r = zero;
}
return Policy::disjoint();
}
robustness_handle_meeting(a, b, sides, dx_a, dy_a, wx, wy, d, r);
if (robustness_verify_disjoint_at_one_collinear(a, b, sides))
{
return Policy::disjoint();
}
}
}
if(collinear)
{
if (math::abs(dx_a) + math::abs(dx_b) < math::abs(dy_a) + math::abs(dy_b))
if (collinear_use_first)
{
// Y direction contains larger segments (maybe dx is zero)
return relate_collinear<1>(a, b);
return relate_collinear<0>(a, b);
}
else
{
return relate_collinear<0>(a, b);
// Y direction contains larger segments (maybe dx is zero)
return relate_collinear<1>(a, b);
}
}
@@ -291,8 +230,210 @@ struct relate_cartesian_segments
private :
// Ratio should lie between 0 and 1
// Also these three conditions might be of FP imprecision, the segments were actually (nearly) collinear
template <typename T>
static inline bool robustness_verify_r(
segment_type1 const& a, segment_type2 const& b,
T& r)
{
T const zero = 0;
T const one = 1;
if (r < zero || r > one)
{
if (verify_disjoint<0>(a, b) || verify_disjoint<1>(a, b))
{
// Can still be disjoint (even if not one is left or right from another)
// This is e.g. in case #snake4 of buffer test.
return false;
}
//std::cout << "ROBUSTNESS: correction of r " << r << std::endl;
// sides.debug();
// ROBUSTNESS: the r value can in epsilon-cases much larger than 1, while (with perfect arithmetic)
// it should be one. It can be 1.14 or even 1.98049 or 2 (while still intersecting)
// If segments are crossing (we can see that with the sides)
// and one is inside the other, there must be an intersection point.
// We correct for that.
// This is (only) in case #ggl_list_20110820_christophe in unit tests
// If segments are touching (two sides zero), of course they should intersect
// This is (only) in case #buffer_rt_i in the unit tests)
// If one touches in the middle, they also should intersect (#buffer_rt_j)
// Note that even for ttmath r is occasionally > 1, e.g. 1.0000000000000000000000036191231203575
if (r > one)
{
r = one;
}
else if (r < zero)
{
r = zero;
}
}
return true;
}
static inline void robustness_verify_collinear(
segment_type1 const& a, segment_type2 const& b,
side_info& sides,
bool& collinear)
{
if ((sides.zero<0>() && ! sides.zero<1>()) || (sides.zero<1>() && ! sides.zero<0>()))
{
// If one of the segments is collinear, the other must be as well.
// So handle it as collinear.
// (In float/double epsilon margins it can easily occur that one or two of them are -1/1)
// sides.debug();
sides.set<0>(0,0);
sides.set<1>(0,0);
collinear = true;
}
}
static inline void robustness_verify_meeting(
segment_type1 const& a, segment_type2 const& b,
side_info& sides,
bool& collinear, bool& collinear_use_first)
{
if (sides.meeting())
{
// If two segments meet each other at their segment-points, two sides are zero,
// the other two are not (unless collinear but we don't mean those here).
// However, in near-epsilon ranges it can happen that two sides are zero
// but they do not meet at their segment-points.
// In that case they are nearly collinear and handled as such.
if (! point_equals
(
select(sides.zero_index<0>(), a),
select(sides.zero_index<1>(), b)
)
)
{
sides.set<0>(0,0);
sides.set<1>(0,0);
collinear = true;
if (collinear_use_first && analyse_equal<0>(a, b))
{
collinear_use_first = false;
}
else if (! collinear_use_first && analyse_equal<1>(a, b))
{
collinear_use_first = true;
}
}
}
}
// Verifies and if necessary correct missed touch because of robustness
// This is the case at multi_polygon_buffer unittest #rt_m
static inline bool robustness_verify_same_side(
segment_type1 const& a, segment_type2 const& b,
side_info& sides)
{
int corrected = 0;
if (sides.one_touching<0>())
{
if (point_equals(
select(sides.zero_index<0>(), a),
select(0, b)
))
{
sides.correct_to_zero<1, 0>();
corrected = 1;
}
if (point_equals
(
select(sides.zero_index<0>(), a),
select(1, b)
))
{
sides.correct_to_zero<1, 1>();
corrected = 2;
}
}
else if (sides.one_touching<1>())
{
if (point_equals(
select(sides.zero_index<1>(), b),
select(0, a)
))
{
sides.correct_to_zero<0, 0>();
corrected = 3;
}
if (point_equals
(
select(sides.zero_index<1>(), b),
select(1, a)
))
{
sides.correct_to_zero<0, 1>();
corrected = 4;
}
}
return corrected == 0;
}
static inline bool robustness_verify_disjoint_at_one_collinear(
segment_type1 const& a, segment_type2 const& b,
side_info const& sides)
{
if (sides.one_of_all_zero())
{
if (verify_disjoint<0>(a, b) || verify_disjoint<1>(a, b))
{
return true;
}
}
return false;
}
// If r is one, or zero, segments should meet and their endpoints.
// Robustness issue: check if this is really the case.
// It turns out to be no problem, see buffer test #rt_s1 (and there are many cases generated)
// It generates an "ends in the middle" situation which is correct.
template <typename T, typename R>
static inline void robustness_handle_meeting(segment_type1 const& a, segment_type2 const& b,
side_info& sides,
T const& dx_a, T const& dy_a, T const& wx, T const& wy,
T const& d, R const& r)
{
return;
T const db = geometry::detail::determinant<T>(dx_a, dy_a, wx, wy);
R const zero = 0;
R const one = 1;
if (math::equals(r, zero) || math::equals(r, one))
{
R rb = db / d;
if (rb <= 0 || rb >= 1 || math::equals(rb, 0) || math::equals(rb, 1))
{
if (sides.one_zero<0>() && ! sides.one_zero<1>()) // or vice versa
{
#if defined(BOOST_GEOMETRY_COUNT_INTERSECTION_EQUAL)
extern int g_count_intersection_equal;
g_count_intersection_equal++;
#endif
sides.debug();
std::cout << "E r=" << r << " r.b=" << rb << " ";
}
}
}
}
template <std::size_t Dimension>
static inline bool double_check_disjoint(segment_type1 const& a,
static inline bool verify_disjoint(segment_type1 const& a,
segment_type2 const& b)
{
coordinate_type a_1, a_2, b_1, b_2;
@@ -321,7 +462,30 @@ private :
;
}
// We cannot use geometry::equals here. Besides that this will be changed
// to compare segment-coordinate-values directly (not necessary to retrieve point first)
template <typename Point1, typename Point2>
static inline bool point_equality(Point1 const& point1, Point2 const& point2,
bool& equals_0, bool& equals_1)
{
equals_0 = math::equals(get<0>(point1), get<0>(point2));
equals_1 = math::equals(get<1>(point1), get<1>(point2));
return equals_0 && equals_1;
}
template <std::size_t Dimension>
static inline bool analyse_equal(segment_type1 const& a, segment_type2 const& b)
{
coordinate_type const a_1 = geometry::get<0, Dimension>(a);
coordinate_type const a_2 = geometry::get<1, Dimension>(a);
coordinate_type const b_1 = geometry::get<0, Dimension>(b);
coordinate_type const b_2 = geometry::get<1, Dimension>(b);
return math::equals(a_1, b_1)
|| math::equals(a_2, b_1)
|| math::equals(a_1, b_2)
|| math::equals(a_2, b_2)
;
}
template <std::size_t Dimension>
static inline return_type relate_collinear(segment_type1 const& a,
@@ -367,25 +531,28 @@ private :
// Handle "equal", in polygon neighbourhood comparisons a common case
// Check if segments are equal...
bool const a1_eq_b1 = math::equals(get<0, 0>(a), get<0, 0>(b))
&& math::equals(get<0, 1>(a), get<0, 1>(b));
bool const a2_eq_b2 = math::equals(get<1, 0>(a), get<1, 0>(b))
&& math::equals(get<1, 1>(a), get<1, 1>(b));
if (a1_eq_b1 && a2_eq_b2)
bool const opposite = a_swapped ^ b_swapped;
bool const both_swapped = a_swapped && b_swapped;
// Check if segments are equal or opposite equal...
bool const swapped_a1_eq_b1 = math::equals(a_1, b_1);
bool const swapped_a2_eq_b2 = math::equals(a_2, b_2);
if (swapped_a1_eq_b1 && swapped_a2_eq_b2)
{
return Policy::segment_equal(a, false);
return Policy::segment_equal(a, opposite);
}
// ... or opposite equal
bool const a1_eq_b2 = math::equals(get<0, 0>(a), get<1, 0>(b))
&& math::equals(get<0, 1>(a), get<1, 1>(b));
bool const a2_eq_b1 = math::equals(get<1, 0>(a), get<0, 0>(b))
&& math::equals(get<1, 1>(a), get<0, 1>(b));
if (a1_eq_b2 && a2_eq_b1)
{
return Policy::segment_equal(a, true);
}
bool const swapped_a2_eq_b1 = math::equals(a_2, b_1);
bool const swapped_a1_eq_b2 = math::equals(a_1, b_2);
bool const a1_eq_b1 = both_swapped ? swapped_a2_eq_b2 : a_swapped ? swapped_a2_eq_b1 : b_swapped ? swapped_a1_eq_b2 : swapped_a1_eq_b1;
bool const a2_eq_b2 = both_swapped ? swapped_a1_eq_b1 : a_swapped ? swapped_a1_eq_b2 : b_swapped ? swapped_a2_eq_b1 : swapped_a2_eq_b2;
bool const a1_eq_b2 = both_swapped ? swapped_a2_eq_b1 : a_swapped ? swapped_a2_eq_b2 : b_swapped ? swapped_a1_eq_b1 : swapped_a1_eq_b2;
bool const a2_eq_b1 = both_swapped ? swapped_a1_eq_b2 : a_swapped ? swapped_a1_eq_b1 : b_swapped ? swapped_a2_eq_b2 : swapped_a2_eq_b1;
// The rest below will return one or two intersections.
@@ -393,7 +560,7 @@ private :
// For IM it is important to know which relates to which. So this information is given,
// without performance penalties to intersection calculation
bool const has_common_points = a1_eq_b1 || a1_eq_b2 || a2_eq_b1 || a2_eq_b2;
bool const has_common_points = swapped_a1_eq_b1 || swapped_a1_eq_b2 || swapped_a2_eq_b1 || swapped_a2_eq_b2;
// "Touch" -> one intersection point -> one but not two common points
@@ -413,7 +580,7 @@ private :
// #4: a2<----a1 b1--->b2 (no arrival at all)
// Where the arranged forms have two forms:
// a_1-----a_2/b_1-------b_2 or reverse (B left of A)
if (has_common_points && (math::equals(a_2, b_1) || math::equals(b_2, a_1)))
if ((swapped_a2_eq_b1 || swapped_a1_eq_b2) && ! swapped_a1_eq_b1 && ! swapped_a2_eq_b2)
{
if (a2_eq_b1) return Policy::collinear_touch(get<1, 0>(a), get<1, 1>(a), 0, -1);
if (a1_eq_b2) return Policy::collinear_touch(get<0, 0>(a), get<0, 1>(a), -1, 0);
@@ -473,7 +640,6 @@ private :
if (a1_eq_b1) return Policy::collinear_interior_boundary_intersect(a_in_b ? a : b, a_in_b, arrival_a, -arrival_a, false);
}
bool const opposite = a_swapped ^ b_swapped;
// "Inside", a completely within b or b completely within a
@@ -535,7 +701,6 @@ private :
the picture might seem wrong but it (supposed to be) is right.
*/
bool const both_swapped = a_swapped && b_swapped;
if (b_1 < a_2 && a_2 < b_2)
{
// Left column, from bottom to top
@@ -557,7 +722,7 @@ private :
;
}
// Nothing should goes through. If any we have made an error
// Robustness: it can occur here...
// std::cout << "Robustness issue, non-logical behaviour" << std::endl;
return Policy::error("Robustness issue, non-logical behaviour");
}
};

32
include/boost/geometry/strategies/side_info.hpp
View File

@@ -44,6 +44,19 @@ public :
sides[Which].second = second;
}
template <int Which, int Index>
inline void correct_to_zero()
{
if (Index == 0)
{
sides[Which].first = 0;
}
else
{
sides[Which].second = 0;
}
}
template <int Which, int Index>
inline int get() const
{
@@ -81,6 +94,15 @@ public :
&& sides[1].second == 0 && sides[0].second == 0);
}
template <int Which>
inline bool one_touching() const
{
// This is normally a situation which can't occur:
// If one is completely left or right, the other cannot touch
return one_zero<Which>()
&& sides[1 - Which].first * sides[1 - Which].second == 1;
}
inline bool meeting() const
{
// Two of them (in each segment) zero, two not
@@ -100,6 +122,16 @@ public :
|| (sides[Which].first != 0 && sides[Which].second == 0);
}
inline bool one_of_all_zero() const
{
int const sum = std::abs(sides[0].first)
+ std::abs(sides[0].second)
+ std::abs(sides[1].first)
+ std::abs(sides[1].second);
return sum == 3;
}
template <int Which>
inline int zero_index() const
{