Correct calculations of overlap and contains operations over polygons.

diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c
index 24157a53c59955f70abb2669df8f3fb57707d3da..ea7279e8c7c14d4cc8c0ff177508e15c61d016a4 100644 (file)
--- a/src/backend/utils/adt/geo_ops.c
+++ b/src/backend/utils/adt/geo_ops.c
@@ -8,7 +8,7 @@
  *
  *
  * IDENTIFICATION
- *       $PostgreSQL: pgsql/src/backend/utils/adt/geo_ops.c,v 1.102 2009/06/23 16:25:02 tgl Exp $
+ *       $PostgreSQL: pgsql/src/backend/utils/adt/geo_ops.c,v 1.103 2009/07/28 09:47:59 teodor Exp $
  *
  *-------------------------------------------------------------------------
  */
@@ -66,6 +66,8 @@ static bool has_interpt_sl(LSEG *lseg, LINE *line);
 static double dist_pl_internal(Point *pt, LINE *line);
 static double dist_ps_internal(Point *pt, LSEG *lseg);
 static Point *line_interpt_internal(LINE *l1, LINE *l2);
+static bool lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start);
+static Point* lseg_interpt_internal(LSEG *l1, LSEG *l2);
 
 
 /*
@@ -2352,15 +2354,9 @@ lseg_center(PG_FUNCTION_ARGS)
        PG_RETURN_POINT_P(result);
 }
 
-
-/* lseg_interpt -
- *             Find the intersection point of two segments (if any).
- */
-Datum
-lseg_interpt(PG_FUNCTION_ARGS)
+static Point*
+lseg_interpt_internal(LSEG *l1, LSEG *l2)
 {
-       LSEG       *l1 = PG_GETARG_LSEG_P(0);
-       LSEG       *l2 = PG_GETARG_LSEG_P(1);
        Point      *result;
        LINE            tmp1,
                                tmp2;
@@ -2372,7 +2368,7 @@ lseg_interpt(PG_FUNCTION_ARGS)
        line_construct_pts(&tmp2, &l2->p[0], &l2->p[1]);
        result = line_interpt_internal(&tmp1, &tmp2);
        if (!PointerIsValid(result))
-               PG_RETURN_NULL();
+               return NULL;
 
        /*
         * If the line intersection point isn't within l1 (or equivalently l2),
@@ -2380,7 +2376,10 @@ lseg_interpt(PG_FUNCTION_ARGS)
         */
        if (!on_ps_internal(result, l1) ||
                !on_ps_internal(result, l2))
-               PG_RETURN_NULL();
+       {
+               pfree(result);
+               return NULL;
+       }
 
        /*
         * If there is an intersection, then check explicitly for matching
@@ -2400,6 +2399,23 @@ lseg_interpt(PG_FUNCTION_ARGS)
                result->y = l1->p[1].y;
        }
 
+       return result;
+}
+
+/* lseg_interpt -
+ *             Find the intersection point of two segments (if any).
+ */
+Datum
+lseg_interpt(PG_FUNCTION_ARGS)
+{
+       LSEG       *l1 = PG_GETARG_LSEG_P(0);
+       LSEG       *l2 = PG_GETARG_LSEG_P(1);
+       Point      *result;
+       
+       result = lseg_interpt_internal(l1, l2);
+       if (!PointerIsValid(result))
+               PG_RETURN_NULL();
+
        PG_RETURN_POINT_P(result);
 }
 
@@ -3742,10 +3758,7 @@ poly_same(PG_FUNCTION_ARGS)
 }
 
 /*-----------------------------------------------------------------
- * Determine if polygon A overlaps polygon B by determining if
- * their bounding boxes overlap.
- *
- * XXX ought to do a more correct check!
+ * Determine if polygon A overlaps polygon B 
  *-----------------------------------------------------------------*/
 Datum
 poly_overlap(PG_FUNCTION_ARGS)
@@ -3754,7 +3767,54 @@ poly_overlap(PG_FUNCTION_ARGS)
        POLYGON    *polyb = PG_GETARG_POLYGON_P(1);
        bool            result;
 
-       result = box_ov(&polya->boundbox, &polyb->boundbox);
+       /* Quick check by bounding box */
+       result = (polya->npts > 0 && polyb->npts > 0 && 
+                       box_ov(&polya->boundbox, &polyb->boundbox)) ? true : false;
+
+       /*
+        * Brute-force algorithm - try to find intersected edges,
+        * if so then polygons are overlapped else check is one 
+        * polygon inside other or not by testing single point 
+        * of them.
+        */
+       if (result)
+       {
+               int             ia, ib;
+               LSEG    sa, sb;
+
+               /* Init first of polya's edge with last point */
+               sa.p[0] = polya->p[polya->npts - 1];
+               result = false;
+
+               for(ia=0; ia<polya->npts && result == false; ia++)
+               {
+                       /* Second point of polya's edge is a current one */ 
+                       sa.p[1] = polya->p[ia];
+
+                       /* Init first of polyb's edge with last point */
+                       sb.p[0] = polyb->p[polyb->npts - 1];
+
+                       for(ib=0; ib<polyb->npts && result == false; ib++)
+                       {
+                               sb.p[1] = polyb->p[ib];
+                               result = lseg_intersect_internal(&sa, &sb);
+                               sb.p[0] = sb.p[1];
+                       }
+
+                       /* 
+                        * move current endpoint to the first point
+                        * of next edge
+                        */
+                       sa.p[0] = sa.p[1];
+               }
+
+               if (result==false)
+               {
+                       result = (      point_inside(polya->p, polyb->npts, polyb->p)
+                                               ||
+                                               point_inside(polyb->p, polya->npts, polya->p) );
+               }
+       }
 
        /*
         * Avoid leaking memory for toasted inputs ... needed for rtree indexes
@@ -3765,6 +3825,119 @@ poly_overlap(PG_FUNCTION_ARGS)
        PG_RETURN_BOOL(result);
 }
 
+/*
+ * Tests special kind of segment for in/out of polygon.
+ * Special kind means:
+ *  - point a should be on segment s
+ *  - segment (a,b) should not be contained by s
+ * Returns true if:
+ *  - segment (a,b) is collinear to s and (a,b) is in polygon
+ *  - segment (a,b) s not collinear to s. Note: that doesn't
+ *    mean that segment is in polygon! 
+ */ 
+
+static bool
+touched_lseg_inside_poly(Point *a, Point *b, LSEG *s, POLYGON *poly, int start)
+{
+       /* point a is on s, b is not */
+       LSEG t;
+
+       t.p[0] = *a;
+       t.p[1] = *b;
+       
+#define POINTEQ(pt1, pt2)   (FPeq((pt1)->x, (pt2)->x) && FPeq((pt1)->y, (pt2)->y))
+       if ( POINTEQ(a, s->p) )
+       {
+               if ( on_ps_internal(s->p+1, &t) )
+                       return lseg_inside_poly(b, s->p+1, poly, start);
+       }
+       else if (POINTEQ(a, s->p+1))
+       {
+               if ( on_ps_internal(s->p, &t) )
+                       return lseg_inside_poly(b, s->p, poly, start);
+       }
+       else if ( on_ps_internal(s->p, &t) )
+       {
+               return lseg_inside_poly(b, s->p, poly, start);
+       }
+       else if ( on_ps_internal(s->p+1, &t) )
+       {
+               return lseg_inside_poly(b, s->p+1, poly, start);
+       }
+
+       return true; /* may be not true, but that will check later */
+}
+
+/*
+ * Returns true if segment (a,b) is in polygon, option
+ * start is used for optimization - function checks 
+ * polygon's edges started from start
+ */
+static bool
+lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start)
+{
+       LSEG    s,
+                       t;
+       int     i;
+       bool    res = true,
+                       intersection = false;
+
+       t.p[0] = *a;
+       t.p[1] = *b;
+       s.p[0] = poly->p[( start == 0) ? (poly->npts - 1) : (start - 1)];
+
+       for(i=start; i<poly->npts && res == true; i++)
+       {
+               Point   *interpt;
+
+               s.p[1] = poly->p[i];
+
+               if ( on_ps_internal(t.p, &s) )
+               {
+                       if ( on_ps_internal(t.p+1, &s) )
+                               return true; /* t is contained by s */
+
+                       /* Y-cross */
+                       res = touched_lseg_inside_poly(t.p, t.p+1, &s, poly, i+1);
+               } 
+               else if ( on_ps_internal(t.p+1, &s) )
+               {
+                       /* Y-cross */
+                       res = touched_lseg_inside_poly(t.p+1, t.p, &s, poly, i+1);
+               }
+               else if ( (interpt = lseg_interpt_internal(&t, &s)) != NULL )
+               {
+                       /*
+                        * segments are X-crossing, go to check each subsegment
+                        */
+
+                       intersection = true;
+                       res = lseg_inside_poly(t.p, interpt, poly, i+1);
+                       if (res)
+                               res = lseg_inside_poly(t.p+1, interpt, poly, i+1);
+                       pfree(interpt);
+               }
+
+               s.p[0] = s.p[1];
+       }
+
+       if (res && !intersection)
+       {
+               Point p;
+
+               /*
+                * if X-intersection wasn't found  then check central point
+                * of tested segment. In opposite case we already check all 
+                * subsegments
+                */
+               p.x = (t.p[0].x + t.p[1].x) / 2.0; 
+               p.y = (t.p[0].y + t.p[1].y) / 2.0;
+
+               res = point_inside(&p, poly->npts, poly->p); 
+       }
+
+       return res;
+}
 
 /*-----------------------------------------------------------------
  * Determine if polygon A contains polygon B.
@@ -3775,49 +3948,30 @@ poly_contain(PG_FUNCTION_ARGS)
        POLYGON    *polya = PG_GETARG_POLYGON_P(0);
        POLYGON    *polyb = PG_GETARG_POLYGON_P(1);
        bool            result;
-       int                     i;
 
        /*
         * Quick check to see if bounding box is contained.
         */
-       if (DatumGetBool(DirectFunctionCall2(box_contain,
-                                                                                BoxPGetDatum(&polya->boundbox),
-                                                                                BoxPGetDatum(&polyb->boundbox))))
+       if (polya->npts > 0 && polyb->npts > 0 &&
+                       DatumGetBool(DirectFunctionCall2(box_contain,
+                                                                                       BoxPGetDatum(&polya->boundbox),
+                                                                                       BoxPGetDatum(&polyb->boundbox))))
        {
-               result = true;                  /* assume true for now */
-               for (i = 0; i < polyb->npts; i++)
-               {
-                       if (point_inside(&(polyb->p[i]), polya->npts, &(polya->p[0])) == 0)
-                       {
-#ifdef GEODEBUG
-                               printf("poly_contain- point (%f,%f) not in polygon\n", polyb->p[i].x, polyb->p[i].y);
-#endif
-                               result = false;
-                               break;
-                       }
-               }
-               if (result)
+               int             i;
+               LSEG    s;
+
+               s.p[0] = polyb->p[polyb->npts - 1];
+               result = true;
+
+               for(i=0; i<polyb->npts && result == true; i++)
                {
-                       for (i = 0; i < polya->npts; i++)
-                       {
-                               if (point_inside(&(polya->p[i]), polyb->npts, &(polyb->p[0])) == 1)
-                               {
-#ifdef GEODEBUG
-                                       printf("poly_contain- point (%f,%f) in polygon\n", polya->p[i].x, polya->p[i].y);
-#endif
-                                       result = false;
-                                       break;
-                               }
-                       }
+                       s.p[1] = polyb->p[i];
+                       result = lseg_inside_poly(s.p, s.p+1, polya, 0);
+                       s.p[0] = s.p[1];
                }
        }
        else
        {
-#ifdef GEODEBUG
-               printf("poly_contain- bound box ((%f,%f),(%f,%f)) not inside ((%f,%f),(%f,%f))\n",
-                          polyb->boundbox.low.x, polyb->boundbox.low.y, polyb->boundbox.high.x, polyb->boundbox.high.y,
-                          polya->boundbox.low.x, polya->boundbox.low.y, polya->boundbox.high.x, polya->boundbox.high.y);
-#endif
                result = false;
        }
 

diff --git a/src/test/regress/expected/create_index.out b/src/test/regress/expected/create_index.out
index 7b37cbf40bc3c2f26d1bd6dc8705141768cde989..6e06be5ed534f1c5f1a5ce121a76f8cda01ffbdc 100644 (file)
--- a/src/test/regress/expected/create_index.out
+++ b/src/test/regress/expected/create_index.out
@@ -53,6 +53,8 @@ CREATE INDEX gpolygonind ON polygon_tbl USING gist (f1);
 CREATE INDEX gcircleind ON circle_tbl USING gist (f1);
 CREATE TEMP TABLE gpolygon_tbl AS
     SELECT polygon(home_base) AS f1 FROM slow_emp4000;
+INSERT INTO gpolygon_tbl VALUES ( '(1000,0,0,1000)' ); 
+INSERT INTO gpolygon_tbl VALUES ( '(0,1000,1000,1000)' ); 
 CREATE TEMP TABLE gcircle_tbl AS
     SELECT circle(home_base) AS f1 FROM slow_emp4000;
 CREATE INDEX ggpolygonind ON gpolygon_tbl USING gist (f1);

diff --git a/src/test/regress/expected/polygon.out b/src/test/regress/expected/polygon.out
index 56b72aaa8e12b9a436a657c30bbf75913a4fea9b..7e0ae242664dbe11dd314f0e56e2b9eb83d18be3 100644 (file)
--- a/src/test/regress/expected/polygon.out
+++ b/src/test/regress/expected/polygon.out
@@ -180,6 +180,59 @@ SELECT polygon '(2.0,0.0),(2.0,4.0),(0.0,0.0)' @> polygon '(3.0,1.0),(3.0,3.0),(
  f
 (1 row)
 
+--     +------------------------+
+--     |    *---*               1
+--     |  + |   | 
+--     |  2 *---*
+--     +------------------------+
+--                              3
+--     endpoints '+' is ofr one polygon, '*' - for another
+--     Edges 1-2, 2-3 are not shown on picture
+SELECT '((0,4),(6,4),(1,2),(6,0),(0,0))'::polygon @> '((2,1),(2,3),(3,3),(3,1))'::polygon AS "false";
+ false 
+-------
+ f
+(1 row)
+
+--     +-----------+
+--     |    *---* /      
+--     |    |   |/ 
+--     |    |   +  
+--     |    |   |\ 
+--     |    *---* \
+--     +-----------+
+SELECT '((0,4),(6,4),(3,2),(6,0),(0,0))'::polygon @> '((2,1),(2,3),(3,3),(3,1))'::polygon AS "true";
+ true 
+------
+ t
+(1 row)
+
+--     +-----------------+
+--     |                 |
+--     |   +---*---*-----+
+--     |   |   |   |
+--     |   +---*---*-----+
+--     |                 |
+--     +-----------------+
+SELECT '((1,1),(1,4),(5,4),(5,3),(2,3),(2,2),(5,2),(5,1))'::polygon @> '((3,2),(3,3),(4,3),(4,2))'::polygon AS "false";
+ false 
+-------
+ f
+(1 row)
+
+--     +---------+
+--     |         |
+--     |    *----*
+--     |    |    |
+--     |    *----*
+--     |         |
+--     +---------+
+SELECT '((0,0),(0,3),(3,3),(3,0))'::polygon @> '((2,1),(2,2),(3,2),(3,1))'::polygon AS "true";
+ true 
+------
+ t
+(1 row)
+
 -- same 
 SELECT polygon '(2.0,0.0),(2.0,4.0),(0.0,0.0)' ~= polygon '(3.0,1.0),(3.0,3.0),(1.0,0.0)' AS false;
  false 
@@ -194,3 +247,34 @@ SELECT polygon '(2.0,0.0),(2.0,4.0),(0.0,0.0)' && polygon '(3.0,1.0),(3.0,3.0),(
  t
 (1 row)
 
+--     +--------------------+
+--     |    *---*              1
+--     |  + |   | 
+--     |  2 *---*
+--     +--------------------+
+--                             3
+--     Edges 1-2, 2-3 are not shown on picture
+SELECT '((0,4),(6,4),(1,2),(6,0),(0,0))'::polygon && '((2,1),(2,3),(3,3),(3,1))'::polygon AS "true";
+ true 
+------
+ t
+(1 row)
+
+--     +--+ *--*
+--     |  | |  |
+--     |  | *--*
+--     |  +----+
+--     |       |
+--     +-------+
+SELECT '((1,4),(1,1),(4,1),(4,2),(2,2),(2,4),(1,4))'::polygon && '((3,3),(4,3),(4,4),(3,4),(3,3))'::polygon AS "false";
+ false 
+-------
+ f
+(1 row)
+
+SELECT '((200,800),(800,800),(800,200),(200,200))' &&  '(1000,1000,0,0)'::polygon AS "true";
+ true 
+------
+ t
+(1 row)
+

diff --git a/src/test/regress/sql/create_index.sql b/src/test/regress/sql/create_index.sql
index 66ab0f9f02fc01d8016d7be2c30cdd2568a3306a..9527ab7a7bff767d38f93d8bafdaa422aa13eba8 100644 (file)
--- a/src/test/regress/sql/create_index.sql
+++ b/src/test/regress/sql/create_index.sql
@@ -78,6 +78,8 @@ CREATE INDEX gcircleind ON circle_tbl USING gist (f1);
 
 CREATE TEMP TABLE gpolygon_tbl AS
     SELECT polygon(home_base) AS f1 FROM slow_emp4000;
+INSERT INTO gpolygon_tbl VALUES ( '(1000,0,0,1000)' ); 
+INSERT INTO gpolygon_tbl VALUES ( '(0,1000,1000,1000)' ); 
 
 CREATE TEMP TABLE gcircle_tbl AS
     SELECT circle(home_base) AS f1 FROM slow_emp4000;

diff --git a/src/test/regress/sql/polygon.sql b/src/test/regress/sql/polygon.sql
index 1f45de0a6d96914970e7baee1f77511a9439c8ae..743ca22f536a29133893ee923174a79c99e4c796 100644 (file)
--- a/src/test/regress/sql/polygon.sql
+++ b/src/test/regress/sql/polygon.sql
@@ -111,9 +111,64 @@ SELECT polygon '(2.0,0.0),(2.0,4.0),(0.0,0.0)' <@ polygon '(3.0,1.0),(3.0,3.0),(
 -- contains 
 SELECT polygon '(2.0,0.0),(2.0,4.0),(0.0,0.0)' @> polygon '(3.0,1.0),(3.0,3.0),(1.0,0.0)' AS false;
 
+--     +------------------------+
+--     |    *---*               1
+--     |  + |   | 
+--     |  2 *---*
+--     +------------------------+
+--                              3
+--     endpoints '+' is ofr one polygon, '*' - for another
+--     Edges 1-2, 2-3 are not shown on picture
+SELECT '((0,4),(6,4),(1,2),(6,0),(0,0))'::polygon @> '((2,1),(2,3),(3,3),(3,1))'::polygon AS "false";
+
+--     +-----------+
+--     |    *---* /      
+--     |    |   |/ 
+--     |    |   +  
+--     |    |   |\ 
+--     |    *---* \
+--     +-----------+
+SELECT '((0,4),(6,4),(3,2),(6,0),(0,0))'::polygon @> '((2,1),(2,3),(3,3),(3,1))'::polygon AS "true";
+
+--     +-----------------+
+--     |                 |
+--     |   +---*---*-----+
+--     |   |   |   |
+--     |   +---*---*-----+
+--     |                 |
+--     +-----------------+
+SELECT '((1,1),(1,4),(5,4),(5,3),(2,3),(2,2),(5,2),(5,1))'::polygon @> '((3,2),(3,3),(4,3),(4,2))'::polygon AS "false";
+
+--     +---------+
+--     |         |
+--     |    *----*
+--     |    |    |
+--     |    *----*
+--     |         |
+--     +---------+
+SELECT '((0,0),(0,3),(3,3),(3,0))'::polygon @> '((2,1),(2,2),(3,2),(3,1))'::polygon AS "true";
+
 -- same 
 SELECT polygon '(2.0,0.0),(2.0,4.0),(0.0,0.0)' ~= polygon '(3.0,1.0),(3.0,3.0),(1.0,0.0)' AS false;
 
 -- overlap 
 SELECT polygon '(2.0,0.0),(2.0,4.0),(0.0,0.0)' && polygon '(3.0,1.0),(3.0,3.0),(1.0,0.0)' AS true;
 
+--     +--------------------+
+--     |    *---*              1
+--     |  + |   | 
+--     |  2 *---*
+--     +--------------------+
+--                             3
+--     Edges 1-2, 2-3 are not shown on picture
+SELECT '((0,4),(6,4),(1,2),(6,0),(0,0))'::polygon && '((2,1),(2,3),(3,3),(3,1))'::polygon AS "true";
+
+--     +--+ *--*
+--     |  | |  |
+--     |  | *--*
+--     |  +----+
+--     |       |
+--     +-------+
+SELECT '((1,4),(1,1),(4,1),(4,2),(2,2),(2,4),(1,4))'::polygon && '((3,3),(4,3),(4,4),(3,4),(3,3))'::polygon AS "false";
+SELECT '((200,800),(800,800),(800,200),(200,200))' &&  '(1000,1000,0,0)'::polygon AS "true";
+