[relate] add support for no-IP casesin A/A, fix for_each_disjoint_geometry_if for Multi Geometries.

Fix a bug in for_each_disjoint_geometry_if resulting in assert failure if no turns are generated. This affects L/L, L/A and A/A.
Enable all combinations of Areal geometries.
Add ignoring of empty disjoint geometries in L/A and A/A.
Add tests for disjoint geometries.

include/boost/geometry/algorithms/detail/relate/areal_areal.hpp

@@ -31,29 +31,58 @@ namespace boost { namespace geometry
 #ifndef DOXYGEN_NO_DETAIL
 namespace detail { namespace relate {
     
+// WARNING!
+// TODO: In the worst case calling this Pred in a loop for MultiPolygon/MultiPolygon may take O(NM)
+// Use the rtree in this case!
+
 // may be used to set EI and EB for an Areal geometry for which no turns were generated
-template <typename Result, bool TransposeResult>
+template <typename OtherAreal, typename Result, bool TransposeResult>
 class no_turns_aa_pred
 {
 public:
-    no_turns_aa_pred(Result & res)
+    no_turns_aa_pred(OtherAreal const& other_areal, Result & res)
         : m_result_ptr(boost::addressof(res))
+        , m_other_areal_ptr(boost::addressof(other_areal))
+        , m_flags(0)
     {}
 
     template <typename Areal>
     bool operator()(Areal const& areal)
     {
-        // TODO:
-        // add an assertion for empty/invalid geometries
+        typedef typename geometry::point_type<Areal>::type point_type;
+        point_type pt;
+        bool ok = boost::geometry::point_on_border(pt, areal);
 
-        update<interior, exterior, '2', TransposeResult>(*m_result_ptr);
-        update<boundary, exterior, '1', TransposeResult>(*m_result_ptr);
+        // TODO: for now ignore, later throw an exception?
+        if ( !ok )
+            return true;
+
+        // check if the areal is inside the other_areal
+        int pig = detail::within::point_in_geometry(pt, *m_other_areal_ptr);
+        //BOOST_ASSERT( pig != 0 );
+        
+        // inside
+        if ( pig > 0 )
+        {
+            update<interior, interior, '2', TransposeResult>(*m_result_ptr);
+            update<boundary, interior, '1', TransposeResult>(*m_result_ptr);
+            m_flags |= 1;
+        }
+        // outside
+        else
+        {
+            update<interior, exterior, '2', TransposeResult>(*m_result_ptr);
+            update<boundary, exterior, '1', TransposeResult>(*m_result_ptr);
+            m_flags |= 2;
+        }
                     
-        return false;
+        return m_flags != 3;
     }
 
 private:
     Result * m_result_ptr;
+    const OtherAreal * m_other_areal_ptr;
+    int m_flags;
 };
 
 // The implementation of an algorithm calculating relate() for A/A
@@ -90,32 +119,20 @@ struct areal_areal
         turns::get_turns<Geometry1, Geometry2>::apply(turns, geometry1, geometry2, interrupt_policy);
         if ( result.interrupt )
             return;
-//
-//        boundary_checker<Geometry1> boundary_checker1(geometry1);
-//        no_turns_la_linestring_pred
-//            <
-//                Geometry2,
-//                Result,
-//                boundary_checker<Geometry1>,
-//                TransposeResult
-//            > pred1(geometry2, result, boundary_checker1);
-//        for_each_disjoint_geometry_if<0, Geometry1>::apply(turns.begin(), turns.end(), geometry1, pred1);
-//        if ( result.interrupt )
-//            return;
-//
-//        no_turns_la_areal_pred<Result, !TransposeResult> pred2(result);
-//        for_each_disjoint_geometry_if<1, Geometry2>::apply(turns.begin(), turns.end(), geometry2, pred2);
-//        if ( result.interrupt )
-//            return;
-//        
-//        if ( turns.empty() )
-//            return;
-//
-//        // This is set here because in the case if empty Areal geometry were passed
-//        // those shouldn't be set
-//        set<exterior, interior, '2', TransposeResult>(result);// FFFFFF2Fd
-//        if ( result.interrupt )
-//            return;
+
+        no_turns_aa_pred<Geometry2, Result, false> pred1(geometry2, result);
+        for_each_disjoint_geometry_if<0, Geometry1>::apply(turns.begin(), turns.end(), geometry1, pred1);
+        if ( result.interrupt )
+            return;
+
+        no_turns_aa_pred<Geometry1, Result, true> pred2(geometry1, result);
+        for_each_disjoint_geometry_if<1, Geometry2>::apply(turns.begin(), turns.end(), geometry2, pred2);
+        if ( result.interrupt )
+            return;
+        
+        if ( turns.empty() )
+            return;
+
 //
 //        {
 //            // for different multi or same ring id: x, u, i, c

include/boost/geometry/algorithms/detail/relate/follow_helpers.hpp

@@ -14,6 +14,8 @@
 #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_RELATE_FOLLOW_HELPERS_HPP
 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_RELATE_FOLLOW_HELPERS_HPP
 
+#include <boost/geometry/util/range.hpp>
+
 namespace boost { namespace geometry
 {
 
@@ -55,6 +57,39 @@ struct for_each_disjoint_geometry_if<OpId, Geometry, Tag, true>
     static inline bool apply(TurnIt first, TurnIt last,
                              Geometry const& geometry,
                              Pred & pred)
+    {
+        if ( first != last )
+            return for_turns(first, last, geometry, pred);
+        else
+            return for_empty(geometry, pred);
+    }
+
+    template <typename Pred>
+    static inline bool for_empty(Geometry const& geometry,
+                                 Pred & pred)
+    {
+        typedef typename boost::range_iterator<Geometry const>::type iterator;
+
+        bool found = false;
+
+        // O(N)
+        // check predicate for each contained geometry without generated turn
+        for ( iterator it = boost::begin(geometry) ;
+              it != boost::end(geometry) ; ++it )
+        {
+            found = true;
+            bool cont = pred(*it);
+            if ( !cont )
+                break;
+        }
+        
+        return true;
+    }
+
+    template <typename TurnIt, typename Pred>
+    static inline bool for_turns(TurnIt first, TurnIt last,
+                                 Geometry const& geometry,
+                                 Pred & pred)
     {
         BOOST_ASSERT(first != last);
 
@@ -84,9 +119,8 @@ struct for_each_disjoint_geometry_if<OpId, Geometry, Tag, true>
             if ( *it == false )
             {
                 found = true;
-                bool cont = pred(
-                                *(boost::begin(geometry)
-                                    + std::distance(detected_intersections.begin(), it)));
+                bool cont = pred(range::at(geometry,
+                                           std::distance(detected_intersections.begin(), it)));
                 if ( !cont )
                     break;
             }

include/boost/geometry/algorithms/detail/relate/linear_areal.hpp

@@ -31,7 +31,8 @@ namespace boost { namespace geometry
 #ifndef DOXYGEN_NO_DETAIL
 namespace detail { namespace relate {
 
-// TODO: In the worst case for MultiLinestring/MultiPolygon this is O(NM)
+// WARNING!
+// TODO: In the worst case calling this Pred in a loop for MultiLinestring/MultiPolygon may take O(NM)
 // Use the rtree in this case!
 
 // may be used to set IE and BE for a Linear geometry for which no turns were generated
@@ -118,7 +119,15 @@ public:
     bool operator()(Areal const& areal)
     {
         // TODO:
-        // add an assertion for empty/invalid geometries
+        // handle empty/invalid geometries in a different way than this:
+
+        typedef typename geometry::point_type<Areal>::type point_type;
+        point_type dummy;
+        bool ok = boost::geometry::point_on_border(dummy, areal);
+
+        // TODO: for now ignore, later throw an exception?
+        if ( !ok )
+            return true;
 
         update<interior, exterior, '2', TransposeResult>(*m_result_ptr);
         update<boundary, exterior, '1', TransposeResult>(*m_result_ptr);

include/boost/geometry/algorithms/detail/relate/relate.hpp

@@ -102,6 +102,7 @@ struct relate<Geometry, Point, Tag1, point_tag>
     : detail::relate::geometry_point<Geometry, Point>
 {};
 
+
 template <typename Linestring1, typename Linestring2>
 struct relate<Linestring1, Linestring2, linestring_tag, linestring_tag>
     : detail::relate::linear_linear<Linestring1, Linestring2>
@@ -122,6 +123,7 @@ struct relate<MultiLinestring1, MultiLinestring2, multi_linestring_tag, multi_li
     : detail::relate::linear_linear<MultiLinestring1, MultiLinestring2>
 {};
 
+
 template <typename Linestring, typename Polygon>
 struct relate<Linestring, Polygon, linestring_tag, polygon_tag>
     : detail::relate::linear_areal<Linestring, Polygon>
@@ -178,11 +180,52 @@ struct relate<MultiPolygon, MultiLinestring, multi_polygon_tag, multi_linestring
     : detail::relate::areal_linear<MultiPolygon, MultiLinestring>
 {};
 
+
 template <typename Polygon1, typename Polygon2>
 struct relate<Polygon1, Polygon2, polygon_tag, polygon_tag>
     : detail::relate::areal_areal<Polygon1, Polygon2>
 {};
 
+template <typename Ring, typename Polygon>
+struct relate<Ring, Polygon, ring_tag, polygon_tag>
+    : detail::relate::areal_areal<Ring, Polygon>
+{};
+
+template <typename Polygon, typename Ring>
+struct relate<Polygon, Ring, polygon_tag, ring_tag>
+    : detail::relate::areal_areal<Polygon, Ring>
+{};
+
+template <typename Ring1, typename Ring2>
+struct relate<Ring1, Ring2, ring_tag, ring_tag>
+    : detail::relate::areal_areal<Ring1, Ring2>
+{};
+
+template <typename Polygon, typename MultiPolygon>
+struct relate<Polygon, MultiPolygon, polygon_tag, multi_polygon_tag>
+    : detail::relate::areal_areal<Polygon, MultiPolygon>
+{};
+
+template <typename MultiPolygon, typename Polygon>
+struct relate<MultiPolygon, Polygon, multi_polygon_tag, polygon_tag>
+    : detail::relate::areal_areal<MultiPolygon, Polygon>
+{};
+
+template <typename Ring, typename MultiPolygon>
+struct relate<Ring, MultiPolygon, ring_tag, multi_polygon_tag>
+    : detail::relate::areal_areal<Ring, MultiPolygon>
+{};
+
+template <typename MultiPolygon, typename Ring>
+struct relate<MultiPolygon, Ring, multi_polygon_tag, ring_tag>
+    : detail::relate::areal_areal<MultiPolygon, Ring>
+{};
+
+template <typename MultiPolygon1, typename MultiPolygon2>
+struct relate<MultiPolygon1, MultiPolygon2, multi_polygon_tag, multi_polygon_tag>
+    : detail::relate::areal_areal<MultiPolygon1, MultiPolygon2>
+{};
+
 }} // namespace detail_dispatch::relate
 
 namespace detail { namespace relate {

test/algorithms/relate.cpp

@@ -615,6 +615,11 @@ void test_multi_linestring_multi_polygon()
     test_geometry<mls, mpoly>("MULTILINESTRING((0 0,10 0,10 10),(10 10,0 10,0 0),(20 20,50 50,20 80,20 20))",
                               "MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))",
                               "F11FFF2F2");
+
+    // disjoint
+    test_geometry<mls, mpoly>("MULTILINESTRING((20 20,30 30),(30 30,40 40))",
+                              "MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))",
+                              "FF1FF0212");
 }
 
 template <typename P>
@@ -640,6 +645,9 @@ void polygon_polygon()
     test_geometry<poly, poly>("POLYGON((0 0,0 10,10 10,10 0,0 0))",
                               "POLYGON((5 5,5 10,6 5,5 5))",
                               "212F01FF2");
+    test_geometry<poly, poly>("POLYGON((0 0,0 10,10 10,10 0,0 0))",
+                              "POLYGON((5 5,5 6,6 6,6 5,5 5))",
+                              "212FF1FF2");
 
     // overlapping
     test_geometry<poly, poly>("POLYGON((0 0,0 10,10 10,10 0,0 0))",
@@ -653,6 +661,32 @@ void polygon_polygon()
     test_geometry<poly, poly>("POLYGON((0 0,0 10,10 10,10 0,0 0))",
                               "POLYGON((10 10,10 5,10 0,5 0,0 0,0 10,5 10,10 10))",
                               "2FFF1FFF2");
+
+    // disjoint
+    test_geometry<poly, poly>("POLYGON((0 0,0 10,10 10,10 0,0 0))",
+                              "POLYGON((0 20,0 30,10 30,10 20,0 20))",
+                              "FF2FF1212");
+}
+
+template <typename P>
+void multi_polygon_multi_polygon()
+{
+    typedef bg::model::polygon<P> poly;
+    typedef bg::model::multi_polygon<poly> mpoly;
+
+    test_geometry<mpoly, mpoly>("MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))",
+                                "MULTIPOLYGON(((5 5,5 10,6 10,6 5,5 5)),((0 20,0 30,10 30,10 20,0 20)))",
+                                "212F11212");
+    test_geometry<mpoly, mpoly>("MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)),((0 20,0 30,10 30,10 20,0 20)))",
+                                "MULTIPOLYGON(((5 5,5 10,6 10,6 5,5 5)))",
+                                "212F11FF2");
+
+    test_geometry<mpoly, mpoly>("MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))",
+                                "MULTIPOLYGON(((5 5,5 6,6 6,6 5,5 5)),((0 20,0 30,10 30,10 20,0 20)))",
+                                "212FF1212");
+    test_geometry<mpoly, mpoly>("MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)),((0 20,0 30,10 30,10 20,0 20)))",
+                                "MULTIPOLYGON(((5 5,5 6,6 6,6 5,5 5)))",
+                                "212FF1FF2");
 }
 
 template <typename P>
@@ -669,6 +703,7 @@ void test_all()
     test_linestring_multi_polygon<P>();
     test_multi_linestring_multi_polygon<P>();
     polygon_polygon<P>();
+    multi_polygon_multi_polygon<P>();
 }
 
 int test_main( int , char* [] )