[Boost-commit] svn:boost r77988 - in trunk/boost/geometry: algorithms/detail/overlay iterators strategies strategies/cartesian

Subject: [Boost-commit] svn:boost r77988 - in trunk/boost/geometry: algorithms/detail/overlay iterators strategies strategies/cartesian
From: barend.gehrels_at_[hidden]
Date: 2012-04-15 07:44:16

Author: barendgehrels
Date: 2012-04-15 07:44:15 EDT (Sun, 15 Apr 2012)
New Revision: 77988
URL: http://svn.boost.org/trac/boost/changeset/77988

Log:
[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)
Text files modified: 
   trunk/boost/geometry/algorithms/detail/overlay/get_turn_info.hpp |    73 ++++---                                 
   trunk/boost/geometry/iterators/ever_circling_iterator.hpp        |   126 +++++++++----                           
   trunk/boost/geometry/strategies/cartesian/cart_intersect.hpp     |   375 ++++++++++++++++++++++++++++----------- 
   trunk/boost/geometry/strategies/side_info.hpp                    |    32 +++                                     
   4 files changed, 428 insertions(+), 178 deletions(-)

Modified: trunk/boost/geometry/algorithms/detail/overlay/get_turn_info.hpp
==============================================================================
--- trunk/boost/geometry/algorithms/detail/overlay/get_turn_info.hpp	(original)
+++ trunk/boost/geometry/algorithms/detail/overlay/get_turn_info.hpp	2012-04-15 07:44:15 EDT (Sun, 15 Apr 2012)
@@ -571,6 +571,9 @@
             : 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 @@
 
         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

Modified: trunk/boost/geometry/iterators/ever_circling_iterator.hpp
==============================================================================
--- trunk/boost/geometry/iterators/ever_circling_iterator.hpp	(original)
+++ trunk/boost/geometry/iterators/ever_circling_iterator.hpp	2012-04-15 07:44:15 EDT (Sun, 15 Apr 2012)
@@ -95,67 +95,115 @@
     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)
+    {}
+
+    /// Default constructor
+    explicit inline ever_circling_range_iterator()
+        : m_range(NULL)
+        , m_size(0)
+        , m_index(0)
+    {}
 
-    explicit inline ever_circling_range_iterator(Range& range,
-            bool skip_first = false)
-      : m_range(range)
-      , m_skip_first(skip_first)
+    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)
+    typedef std::ptrdiff_t difference_type;
+
+private:
+    friend class boost::iterator_core_access;
+
+    inline typename boost::range_value<Range>::type const& dereference() const
     {
-        this->base_reference() = start;
+        return *m_iterator;
     }
 
-    /// Navigate to a certain position, should be in [start .. end], if at end
-    /// it will circle again.
-    inline void moveto(iterator_type it)
+    inline difference_type distance_to(ever_circling_range_iterator<Range> const& other) const
     {
-        this->base_reference() = it;
-        check_end();
+        return other.m_index - this->m_index;
     }
 
-private:
+    inline bool equal(ever_circling_range_iterator<Range> const& other) const
+    {
+        return this->m_range == other.m_range
+            && this->m_index == other.m_index;
+    }
 
-    friend class boost::iterator_core_access;
+    inline void increment()
+    {
+        ++m_index;
+        if (m_index >= 0 && m_index < m_size)
+        {
+            ++m_iterator;
+        }
+        else
+        {
+            update_iterator();
+        }
+    }
+
+    inline void decrement()
+    {
+        --m_index;
+        if (m_index >= 0 && m_index < m_size)
+        {
+            --m_iterator;
+        }
+        else
+        {
+            update_iterator();
+        }
+    }
 
-    inline void increment(bool possibly_skip = true)
+    inline void advance(difference_type n)
     {
-        (this->base_reference())++;
-        check_end(possibly_skip);
+        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 check_end(bool possibly_skip = true)
+    inline void update_iterator()
     {
-        if (this->base_reference() == boost::end(m_range))
+        while (m_index < 0)
         {
-            this->base_reference() = boost::begin(m_range);
-            if (m_skip_first && possibly_skip)
-            {
-                increment(false);
-            }
+            m_index += m_size;
         }
+		m_index = m_index % m_size;
+        this->m_iterator = boost::begin(*m_range) + m_index;
     }
 
-    Range& m_range;
-    bool m_skip_first;
+    Range* m_range;
+    typename boost::range_iterator<Range>::type m_iterator;
+    difference_type m_size;
+    difference_type m_index;
 };
 
 

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-04-15 07:44:15 EDT (Sun, 15 Apr 2012)
@@ -120,6 +120,8 @@
         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 @@
 
         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 @@
                         {
                                 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 (! robustness_verify_r(a, b, r))
+                {
+                    return Policy::disjoint();
+                }
+
+                robustness_handle_meeting(a, b, sides, dx_a, dy_a, wx, wy, d, r);
 
-                if (r < zero || r > one)
+                if (robustness_verify_disjoint_at_one_collinear(a, b, sides))
                 {
-					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();
                 }
+
                         }
                 }
 
         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 @@
 
 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 @@
             ;
     }
 
+    // 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,33 +531,36 @@
 
         // 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)
-        {
-            return Policy::segment_equal(a, false);
-        }
-
-        // ... 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)
+        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, true);
+            return Policy::segment_equal(a, opposite);
         }
 
+		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.
         // The delegated class can decide which is the intersection point, or two, build the Intersection Matrix (IM)
         // 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 @@
         // #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 @@
             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 @@
           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 @@
                 ;
         }
         // 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");
     }
 };

Modified: trunk/boost/geometry/strategies/side_info.hpp
==============================================================================
--- trunk/boost/geometry/strategies/side_info.hpp	(original)
+++ trunk/boost/geometry/strategies/side_info.hpp	2012-04-15 07:44:15 EDT (Sun, 15 Apr 2012)
@@ -45,6 +45,19 @@
     }
 
     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
     {
         return Index == 0 ? sides[Which].first : sides[Which].second;
@@ -81,6 +94,15 @@
             && 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 @@
             || (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
     {