Correct calculations of overlap and contains operations over polygons.
This commit is contained in:
Teodor Sigaev
parent
1f4b046c18
commit
49475aab8d
|
|
@ -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++)
|
||||
int i;
|
||||
LSEG s;
|
||||
|
||||
s.p[0] = polyb->p[polyb->npts - 1];
|
||||
result = true;
|
||||
|
||||
for(i=0; i<polyb->npts && result == true; 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)
|
||||
{
|
||||
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;
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -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";
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue