diff options
Diffstat (limited to 'src/backend/utils/adt/geo_ops.c')
| -rw-r--r-- | src/backend/utils/adt/geo_ops.c | 252 |
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; } |