From 49475aab8d0da9077a5828d404ba1e15c3e093c2 Mon Sep 17 00:00:00 2001 From: Teodor Sigaev Date: Tue, 28 Jul 2009 09:48:00 +0000 Subject: [PATCH] Correct calculations of overlap and contains operations over polygons. --- src/backend/utils/adt/geo_ops.c | 252 +++++++++++++++++---- src/test/regress/expected/create_index.out | 2 + src/test/regress/expected/polygon.out | 84 +++++++ src/test/regress/sql/create_index.sql | 2 + src/test/regress/sql/polygon.sql | 55 +++++ 5 files changed, 346 insertions(+), 49 deletions(-) diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c index 24157a53c5..ea7279e8c7 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; ianpts && 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; ibnpts && 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; inpts && 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; inpts && 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 7b37cbf40b..6e06be5ed5 100644 --- 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 56b72aaa8e..7e0ae24266 100644 --- 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 66ab0f9f02..9527ab7a7b 100644 --- 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 1f45de0a6d..743ca22f53 100644 --- 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"; + -- 2.40.0