geo/geomfn: implement ST_MaxDistance / ST_DFullyWithin

docs/generated/sql/functions.md

@@ -750,6 +750,8 @@ has no relationship with the commit order of concurrent transactions.</p>
 <p>This function utilizes the GEOS module.</p>
 <p>This function will automatically use any available index.</p>
 </span></td></tr>
+<tr><td><a name="st_dfullywithin"></a><code>st_dfullywithin(geometry_a: geometry, geometry_b: geometry, distance: <a href="float.html">float</a>) &rarr; <a href="bool.html">bool</a></code></td><td><span class="funcdesc"><p>Returns true if every pair of points comprising geometry_a and geometry_b are within distance units. In other words, the ST_MaxDistance between geometry_a and geometry_b is less than or equal to distance units.</p>
+</span></td></tr>
 <tr><td><a name="st_distance"></a><code>st_distance(geography_a: geography, geography_b: geography) &rarr; <a href="float.html">float</a></code></td><td><span class="funcdesc"><p>Returns the distance in meters between geography_a and geography_b.  Uses a spheroid to perform the operation.&quot;\n\nWhen operating on a spheroid, this function will use the sphere to calculate the closest two points using S2. The spheroid distance between these two points is calculated using GeographicLib. This follows observed PostGIS behavior.</p>
 <p>This function utilizes the GeographicLib library for spheroid calculations.</p>
 <p>This function will automatically use any available index.</p>
@@ -865,6 +867,8 @@ has no relationship with the commit order of concurrent transactions.</p>
 </span></td></tr>
 <tr><td><a name="st_makepoint"></a><code>st_makepoint(x: <a href="float.html">float</a>, y: <a href="float.html">float</a>) &rarr; geometry</code></td><td><span class="funcdesc"><p>Returns a new Point with the given X and Y coordinates.</p>
 </span></td></tr>
+<tr><td><a name="st_maxdistance"></a><code>st_maxdistance(geometry_a: geometry, geometry_b: geometry) &rarr; <a href="float.html">float</a></code></td><td><span class="funcdesc"><p>Returns the maximum distance across every pair of points comprising the given geometries. Note if the geometries are the same, it will return the maximum distance between the geometry’s vertexes.</p>
+</span></td></tr>
 <tr><td><a name="st_mlinefromtext"></a><code>st_mlinefromtext(str: <a href="string.html">string</a>, srid: <a href="int.html">int</a>) &rarr; geometry</code></td><td><span class="funcdesc"><p>Returns the Geometry from a WKT or EWKT representation with an SRID. If the shape underneath is not MultiLineString, NULL is returned. If the SRID is present in both the EWKT and the argument, the argument value is used.</p>
 </span></td></tr>
 <tr><td><a name="st_mlinefromtext"></a><code>st_mlinefromtext(val: <a href="string.html">string</a>) &rarr; geometry</code></td><td><span class="funcdesc"><p>Returns the Geometry from a WKT or EWKT representation. If the shape underneath is not MultiLineString, NULL is returned.</p>

pkg/geo/geodist/geodist.go

@@ -78,6 +78,8 @@ type DistanceUpdater interface {
 	OnIntersects() bool
 	// Distance returns the distance to return so far.
 	Distance() float64
+	// IsMaxDistance returns whether the updater is looking for maximum distance.
+	IsMaxDistance() bool
 }
 
 // EdgeCrosser is a provided hook that calculates whether edges intersect.
@@ -150,20 +152,26 @@ func ShapeDistance(c DistanceCalculator, a Shape, b Shape) (bool, error) {
 	return false, errors.Newf("unknown shape: %T", a)
 }
 
-// onPointToEdgesExceptFirstEdgeStart updates the distance using the edges between a point and a shape.
-// It will only check the ends of the edges, and assumes the check against .Edge(0).V0 is not required.
+// onPointToEdgesExceptFirstEdgeStart updates the distance against the edges of a shape and a point.
+// It will only check the V1 of each edge and assumes the first edge start does not need the distance
+// to be computed.
 func onPointToEdgesExceptFirstEdgeStart(c DistanceCalculator, a Point, b shapeWithEdges) bool {
 	for edgeIdx := 0; edgeIdx < b.NumEdges(); edgeIdx++ {
 		edge := b.Edge(edgeIdx)
 		// Check against all V1 of every edge.
 		if c.DistanceUpdater().Update(a, edge.V1) {
 			return true
 		}
-		// Also project the point to the infinite line of the edge, and compare if the closestPoint
-		// lies on the edge.
-		if closestPoint, ok := c.ClosestPointToEdge(edge, a); ok {
-			if c.DistanceUpdater().Update(a, closestPoint) {
-				return true
+		// The max distance between a point and the set of points representing an edge is the
+		// maximum distance from the point and the pair of end-points of the edge, so we don't
+		// need to update the distance using the projected point.
+		if !c.DistanceUpdater().IsMaxDistance() {
+			// Also project the point to the infinite line of the edge, and compare if the closestPoint
+			// lies on the edge.
+			if closestPoint, ok := c.ClosestPointToEdge(edge, a); ok {
+				if c.DistanceUpdater().Update(a, closestPoint) {
+					return true
+				}
 			}
 		}
 	}
@@ -183,9 +191,20 @@ func onPointToLineString(c DistanceCalculator, a Point, b LineString) bool {
 // onPointToPolygon updates the distance between a point and a polygon.
 // Returns true if the calling function should early exit.
 func onPointToPolygon(c DistanceCalculator, a Point, b Polygon) bool {
-	// If the exterior ring does not contain the point, we just need to calculate the distance to
+	// MaxDistance: When computing the maximum distance, the cases are:
+	//   - The point P is not contained in the exterior of the polygon G.
+	//     Say vertex V is the vertex of the exterior of the polygon that is
+	//     furthest away from point P (among all the exterior vertices).
+	//   - One can prove that any vertex of the holes will be closer to point P than vertex V.
+	//     Similarly we can prove that any point in the interior of the polygin is closer to P than vertex V.
+	//     Therefore we only need to compare with the exterior.
+	//   - The point P is contained in the exterior and inside a hole of polygon G.
+	//     One can again prove that the furthest point in the polygon from P is one of the vertices of the exterior.
+	//   - The point P is contained in the polygon. One can again prove the same property.
+	//   So we only need to compare with the exterior ring.
+	// MinDistance: If the exterior ring does not contain the point, we just need to calculate the distance to
 	// the exterior ring.
-	if !c.PointInLinearRing(a, b.LinearRing(0)) {
+	if c.DistanceUpdater().IsMaxDistance() || !c.PointInLinearRing(a, b.LinearRing(0)) {
 		return onPointToEdgesExceptFirstEdgeStart(c, a, b.LinearRing(0))
 	}
 	// At this point it may be inside a hole.
@@ -210,9 +229,13 @@ func onShapeEdgesToShapeEdges(c DistanceCalculator, a shapeWithEdges, b shapeWit
 		crosser := c.NewEdgeCrosser(aEdge, b.Edge(0).V0)
 		for bEdgeIdx := 0; bEdgeIdx < b.NumEdges(); bEdgeIdx++ {
 			bEdge := b.Edge(bEdgeIdx)
-			// If the edges cross, the distance is 0.
-			if crosser.ChainCrossing(bEdge.V1) {
-				return c.DistanceUpdater().OnIntersects()
+			// Max distance between 2 edges is the maximum of the distance across pairs of vertices chosen from each edge.
+			// It does not matter whether the edges cross, so we skip this check.
+			if !c.DistanceUpdater().IsMaxDistance() {
+				// If the edges cross, the distance is 0.
+				if crosser.ChainCrossing(bEdge.V1) {
+					return c.DistanceUpdater().OnIntersects()
+				}
 			}
 
 			// Compare each vertex against the edge of the other.
@@ -230,10 +253,12 @@ func onShapeEdgesToShapeEdges(c DistanceCalculator, a shapeWithEdges, b shapeWit
 					c.DistanceUpdater().Update(toCheck.vertex, toCheck.edge.V1) {
 					return true
 				}
-				// Also check the projection of the vertex onto the edge.
-				if closestPoint, ok := c.ClosestPointToEdge(toCheck.edge, toCheck.vertex); ok {
-					if c.DistanceUpdater().Update(toCheck.vertex, closestPoint) {
-						return true
+				if !c.DistanceUpdater().IsMaxDistance() {
+					// Also check the projection of the vertex onto the edge.
+					if closestPoint, ok := c.ClosestPointToEdge(toCheck.edge, toCheck.vertex); ok {
+						if c.DistanceUpdater().Update(toCheck.vertex, closestPoint) {
+							return true
+						}
 					}
 				}
 			}
@@ -245,14 +270,17 @@ func onShapeEdgesToShapeEdges(c DistanceCalculator, a shapeWithEdges, b shapeWit
 // onLineStringToPolygon updates the distance between a polyline and a polygon.
 // Returns true if the calling function should early exit.
 func onLineStringToPolygon(c DistanceCalculator, a LineString, b Polygon) bool {
-	// If we know at least one point is outside the exterior ring, then there are two cases:
-	// * the line is always outside the exterior ring. We only need to compare the line
-	//   against the exterior ring.
-	// * the line intersects with the exterior ring.
-	// In both these cases, we can defer to the edge to edge comparison between the line
-	// and the exterior ring.
-	// We use the first point of the linestring for this check.
-	if !c.PointInLinearRing(a.Vertex(0), b.LinearRing(0)) {
+	// MinDistance: If we know at least one point is outside the exterior ring, then there are two cases:
+	//   * the line is always outside the exterior ring. We only need to compare the line
+	//     against the exterior ring.
+	//   * the line intersects with the exterior ring.
+	//   In both these cases, we can defer to the edge to edge comparison between the line
+	//   and the exterior ring.
+	//   We use the first point of the linestring for this check.
+	// MaxDistance: the furthest distance from a LineString to a Polygon is always against the
+	//   exterior ring. This follows the reasoning under "onPointToPolygon", but we must now
+	//   check each point in the LineString.
+	if c.DistanceUpdater().IsMaxDistance() || !c.PointInLinearRing(a.Vertex(0), b.LinearRing(0)) {
 		return onShapeEdgesToShapeEdges(c, a, b.LinearRing(0))
 	}
 
@@ -291,6 +319,7 @@ func onPolygonToPolygon(c DistanceCalculator, a Polygon, b Polygon) bool {
 	aFirstPoint := a.LinearRing(0).Vertex(0)
 	bFirstPoint := b.LinearRing(0).Vertex(0)
 
+	// MinDistance:
 	// If there is at least one point on the the exterior ring of B that is outside the exterior ring
 	// of A, then we have one of these two cases:
 	// * The exterior rings of A and B intersect. The distance can always be found by comparing
@@ -304,7 +333,14 @@ func onPolygonToPolygon(c DistanceCalculator, a Polygon, b Polygon) bool {
 	// that is outside the exterior ring of B.
 	//
 	// As such, we only need to compare the exterior rings if we detect this.
-	if !c.PointInLinearRing(bFirstPoint, a.LinearRing(0)) && !c.PointInLinearRing(aFirstPoint, b.LinearRing(0)) {
+	//
+	// MaxDistance:
+	// The furthest distance between two polygons is always against the exterior rings of each other.
+	// This closely follows the reasoning pointed out in "onPointToPolygon". Holes are always located
+	// inside the exterior ring of a polygon, so the exterior ring will always contain a point
+	// with a larger max distance.
+	if c.DistanceUpdater().IsMaxDistance() ||
+		(!c.PointInLinearRing(bFirstPoint, a.LinearRing(0)) && !c.PointInLinearRing(aFirstPoint, b.LinearRing(0))) {
 		return onShapeEdgesToShapeEdges(c, a.LinearRing(0), b.LinearRing(0))
 	}
 

pkg/geo/geogfn/distance.go

@@ -276,6 +276,11 @@ func (u *spheroidMinDistanceUpdater) OnIntersects() bool {
 	return true
 }
 
+// IsMaxDistance implements the geodist.DistanceUpdater interface.
+func (u *spheroidMinDistanceUpdater) IsMaxDistance() bool {
+	return false
+}
+
 // spheroidDistanceCalculator implements geodist.DistanceCalculator
 type spheroidDistanceCalculator struct {
 	updater *spheroidMinDistanceUpdater

pkg/geo/geomfn/distance.go

@@ -28,6 +28,15 @@ func MinDistance(a *geo.Geometry, b *geo.Geometry) (float64, error) {
 	return minDistanceInternal(a, b, 0)
 }
 
+// MaxDistance returns the maximum distance across every pair of points comprising
+// geometries A and B.
+func MaxDistance(a *geo.Geometry, b *geo.Geometry) (float64, error) {
+	if a.SRID() != b.SRID() {
+		return 0, geo.NewMismatchingSRIDsError(a, b)
+	}
+	return maxDistanceInternal(a, b, math.MaxFloat64)
+}
+
 // DWithin determines if any part of geometry A is within D units of geometry B.
 func DWithin(a *geo.Geometry, b *geo.Geometry, d float64) (bool, error) {
 	if a.SRID() != b.SRID() {
@@ -43,6 +52,31 @@ func DWithin(a *geo.Geometry, b *geo.Geometry, d float64) (bool, error) {
 	return dist <= d, nil
 }
 
+// DFullyWithin determines whether the maximum distance across every pair of points
+// comprising geometries A and B is within D units.
+func DFullyWithin(a *geo.Geometry, b *geo.Geometry, d float64) (bool, error) {
+	if a.SRID() != b.SRID() {
+		return false, geo.NewMismatchingSRIDsError(a, b)
+	}
+	if d < 0 {
+		return false, errors.Newf("dwithin distance cannot be less than zero")
+	}
+	dist, err := maxDistanceInternal(a, b, d)
+	if err != nil {
+		return false, err
+	}
+	return dist <= d, nil
+}
+
+// maxDistanceInternal finds the maximum distance between two geometries.
+// We can re-use the same algorithm as min-distance, allowing skips of checks that involve
+// the interiors or intersections as those will always be less then the maximum min-distance.
+func maxDistanceInternal(a *geo.Geometry, b *geo.Geometry, stopAfterGT float64) (float64, error) {
+	u := newGeomMaxDistanceUpdater(stopAfterGT)
+	c := &geomDistanceCalculator{updater: u}
+	return distanceInternal(a, b, c)
+}
+
 // minDistanceInternal finds the minimum distance between two geometries.
 // This implementation is done in-house, as compared to using GEOS.
 func minDistanceInternal(a *geo.Geometry, b *geo.Geometry, stopAfterLE float64) (float64, error) {
@@ -273,9 +307,60 @@ func (u *geomMinDistanceUpdater) OnIntersects() bool {
 	return true
 }
 
+// IsMaxDistance implements the geodist.DistanceUpdater interface.
+func (u *geomMinDistanceUpdater) IsMaxDistance() bool {
+	return false
+}
+
+// geomMaxDistanceUpdater finds the maximum distance using geom calculations.
+// Methods will return early if it finds a distance > stopAfterGT.
+type geomMaxDistanceUpdater struct {
+	currentValue float64
+	stopAfterGT  float64
+}
+
+var _ geodist.DistanceUpdater = (*geomMaxDistanceUpdater)(nil)
+
+// newGeomMaxDistanceUpdater returns a new geomMaxDistanceUpdater with the
+// correct arguments set up.
+func newGeomMaxDistanceUpdater(stopAfterGT float64) *geomMaxDistanceUpdater {
+	return &geomMaxDistanceUpdater{
+		currentValue: 0,
+		stopAfterGT:  stopAfterGT,
+	}
+}
+
+// Distance implements the DistanceUpdater interface.
+func (u *geomMaxDistanceUpdater) Distance() float64 {
+	return u.currentValue
+}
+
+// Update implements the geodist.DistanceUpdater interface.
+func (u *geomMaxDistanceUpdater) Update(aInterface geodist.Point, bInterface geodist.Point) bool {
+	a := aInterface.(*geomGeodistPoint).Coord
+	b := bInterface.(*geomGeodistPoint).Coord
+
+	dist := coordNorm(coordSub(a, b))
+	if dist > u.currentValue {
+		u.currentValue = dist
+		return dist > u.stopAfterGT
+	}
+	return false
+}
+
+// OnIntersects implements the geodist.DistanceUpdater interface.
+func (u *geomMaxDistanceUpdater) OnIntersects() bool {
+	return false
+}
+
+// IsMaxDistance implements the geodist.DistanceUpdater interface.
+func (u *geomMaxDistanceUpdater) IsMaxDistance() bool {
+	return true
+}
+
 // geomDistanceCalculator implements geodist.DistanceCalculator
 type geomDistanceCalculator struct {
-	updater *geomMinDistanceUpdater
+	updater geodist.DistanceUpdater
 }
 
 var _ geodist.DistanceCalculator = (*geomDistanceCalculator)(nil)