From d4418204ec99a8c779bae1dcf81685cb6754ab4b Mon Sep 17 00:00:00 2001
From: Adam Wulkiewicz <adam.wulkiewicz@gmail.com>
Date: Thu, 27 Mar 2014 00:30:43 +0100
Subject: [PATCH] [relate] add support for non-simple non-equal Polygons in
 A/A.

Currently simple interrupt policy is used.
This approach works is there are IPs, Polygons are simple and not equal.
---
 .../algorithms/detail/relate/areal_areal.hpp  | 61 +++++++++++--------
 .../algorithms/detail/relate/linear_areal.hpp |  2 +
 .../algorithms/detail/relate/relate.hpp       |  7 ++-
 test/algorithms/relate.cpp                    | 42 ++++++++-----
 4 files changed, 70 insertions(+), 42 deletions(-)

diff --git a/include/boost/geometry/algorithms/detail/relate/areal_areal.hpp b/include/boost/geometry/algorithms/detail/relate/areal_areal.hpp
index b040151df..065cb0cca 100644
--- a/include/boost/geometry/algorithms/detail/relate/areal_areal.hpp
+++ b/include/boost/geometry/algorithms/detail/relate/areal_areal.hpp
@@ -57,7 +57,7 @@ private:
 };
 
 // The implementation of an algorithm calculating relate() for A/A
-template <typename Geometry1, typename Geometry2, bool TransposeResult = false>
+template <typename Geometry1, typename Geometry2>
 struct areal_areal
 {
     // check Linear / Areal
@@ -72,8 +72,10 @@ struct areal_areal
     template <typename Result>
     static inline void apply(Geometry1 const& geometry1, Geometry2 const& geometry2, Result & result)
     {
+// TODO: If Areal geometry may have infinite size, change the following line:
+
         // The result should be FFFFFFFFF
-        set<exterior, exterior, result_dimension<Geometry2>::value, TransposeResult>(result);// FFFFFFFFd, d in [1,9] or T
+        set<exterior, exterior, result_dimension<Geometry2>::value>(result);// FFFFFFFFd, d in [1,9] or T
 
         if ( result.interrupt )
             return;
@@ -83,9 +85,9 @@ struct areal_areal
         typedef typename std::vector<turn_type>::iterator turn_iterator;
         std::vector<turn_type> turns;
 
-        //interrupt_policy_areal_areal<Geometry1, Geometry2, Result> interrupt_policy(geometry1, geometry2, result);
+        interrupt_policy_areal_areal<Geometry1, Geometry2, Result> interrupt_policy(geometry1, geometry2, result);
 
-        turns::get_turns<Geometry1, Geometry2>::apply(turns, geometry1, geometry2/*, interrupt_policy*/);
+        turns::get_turns<Geometry1, Geometry2>::apply(turns, geometry1, geometry2, interrupt_policy);
         if ( result.interrupt )
             return;
 //
@@ -272,35 +274,40 @@ struct areal_areal
             
             for (iterator it = boost::begin(turns) ; it != boost::end(turns) ; ++it)
             {
-                if ( it->operations[0].operation == overlay::operation_intersection )
-                {
-                    bool const no_interior_rings
-                        = geometry::num_interior_rings(
-                                single_geometry(m_areal, it->operations[1].seg_id)) == 0;
-
-                    // WARNING! THIS IS TRUE ONLY IF THE POLYGON IS SIMPLE!
-                    // OR WITHOUT INTERIOR RINGS (AND OF COURSE VALID)
-                    if ( no_interior_rings )
-                        update<interior, interior, '1', TransposeResult>(m_result);
-                }
-                else if ( it->operations[0].operation == overlay::operation_continue )
-                {
-                    update<interior, boundary, '1', TransposeResult>(m_result);
-                    is_boundary_found = true;
-                }
-                else if ( ( it->operations[0].operation == overlay::operation_union
-                         || it->operations[0].operation == overlay::operation_blocked )
-                       && it->operations[0].position == overlay::position_middle )
-                {
-// TODO: here we could also check the boundaries and set BB at this point
-                    update<interior, boundary, '0', TransposeResult>(m_result);
-                }
+                per_turn<0, false>(*it);
+                per_turn<1, true>(*it);
             }
 
             return m_result.interrupt;
         }
 
     private:
+        template <std::size_t Id, bool TransposeResult, typename Turn>
+        inline void per_turn(Turn const& turn)
+        {
+// THIS WON'T WORK FOR NON-SIMPLE GEOMETRIES!
+
+            overlay::operation_type op = turn.operations[Id].operation;
+
+            if ( op == overlay::operation_intersection )
+            {
+                update<interior, interior, '2', TransposeResult>(m_result);
+                update<boundary, interior, '1', TransposeResult>(m_result);
+                update<boundary, boundary, '0', TransposeResult>(m_result);
+            }
+            else if ( op == overlay::operation_continue ||
+                      op == overlay::operation_blocked )
+            {
+                update<boundary, boundary, '1', TransposeResult>(m_result);
+            }
+            else if ( op == overlay::operation_union )
+            {
+                update<boundary, boundary, '0', TransposeResult>(m_result);
+                update<boundary, exterior, '1', TransposeResult>(m_result);
+                update<interior, exterior, '2', TransposeResult>(m_result);
+            }
+        }
+
         Result & m_result;
         Geometry1 const& m_geometry1;
         Geometry2 const& m_geometry2;
diff --git a/include/boost/geometry/algorithms/detail/relate/linear_areal.hpp b/include/boost/geometry/algorithms/detail/relate/linear_areal.hpp
index cf377b097..f2e41a583 100644
--- a/include/boost/geometry/algorithms/detail/relate/linear_areal.hpp
+++ b/include/boost/geometry/algorithms/detail/relate/linear_areal.hpp
@@ -146,6 +146,8 @@ struct linear_areal
     template <typename Result>
     static inline void apply(Geometry1 const& geometry1, Geometry2 const& geometry2, Result & result)
     {
+// TODO: If Areal geometry may have infinite size, change the following line:
+
         // The result should be FFFFFFFFF
         set<exterior, exterior, result_dimension<Geometry2>::value, TransposeResult>(result);// FFFFFFFFd, d in [1,9] or T
 
diff --git a/include/boost/geometry/algorithms/detail/relate/relate.hpp b/include/boost/geometry/algorithms/detail/relate/relate.hpp
index 3a125a41a..07ea10138 100644
--- a/include/boost/geometry/algorithms/detail/relate/relate.hpp
+++ b/include/boost/geometry/algorithms/detail/relate/relate.hpp
@@ -47,7 +47,7 @@
 #include <boost/geometry/algorithms/detail/relate/point_geometry.hpp>
 #include <boost/geometry/algorithms/detail/relate/linear_linear.hpp>
 #include <boost/geometry/algorithms/detail/relate/linear_areal.hpp>
-//#include <boost/geometry/algorithms/detail/relate/areal_areal.hpp>
+#include <boost/geometry/algorithms/detail/relate/areal_areal.hpp>
 
 namespace boost { namespace geometry
 {
@@ -178,6 +178,11 @@ 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>
+{};
+
 }} // namespace detail_dispatch::relate
 
 namespace detail { namespace relate {
diff --git a/test/algorithms/relate.cpp b/test/algorithms/relate.cpp
index 46d2cc8bb..66a6b25fc 100644
--- a/test/algorithms/relate.cpp
+++ b/test/algorithms/relate.cpp
@@ -622,23 +622,37 @@ void polygon_polygon()
 {
     typedef bg::model::polygon<P> poly;
 
-    //to_svg<poly, poly>("POLYGON((0 0,0 10,10 10,10 0,0 0))",
-    //                   "POLYGON((10 0,10 10,20 10,20 0,10 0))",
-    //                   "pp1.svg");
+    // touching
+    test_geometry<poly, poly>("POLYGON((0 0,0 10,10 10,10 0,0 0))",
+                              "POLYGON((10 0,10 10,20 10,20 0,10 0))",
+                              "FF2F11212");
+    test_geometry<poly, poly>("POLYGON((0 0,0 10,10 10,10 0,0 0))",
+                              "POLYGON((0 -10,0 0,10 0,10 -10,0 -10))",
+                              "FF2F11212");
+    test_geometry<poly, poly>("POLYGON((0 0,0 10,10 10,10 0,0 0))",
+                              "POLYGON((10 0,15 10,20 10,20 0,10 0))",
+                              "FF2F01212");
 
-    //to_svg<poly, poly>("POLYGON((0 0,0 10,10 10,15 5,10 0,0 0))",
-    //                   "POLYGON((10 0,5 5,10 10,20 10,20 0,10 0))",
-    //                   "pp11.svg");
+    // containing
+    test_geometry<poly, poly>("POLYGON((0 0,0 10,10 10,10 0,0 0))",
+                              "POLYGON((5 5,5 10,6 10,6 5,5 5))",
+                              "212F11FF2");
+    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");
 
-    //to_svg<poly, poly>("POLYGON((0 0,0 10,10 10,10 0,0 0))",
-    //                   "POLYGON((5 0,5 10,20 10,20 0,5 0))",
-    //                   "pp2.svg");
+    // overlapping
+    test_geometry<poly, poly>("POLYGON((0 0,0 10,10 10,10 0,0 0))",
+                              "POLYGON((5 0,5 10,20 10,20 0,5 0))",
+                              "212111212");
+    test_geometry<poly, poly>("POLYGON((0 0,0 10,10 10,15 5,10 0,0 0))",
+                              "POLYGON((10 0,5 5,10 10,20 10,20 0,10 0))",
+                              "212101212");
 
-    //to_svg<poly, poly>("POLYGON((0 0,0 10,10 10,10 0,0 0))",
-    //                   "POLYGON((0 -10,0 0,10 0,10 -10,0 -10))",
-    //                   "pp3.svg");
-
-    
+    // equal
+    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");
 }
 
 template <typename P>