8000 Correct calculations of overlap and contains operations over polygons. · home201448/postgres@49475aa · GitHub
[go: up one dir, main page]
Skip to content

Commit 49475aa

Browse files
feodorfeodor
committed
Correct calculations of overlap and contains operations over polygons.
1 parent 1f4b046 commit 49475aa

File tree

5 files changed

+346
-49
lines changed

5 files changed

+346
-49
lines changed

‎src/backend/utils/adt/geo_ops.c

Lines changed: 203 additions & 49 deletions
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@
88
*
99
*
1010
* IDENTIFICATION
11-
* $PostgreSQL: pgsql/src/backend/utils/adt/geo_ops.c,v 1.102 2009/06/23 16:25:02 tgl Exp $
11+
* $PostgreSQL: pgsql/src/backend/utils/adt/geo_ops.c,v 1.103 2009/07/28 09:47:59 teodor Exp $
1212
*
1313
*-------------------------------------------------------------------------
1414
*/
@@ -66,6 +66,8 @@ static bool has_interpt_sl(LSEG *lseg, LINE *line);
6666
static double dist_pl_internal(Point *pt, LINE *line);
6767
static double dist_ps_internal(Point *pt, LSEG *lseg);
6868
static Point *line_interpt_internal(LINE *l1, LINE *l2);
69+
static bool lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start);
70+
static Point* lseg_interpt_internal(LSEG *l1, LSEG *l2);
6971

7072

7173
/*
@@ -2352,15 +2354,9 @@ lseg_center(PG_FUNCTION_ARGS)
23522354
PG_RETURN_POINT_P(result);
23532355
}
23542356

2355-
2356-
/* lseg_interpt -
2357-
* Find the intersection point of two segments (if any).
2358-
*/
2359-
Datum
2360-
lseg_interpt(PG_FUNCTION_ARGS)
2357+
static Point*
2358+
lseg_interpt_internal(LSEG *l1, LSEG *l2)
23612359
{
2362-
LSEG *l1 = PG_GETARG_LSEG_P(0);
2363-
LSEG *l2 = PG_GETARG_LSEG_P(1);
23642360
Point *result;
23652361
LINE tmp1,
23662362
tmp2;
@@ -2372,15 +2368,18 @@ lseg_interpt(PG_FUNCTION_ARGS)
23722368
line_construct_pts(&tmp2, &l2->p[0], &l2->p[1]);
23732369
result = line_interpt_internal(&tmp1, &tmp2);
23742370
if (!PointerIsValid(result))
2375-
PG_RETURN_NULL();
2371+
return NULL;
23762372

23772373
/*
23782374
* If the line intersection point isn't within l1 (or equivalently l2),
23792375
* there is no valid segment intersection point at all.
23802376
*/
23812377
if (!on_ps_internal(result, l1) ||
23822378
!on_ps_internal(result, l2))
2383-
PG_RETURN_NULL();
2379+
{
2380+
pfree(result);
2381+
return NULL;
2382+
}
23842383

23852384
/*
23862385
* If there is an intersection, then check explicitly for matching
@@ -2400,6 +2399,23 @@ lseg_interpt(PG_FUNCTION_ARGS)
24002399
result->y = l1->p[1].y;
24012400
}
24022401

2402+
return result;
2403+
}
2404+
2405+
/* lseg_interpt -
2406+
* Find the intersection point of two segments (if any).
2407+
*/
2408+
Datum
2409+
lseg_interpt(PG_FUNCTION_ARGS)
2410+
{
2411+
LSEG *l1 = PG_GETARG_LSEG_P(0);
2412+
LSEG *l2 = PG_GETARG_LSEG_P(1);
2413+
Point *result;
2414+
2415+
result = lseg_interpt_internal(l1, l2);
2416+
if (!PointerIsValid(result))
2417+
PG_RETURN_NULL();
2418+
24032419
PG_RETURN_POINT_P(result);
24042420
}
24052421

@@ -3742,10 +3758,7 @@ poly_same(PG_FUNCTION_ARGS)
37423758
}
37433759

37443760
/*-----------------------------------------------------------------
3745-
* Determine if polygon A overlaps polygon B by determining if
3746-
* their bounding boxes overlap.
3747-
*
3748-
* XXX ought to do a more correct check!
3761+
* Determine if polygon A overlaps polygon B
37493762
*-----------------------------------------------------------------*/
37503763
Datum
37513764
poly_overlap(PG_FUNCTION_ARGS)
@@ -3754,7 +3767,54 @@ poly_overlap(PG_FUNCTION_ARGS)
37543767
POLYGON *polyb = PG_GETARG_POLYGON_P(1);
37553768
bool result;
37563769

3757-
result = box_ov(&polya->boundbox, &polyb->boundbox);
3770+
/* Quick check by bounding box */
3771+
result = (polya->npts > 0 && polyb->npts > 0 &&
3772+
box_ov(&polya->boundbox, &polyb->boundbox)) ? true : false;
3773+
3774+
/*
3775+
* Brute-force algorithm - try to find intersected edges,
3776+
* if so then polygons are overlapped else check is one
3777+
* polygon inside other or not by testing single point
3778+
* of them.
3779+
*/
3780+
if (result)
3781+
{
3782+
int ia, ib;
3783+
LSEG sa, sb;
3784+
3785+
/* Init first of polya's edge with last point */
3786+
sa.p[0] = polya->p[polya->npts - 1];
3787+
result = false;
3788+
3789+
for(ia=0; ia<polya->npts && result == false; ia++)
3790+
{
3791+
/* Second point of polya's edge is a current one */
3792+
sa.p[1] = polya->p[ia];
3793+
3794+
/* Init first of polyb's edge with last point */
3795+
sb.p[0] = polyb->p[polyb->npts - 1];
3796+
3797+
for(ib=0; ib<polyb->npts && result == false; ib++)
3798+
{
3799+
sb.p[1] = polyb->p[ib];
3800+
result = lseg_intersect_internal(&sa, &sb);
3801+
sb.p[0] = sb.p[1];
3802+
}
3803+
3804+
/*
3805+
* move current endpoint to the first point
3806+
* of next edge
3807+
*/
3808+
sa.p[0] = sa.p[1];
3809+
}
3810+
3811+
if (result==false)
3812+
{
3813+
result = ( point_inside(polya->p, polyb->npts, polyb->p)
3814+
||
3815+
point_inside(polyb->p, polya->npts, polya->p) );
3816+
}
3817+
}
37583818

37593819
/*
37603820
* Avoid leaking memory for toasted inputs ... needed for rtree indexes
@@ -3765,6 +3825,119 @@ poly_overlap(PG_FUNCTION_ARGS)
37653825
PG_RETURN_BOOL(result);
37663826
}
37673827

3828+
/*
3829+
* Tests special kind of segment for in/out of polygon.
3830+
* Special kind means:
3831+
* - point a should be on segment s
3832+
* - segment (a,b) should not be contained by s
3833+
* Returns true if:
3834+
* - segment (a,b) is collinear to s and (a,b) is in polygon
3835+
* - segment (a,b) s not collinear to s. Note: that doesn't
3836+
* mean that segment is in polygon!
3837+
*/
3838+
3839+
static bool
3840+
touched_lseg_inside_poly(Point *a, Point *b, LSEG *s, POLYGON *poly, int start)
3841+
{
3842+
/* point a is on s, b is not */
3843+
LSEG t;
3844+
3845+
t.p[0] = *a;
3846+
t.p[1] = *b;
3847+
3848+
#define POINTEQ(pt1, pt2) (FPeq((pt1)->x, (pt2)->x) && FPeq((pt1)->y, (pt2)->y))
3849+
if ( POINTEQ(a, s->p) )
3850+
{
3851+
if ( on_ps_internal(s->p+1, &t) )
3852+
return lseg_inside_poly(b, s->p+1, poly, start);
3853+
}
3854+
else if (POINTEQ(a, s->p+1))
3855+
{
3856+
if ( on_ps_internal(s->p, &t) )
3857+
return lseg_inside_poly(b, s->p, poly, start);
3858+
}
3859+
else if ( on_ps_internal(s->p, &t) )
3860+
{
3861+
return lseg_inside_poly(b, s->p, poly, start);
3862+
}
3863+
else if ( on_ps_internal(s->p+1, &t) )
3864+
{
3865+
return lseg_inside_poly(b, s->p+1, poly, start);
3866+
}
3867+
3868+
return true; /* may be not true, but that will check later */
3869+
}
3870+
3871+
/*
3872+
* Returns true if segment (a,b) is in polygon, option
3873+
* start is used for optimization - function checks
3874+
* polygon's edges started from start
3875+
*/
3876+
static bool
3877+
lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start)
3878+
{
3879+
LSEG s,
3880+
t;
3881+
int i;
3882+
bool res = true,
3883+
intersection = false;
3884+
3885+
t.p[0] = *a;
3886+
t.p[1] = *b;
3887+
s.p[0] = poly->p[( start == 0) ? (poly->npts - 1) : (start - 1)];
3888+
3889+
for(i=start; i<poly->npts && res == true; i++)
3890+
{
3891+
Point *interpt;
3892+
3893+
s.p[1] = poly->p[i];
3894+
3895+
if ( on_ps_internal(t.p, &s) )
3896+
{
3897+
if ( on_ps_internal(t.p+1, &s) )
3898+
return true; /* t is contained by s */
3899+
3900+
/* Y-cross */
3901+
res = touched_lseg_inside_poly(t.p, t.p+1, &s, poly, i+1);
3902+
}
3903+
else if ( on_ps_internal(t.p+1, &s) )
< F438 /code>3904+
{
3905+
/* Y-cross */
3906+
res = touched_lseg_inside_poly(t.p+1, t.p, &s, poly, i+1);
3907+
}
3908+
else if ( (interpt = lseg_interpt_internal(&t, &s)) != NULL )
3909+
{
3910+
/*
3911+
* segments are X-crossing, go to check each subsegment
3912+
*/
3913+
3914+
intersection = true;
3915+
res = lseg_inside_poly(t.p, interpt, poly, i+1);
3916+
if (res)
3917+
res = lseg_inside_poly(t.p+1, interpt, poly, i+1);
3918+
pfree(interpt);
3919+
}
3920+
3921+
s.p[0] = s.p[1];
3922+
}
3923+
3924+
if (res && !intersection)
3925+
{
3926+
Point p;
3927+
3928+
/*
3929+
* if X-intersection wasn't found then check central point
3930+
* of tested segment. In opposite case we already check all
3931+
* subsegments
3932+
*/
3933+
p.x = (t.p[0].x + t.p[1].x) / 2.0;
3934+
p.y = (t.p[0].y + t.p[1].y) / 2.0;
3935+
3936+
res = point_inside(&p, poly->npts, poly->p);
3937+
}
3938+
3939+
return res;
3940+
}
37683941

37693942
/*-----------------------------------------------------------------
37703943
* Determine if polygon A contains polygon B.
@@ -3775,49 +3948,30 @@ poly_contain(PG_FUNCTION_ARGS)
37753948
POLYGON *polya = PG_GETARG_POLYGON_P(0);
37763949
POLYGON *polyb = PG_GETARG_POLYGON_P(1);
37773950
bool result;
3778-
int i;
37793951

37803952
/*
37813953
* Quick check to see if bounding box is contained.
37823954
*/
3783-
if (DatumGetBool(DirectFunctionCall2(box_contain,
3784-
BoxPGetDatum(&polya->boundbox),
3785-
BoxPGetDatum(&polyb->boundbox))))
3955+
if (polya->npts > 0 && polyb->npts > 0 &&
3956+
DatumGetBool(DirectFunctionCall2(box_contain,
3957+
BoxPGetDatum(&polya->boundbox),
3958+
BoxPGetDatum(&polyb->boundbox))))
37863959
{
3787-
result = true; /* assume true for now */
3788-
for (i = 0; i < polyb->npts; i++)
3789-
{
3790-
if (point_inside(&(polyb->p[i]), polya->npts, &(polya->p[0])) == 0)
3791-
{
3792-
#ifdef GEODEBUG
3793-
printf("poly_contain- point (%f,%f) not in polygon\n", polyb->p[i].x, polyb->p[i].y);
3794-
#endif
3795-
result = false;
3796-
break;
3797-
}
3798-
}
3799-
if (result)
3960+
int i;
3961+
LSEG s;
3962+
3963+
s.p[0] = polyb->p[polyb->npts - 1];
3964+
result = true;
3965+
3966+
for(i=0; i<polyb->npts && result == true; i++)
38003967
{
3801-
for (i = 0; i < polya->npts; i++)
3802-
{
3803-
if (point_inside(&(polya->p[i]), polyb->npts, &(polyb->p[0])) == 1)
3804-
{
3805-
#ifdef GEODEBUG
3806-
printf("poly_contain- point (%f,%f) in polygon\n", polya->p[i].x, polya->p[i].y);
3807-
#endif
3808-
result = false;
3809-
break;
3810-
}
3811-
}
3968+
s.p[1] = polyb->p[i];
3969+
result = lseg_inside_poly(s.p, s.p+1, polya, 0);
3970+
s.p[0] = s.p[1];
38123971
}
38133972
}
38143973
else
38153974
{
3816-
#ifdef GEODEBUG
3817-
printf("poly_contain- bound box ((%f,%f),(%f,%f)) not inside ((%f,%f),(%f,%f))\n",
3818-
polyb->boundbox.low.x, polyb->boundbox.low.y, polyb->boundbox.high.x, polyb->boundbox.high.y,
3819-
polya->boundbox.low.x, polya->boundbox.low.y, polya->boundbox.high.x, polya->boundbox.high.y);
3820-
#endif
38213975
result = false;
38223976
}
38233977

‎src/test/regress/expected/create_index.out

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -53,6 +53,8 @@ CREATE INDEX gpolygonind ON polygon_tbl USING gist (f1);
5353
CREATE INDEX gcircleind ON circle_tbl USING gist (f1);
5454
CREATE TEMP TABLE gpolygon_tbl AS
5555
SELECT polygon(home_base) AS f1 FROM slow_emp4000;
56+
INSERT INTO gpolygon_tbl VALUES ( '(1000,0,0,1000)' );
57+
INSERT INTO gpolygon_tbl VALUES ( '(0,1000,1000,1000)' );
5658
CREATE TEMP TABLE gcircle_tbl AS
5759
SELECT circle(home_base) AS f1 FROM slow_emp4000;
5860
CREATE INDEX ggpolygonind ON gpolygon_tbl USING gist (f1);

0 commit comments

Comments
 (0)
0