geo/geomfn: implement ST_MaxDistance / ST_DFullyWithin

MaxDistance/DFullyWithin is a specialisation of distance, where it takes
the maximum of distances found instead and disregards intersections and
closest points. The `stopAfterLE` condition becomes `stopAfterGT`.

Release note (sql change): Implemented the ST_MaxDistance and
ST_DFullyWithin function for geometries.

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 all of geometry_a is in distance units of geometry_b. In other words, the max distance between geometry_a and geometry_b is less than 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>
@@ -863,6 +865,8 @@ has no relationship with the commit order of concurrent transactions.</p>
 </span></td></tr>
 <tr><td><a name="st_linestringfromwkb"></a><code>st_linestringfromwkb(wkb: <a href="bytes.html">bytes</a>, srid: <a href="int.html">int</a>) &rarr; geometry</code></td><td><span class="funcdesc"><p>Returns the Geometry from a WKB representation with an SRID. If the shape underneath is not LineString, NULL is returned.</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 between 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.
@@ -159,11 +161,13 @@ func onPointToEdgesExceptFirstEdgeStart(c DistanceCalculator, a Point, b shapeWi
 		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
+		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
+				}
 			}
 		}
 	}
@@ -185,7 +189,8 @@ func onPointToLineString(c DistanceCalculator, a Point, b LineString) bool {
 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
 	// the exterior ring.
-	if !c.PointInLinearRing(a, b.LinearRing(0)) {
+	// Also, if we are just calculating max distance, we only want the distance to the exterior.
+	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 +215,11 @@ 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()
+			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 +237,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
+						}
 					}
 				}
 			}
@@ -252,7 +261,8 @@ func onLineStringToPolygon(c DistanceCalculator, a LineString, b Polygon) bool {
 	// 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)) {
+	// If we are looking for max distance, we only need to check against the outer ring anyway.
+	if c.DistanceUpdater().IsMaxDistance() || !c.PointInLinearRing(a.Vertex(0), b.LinearRing(0)) {
 		return onShapeEdgesToShapeEdges(c, a, b.LinearRing(0))
 	}
 
@@ -304,7 +314,10 @@ 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)) {
+	//
+	// If we are only looking at the max distance, we only want to compare exteriors.
+	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 between geometries A and B.
+// This is determined by the maximum minimum-distance between the two geometries.
+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 if any part of geometry A is fully within D units of geometry B.
+// This is determined by the maximum minimum-distance between the two geometries.
+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)