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);
/*
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;
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),
*/
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
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);
}
}
/*-----------------------------------------------------------------
- * 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)
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
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.
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;
}
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
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)
+
-- 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";
+