Skip to content

Commit

Permalink
Merge #62643
Browse files Browse the repository at this point in the history 
62643: geo/geomfn: add point in polygon optimization to st_contains, st_within r=otan a=andyyang890

This patch improves the performance of st_contains
and st_within for the common use case of testing
whether a (multi)polygon contains a (multi)point.

Release note: None

Co-authored-by: Andy Yang <[email protected]>
  • Loading branch information
craig[bot] and Andy Yang committed Mar 30, 2021
2 parents aec84ca + b3dcc72 commit ed698ae
Show file tree
Hide file tree
Showing 3 changed files with 104 additions and 25 deletions.
28 changes: 26 additions & 2 deletions pkg/geo/geomfn/binary_predicates.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,8 +34,8 @@ func Covers(a geo.Geometry, b geo.Geometry) (bool, error) {
// A point cannot cover a polygon.
return false, nil
case PolygonAndPoint:
// Computing whether a polygon covers a point is the same
// as computing whether a point is covered by the polygon.
// Computing whether a polygon covers a point is equivalent
// to computing whether the point is covered by the polygon.
return PointKindRelatesToPolygonKind(pointKind, polygonKind, PointPolygonCoveredBy)
}

Expand Down Expand Up @@ -72,6 +72,19 @@ func Contains(a geo.Geometry, b geo.Geometry) (bool, error) {
if !a.CartesianBoundingBox().Covers(b.CartesianBoundingBox()) {
return false, nil
}

// Optimization for point in polygon calculations.
pointPolygonPair, pointKind, polygonKind := PointKindAndPolygonKind(a, b)
switch pointPolygonPair {
case PointAndPolygon:
// A point cannot contain a polygon.
return false, nil
case PolygonAndPoint:
// Computing whether a polygon contains a point is equivalent
// to computing whether the point is contained within the polygon.
return PointKindRelatesToPolygonKind(pointKind, polygonKind, PointPolygonWithin)
}

return geos.Contains(a.EWKB(), b.EWKB())
}

Expand Down Expand Up @@ -314,5 +327,16 @@ func Within(a geo.Geometry, b geo.Geometry) (bool, error) {
if !b.CartesianBoundingBox().Covers(a.CartesianBoundingBox()) {
return false, nil
}

// Optimization for point in polygon calculations.
pointPolygonPair, pointKind, polygonKind := PointKindAndPolygonKind(a, b)
switch pointPolygonPair {
case PolygonAndPoint:
// A polygon cannot be contained within a point.
return false, nil
case PointAndPolygon:
return PointKindRelatesToPolygonKind(pointKind, polygonKind, PointPolygonWithin)
}

return geos.Within(a.EWKB(), b.EWKB())
}
3 changes: 3 additions & 0 deletions pkg/geo/geomfn/binary_predicates_test.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -413,6 +413,9 @@ func TestWithin(t *testing.T) {
{leftRectPoint, leftRect, true},
{leftRectMultiPoint, leftRect, true},
{leftRectMultiPoint, leftRectWithHole, false},
{leftRectHoleEdgePoint, leftRectWithHole, false},
{leftRectEdgePoint, leftRectWithHole, false},
{leftRectEdgeMultiPoint, leftRectWithHole, true},
}

for i, tc := range testCases {
Expand Down
98 changes: 75 additions & 23 deletions pkg/geo/geomfn/distance.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -847,6 +847,9 @@ const (
// PointPolygonCoveredBy is the relationship where a (multi)point
// is covered by a (multi)polygon.
PointPolygonCoveredBy
// PointPolygonWithin is the relationship where a (multi)point
// is contained within a (multi)polygon.
PointPolygonWithin
)

// PointKindRelatesToPolygonKind returns whether a (multi)point and
Expand All @@ -865,10 +868,14 @@ func PointKindRelatesToPolygonKind(
pointKindIterator := geo.NewGeomTIterator(pointKindBaseT, geo.EmptyBehaviorOmit)
polygonKindIterator := geo.NewGeomTIterator(polygonKindBaseT, geo.EmptyBehaviorOmit)

// TODO(ayang): Think about how to refactor these nested for loops
// Check whether each point intersects with at least one polygon.
// For Intersects, at least one point must intersect with at least one polygon.
// For CoveredBy, every point must intersect with at least one polygon.
// - For Intersects, at least one point must intersect with at least one polygon.
// - For CoveredBy, every point must intersect with at least one polygon.
// - For Within, every point must intersect with at least one polygon
// and at least one point must be inside at least one polygon.
intersectsOnce := false
insideOnce := false
pointOuterLoop:
for {
point, hasPoint, err := pointKindIterator.Next()
Expand All @@ -880,6 +887,7 @@ pointOuterLoop:
}
// Reset the polygon iterator on each iteration of the point iterator.
polygonKindIterator.Reset()
curIntersects := false
for {
polygon, hasPolygon, err := polygonKindIterator.Next()
if err != nil {
Expand All @@ -888,22 +896,40 @@ pointOuterLoop:
if !hasPolygon {
break
}
intersects, err := pointIntersectsPolygon(point, polygon)
pointSide, err := findPointSideOfPolygon(point, polygon)
if err != nil {
return false, err
}
if intersects {
switch pointSide {
case insideLinearRing:
insideOnce = true
switch relationType {
case PointPolygonWithin:
continue pointOuterLoop
}
fallthrough
case onLinearRing:
intersectsOnce = true
curIntersects = true
switch relationType {
case PointPolygonIntersects:
// A single intersection is sufficient.
return true, nil
case PointPolygonCoveredBy:
// If the current point intersects, check the next point.
continue pointOuterLoop
case PointPolygonWithin:
// We can only skip to the next point if we have already seen a point
// that is inside the (multi)polygon.
if insideOnce {
continue pointOuterLoop
}
default:
return false, errors.Newf("unknown PointPolygonRelationType")
}
case outsideLinearRing:
default:
return false, errors.Newf("findPointSideOfPolygon returned unknown linearRingSide %d", pointSide)
}
}
// Case where a point in the (multi)point does not intersect
Expand All @@ -913,61 +939,87 @@ pointOuterLoop:
// Each point in a (multi)point must intersect a polygon in the
// (multi)point to be covered by it.
return false, nil
case PointPolygonWithin:
if !curIntersects {
return false, nil
}
}
}
return intersectsOnce, nil
switch relationType {
case PointPolygonCoveredBy:
return intersectsOnce, nil
case PointPolygonWithin:
return insideOnce, nil
default:
return false, nil
}
}

// pointIntersectsPolygon returns whether a point intersects with a polygon.
func pointIntersectsPolygon(point geom.T, polygon geom.T) (bool, error) {
// findPointSideOfPolygon returns whether a point intersects with a polygon.
func findPointSideOfPolygon(point geom.T, polygon geom.T) (linearRingSide, error) {
// Convert point from a geom.T to a *geodist.Point.
_, ok := point.(*geom.Point)
if !ok {
return false, errors.Newf("first geometry passed to PointIntersectsPolygon must be a point")
return outsideLinearRing, errors.Newf("first geometry passed to findPointSideOfPolygon must be a point")
}
pointGeodistShape, err := geomToGeodist(point)
if err != nil {
return false, err
return outsideLinearRing, err
}
pointGeodistPoint, ok := pointGeodistShape.(*geodist.Point)
if !ok {
return false, errors.Newf("geomToGeodist failed to convert a *geom.Point to a *geodist.Point")
return outsideLinearRing, errors.Newf("geomToGeodist failed to convert a *geom.Point to a *geodist.Point")
}

// Convert polygon from a geom.T to a geodist.Polygon.
_, ok = polygon.(*geom.Polygon)
if !ok {
return false, errors.Newf("second geometry passed to PointIntersectsPolygon must be a polygon")
return outsideLinearRing, errors.Newf("second geometry passed to findPointSideOfPolygon must be a polygon")
}
polygonGeodistShape, err := geomToGeodist(polygon)
if err != nil {
return false, err
return outsideLinearRing, err
}
polygonGeodistPolygon, ok := polygonGeodistShape.(geodist.Polygon)
if !ok {
return false, errors.Newf("geomToGeodist failed to convert a *geom.Polygon to a geodist.Polygon")
return outsideLinearRing, errors.Newf("geomToGeodist failed to convert a *geom.Polygon to a geodist.Polygon")
}

// Point cannot be inside an empty polygon.
if polygonGeodistPolygon.NumLinearRings() == 0 {
return false, nil
return outsideLinearRing, nil
}

// Check if point is inside main outer boundary of polygon.
// If it isn't, then it's not inside the polygon.
// Find which side the point is relative to the main outer boundary of
// the polygon. If it outside or on the boundary, we can conclude
// that the point is on that side of the overall polygon as well.
mainRing := polygonGeodistPolygon.LinearRing(0)
if findPointSideOfLinearRing(*pointGeodistPoint, mainRing) == outsideLinearRing {
return false, nil
switch pointSide := findPointSideOfLinearRing(*pointGeodistPoint, mainRing); pointSide {
case insideLinearRing:
case outsideLinearRing, onLinearRing:
return pointSide, nil
default:
return outsideLinearRing, errors.Newf("unknown linearRingSide %d", pointSide)
}

// If the point is indeed inside the main outer boundary of the polygon,
// it must also not be in the interior of any holes to be inside the polygon.
// If the point is inside the main outer boundary of the polygon, we must
// determine which side it is relative to every hole in the polygon.
// If it is inside any hole, it is outside the polygon. If it is on the
// ring of any hole, it is on the boundary of the polygon. Otherwise,
// it is inside the polygon.
for ringNum := 1; ringNum < polygonGeodistPolygon.NumLinearRings(); ringNum++ {
polygonHole := polygonGeodistPolygon.LinearRing(ringNum)
if findPointSideOfLinearRing(*pointGeodistPoint, polygonHole) == insideLinearRing {
return false, nil
switch pointSide := findPointSideOfLinearRing(*pointGeodistPoint, polygonHole); pointSide {
case insideLinearRing:
return outsideLinearRing, nil
case onLinearRing:
return onLinearRing, nil
case outsideLinearRing:
continue
default:
return outsideLinearRing, errors.Newf("unknown linearRingSide %d", pointSide)
}
}

return true, nil
return insideLinearRing, nil
}

0 comments on commit ed698ae

Please sign in to comment.