Subject: [Boost-commit] svn:boost r77351 - in trunk/boost/geometry/strategies: . cartesian
From: barend.gehrels_at_[hidden]
Date: 2012-03-16 12:58:27
Author: barendgehrels
Date: 2012-03-16 12:58:26 EDT (Fri, 16 Mar 2012)
New Revision: 77351
URL: http://svn.boost.org/trac/boost/changeset/77351
Log:
[geometry] robustness fixes (all found by buffer robustness tests)
Text files modified:
trunk/boost/geometry/strategies/cartesian/cart_intersect.hpp | 129 ++++++++++++++++++++++++++++-----------
trunk/boost/geometry/strategies/side_info.hpp | 20 ++++++
2 files changed, 112 insertions(+), 37 deletions(-)
Modified: trunk/boost/geometry/strategies/cartesian/cart_intersect.hpp
==============================================================================
--- trunk/boost/geometry/strategies/cartesian/cart_intersect.hpp (original)
+++ trunk/boost/geometry/strategies/cartesian/cart_intersect.hpp 2012-03-16 12:58:26 EDT (Fri, 16 Mar 2012)
@@ -152,6 +152,33 @@
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;
+ }
+ }
+
if (sides.same<0>() || sides.same<1>())
{
// Both points are at same side of other segment, we can leave
@@ -206,51 +233,47 @@
if (r < zero || r > one)
{
- if (sides.crossing() || sides.touching())
+ 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)
{
- // ROBUSTNESS: the r value can in epsilon-cases be 1.14, while (with perfect arithmetic)
- // it should be one. 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)
- // TODO: find more cases
- if (r > one)
- {
- // std::cout << "ROBUSTNESS: correction of r " << r << std::endl;
- r = one;
- }
- else if (r < zero)
- {
- // std::cout << "ROBUSTNESS: correction of r " << r << std::endl;
- r = zero;
- }
+ r = one;
+ }
+ else if (r < zero)
+ {
+ r = zero;
}
-
- if (r < zero)
- {
- if (r < -epsilon)
- {
- return Policy::disjoint();
- }
- r = zero;
- }
- else if (r > one)
- {
- if (r > one + epsilon)
- {
- return Policy::disjoint();
- }
- r = one;
- }
}
}
}
if(collinear)
{
- if (math::abs(dx_a) + math::abs(dx_b) < math::abs(dy_a) + math::abs(dy_a))
+ if (math::abs(dx_a) + math::abs(dx_b) < math::abs(dy_a) + math::abs(dy_b))
{
// Y direction contains larger segments (maybe dx is zero)
return relate_collinear<1>(a, b);
@@ -269,6 +292,38 @@
private :
template <std::size_t Dimension>
+ static inline bool double_check_disjoint(segment_type1 const& a,
+ segment_type2 const& b)
+ {
+ coordinate_type a_1, a_2, b_1, b_2;
+ bool a_swapped = false, b_swapped = false;
+ detail::segment_arrange<Dimension>(a, a_1, a_2, a_swapped);
+ detail::segment_arrange<Dimension>(b, b_1, b_2, b_swapped);
+ return math::smaller(a_2, b_1) || math::larger(a_1, b_2);
+ }
+
+ template <typename Segment>
+ static inline typename point_type<Segment>::type select(int index, Segment const& segment)
+ {
+ return index == 0
+ ? detail::get_from_index<0>(segment)
+ : detail::get_from_index<1>(segment)
+ ;
+ }
+
+ // 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_equals(Point1 const& point1, Point2 const& point2)
+ {
+ return math::equals(get<0>(point1), get<0>(point2))
+ && math::equals(get<1>(point1), get<1>(point2))
+ ;
+ }
+
+
+
+ template <std::size_t Dimension>
static inline return_type relate_collinear(segment_type1 const& a,
segment_type2 const& b)
{
Modified: trunk/boost/geometry/strategies/side_info.hpp
==============================================================================
--- trunk/boost/geometry/strategies/side_info.hpp (original)
+++ trunk/boost/geometry/strategies/side_info.hpp 2012-03-16 12:58:26 EDT (Fri, 16 Mar 2012)
@@ -81,12 +81,32 @@
&& sides[1].second == 0 && sides[0].second == 0);
}
+ inline bool meeting() const
+ {
+ // Two of them (in each segment) zero, two not
+ return one_zero<0>() && one_zero<1>();
+ }
+
template <int Which>
inline bool zero() const
{
return sides[Which].first == 0 && sides[Which].second == 0;
}
+ template <int Which>
+ inline bool one_zero() const
+ {
+ return (sides[Which].first == 0 && sides[Which].second != 0)
+ || (sides[Which].first != 0 && sides[Which].second == 0);
+ }
+
+ template <int Which>
+ inline int zero_index() const
+ {
+ return sides[Which].first == 0 ? 0 : 1;
+ }
+
+
inline void debug() const
{
std::cout << sides[0].first << " "