Correct calculations of overlap and contains operations over polygons.

1 files changed, 203 insertions, 49 deletions

diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c
index 24157a53c59..ea7279e8c7c 100644
--- 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;
 	}