Skip to content

Commit

Permalink
fix ST_Centroid and ST_Buffer for tiny geometries
Browse files Browse the repository at this point in the history 
  backporting from https://github.com/prestodb/presto/pull/13323/files
  • Loading branch information
@cshao239 @wendigo
cshao239 authored and wendigo committed Jan 2, 2024
1 parent 71499d5 commit 4155e9f
Show file tree
Hide file tree
Showing 4 changed files with 48 additions and 135 deletions.
3 changes: 2 additions & 1 deletion docs/src/main/sphinx/functions/geospatial.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -174,7 +174,8 @@ Returns the closure of the combinatorial boundary of this geometry.

:::{function} ST_Buffer(Geometry, distance) -> Geometry
Returns the geometry that represents all points whose distance from the specified geometry
is less than or equal to the specified distance.
is less than or equal to the specified distance. If the points of the geometry are extremely
close together (``delta < 1e-8``), this might return an empty geometry.
:::

:::{function} ST_Difference(Geometry, Geometry) -> Geometry
Expand Down
134 changes: 1 addition & 133 deletions plugin/trino-geospatial/src/main/java/io/trino/plugin/geospatial/GeoFunctions.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,6 @@
import com.esri.core.geometry.ogc.OGCGeometry;
import com.esri.core.geometry.ogc.OGCGeometryCollection;
import com.esri.core.geometry.ogc.OGCLineString;
import com.esri.core.geometry.ogc.OGCMultiPolygon;
import com.esri.core.geometry.ogc.OGCPoint;
import com.esri.core.geometry.ogc.OGCPolygon;
import com.google.common.base.Joiner;
Expand Down Expand Up @@ -433,25 +432,7 @@ public static Slice stCentroid(@SqlType(GEOMETRY_TYPE_NAME) Slice input)
return serialize(createFromEsriGeometry(new Point(), geometry.getEsriSpatialReference()));
}

Point centroid;
switch (geometryType) {
case MULTI_POINT:
centroid = computePointsCentroid((MultiVertexGeometry) geometry.getEsriGeometry());
break;
case LINE_STRING:
case MULTI_LINE_STRING:
centroid = computeLineCentroid((Polyline) geometry.getEsriGeometry());
break;
case POLYGON:
centroid = computePolygonCentroid((Polygon) geometry.getEsriGeometry());
break;
case MULTI_POLYGON:
centroid = computeMultiPolygonCentroid((OGCMultiPolygon) geometry);
break;
default:
throw new TrinoException(INVALID_FUNCTION_ARGUMENT, "Unexpected geometry type: " + geometryType);
}
return serialize(createFromEsriGeometry(centroid, geometry.getEsriSpatialReference()));
return serialize(geometry.centroid());
}

@Description("Returns the minimum convex geometry that encloses all input geometries")
Expand Down Expand Up @@ -1609,119 +1590,6 @@ private static void verifySameSpatialReference(OGCGeometry leftGeometry, OGCGeom
checkArgument(Objects.equals(leftGeometry.getEsriSpatialReference(), rightGeometry.getEsriSpatialReference()), "Input geometries must have the same spatial reference");
}

// Points centroid is arithmetic mean of the input points
private static Point computePointsCentroid(MultiVertexGeometry multiVertex)
{
double xSum = 0;
double ySum = 0;
for (int i = 0; i < multiVertex.getPointCount(); i++) {
Point point = multiVertex.getPoint(i);
xSum += point.getX();
ySum += point.getY();
}
return new Point(xSum / multiVertex.getPointCount(), ySum / multiVertex.getPointCount());
}

// Lines centroid is weighted mean of each line segment, weight in terms of line length
private static Point computeLineCentroid(Polyline polyline)
{
double xSum = 0;
double ySum = 0;
double weightSum = 0;
for (int i = 0; i < polyline.getPathCount(); i++) {
Point startPoint = polyline.getPoint(polyline.getPathStart(i));
Point endPoint = polyline.getPoint(polyline.getPathEnd(i) - 1);
double dx = endPoint.getX() - startPoint.getX();
double dy = endPoint.getY() - startPoint.getY();
double length = sqrt(dx * dx + dy * dy);
weightSum += length;
xSum += (startPoint.getX() + endPoint.getX()) * length / 2;
ySum += (startPoint.getY() + endPoint.getY()) * length / 2;
}
return new Point(xSum / weightSum, ySum / weightSum);
}

// Polygon centroid: area weighted average of centroids in case of holes
private static Point computePolygonCentroid(Polygon polygon)
{
int pathCount = polygon.getPathCount();

if (pathCount == 1) {
return getPolygonSansHolesCentroid(polygon);
}

double xSum = 0;
double ySum = 0;
double areaSum = 0;

for (int i = 0; i < pathCount; i++) {
int startIndex = polygon.getPathStart(i);
int endIndex = polygon.getPathEnd(i);

Polygon sansHoles = getSubPolygon(polygon, startIndex, endIndex);

Point centroid = getPolygonSansHolesCentroid(sansHoles);
double area = sansHoles.calculateArea2D();

xSum += centroid.getX() * area;
ySum += centroid.getY() * area;
areaSum += area;
}

return new Point(xSum / areaSum, ySum / areaSum);
}

private static Polygon getSubPolygon(Polygon polygon, int startIndex, int endIndex)
{
Polyline boundary = new Polyline();
boundary.startPath(polygon.getPoint(startIndex));
for (int i = startIndex + 1; i < endIndex; i++) {
Point current = polygon.getPoint(i);
boundary.lineTo(current);
}

Polygon newPolygon = new Polygon();
newPolygon.add(boundary, false);
return newPolygon;
}

// Polygon sans holes centroid:
// c[x] = (Sigma(x[i] + x[i + 1]) * (x[i] * y[i + 1] - x[i + 1] * y[i]), for i = 0 to N - 1) / (6 * signedArea)
// c[y] = (Sigma(y[i] + y[i + 1]) * (x[i] * y[i + 1] - x[i + 1] * y[i]), for i = 0 to N - 1) / (6 * signedArea)
private static Point getPolygonSansHolesCentroid(Polygon polygon)
{
int pointCount = polygon.getPointCount();
double xSum = 0;
double ySum = 0;
double signedArea = 0;
for (int i = 0; i < pointCount; i++) {
Point current = polygon.getPoint(i);
Point next = polygon.getPoint((i + 1) % polygon.getPointCount());
double ladder = current.getX() * next.getY() - next.getX() * current.getY();
xSum += (current.getX() + next.getX()) * ladder;
ySum += (current.getY() + next.getY()) * ladder;
signedArea += ladder / 2;
}
return new Point(xSum / (signedArea * 6), ySum / (signedArea * 6));
}

// MultiPolygon centroid is weighted mean of each polygon, weight in terms of polygon area
private static Point computeMultiPolygonCentroid(OGCMultiPolygon multiPolygon)
{
double xSum = 0;
double ySum = 0;
double weightSum = 0;
for (int i = 0; i < multiPolygon.numGeometries(); i++) {
Point centroid = computePolygonCentroid((Polygon) multiPolygon.geometryN(i).getEsriGeometry());
Polygon polygon = (Polygon) multiPolygon.geometryN(i).getEsriGeometry();
double weight = polygon.calculateArea2D();
weightSum += weight;
xSum += centroid.getX() * weight;
ySum += centroid.getY() * weight;
}
return new Point(xSum / weightSum, ySum / weightSum);
}

private static boolean envelopes(Slice left, Slice right, EnvelopesPredicate predicate)
{
Envelope leftEnvelope = deserializeEnvelope(left);
Expand Down
2 changes: 1 addition & 1 deletion plugin/trino-geospatial/src/test/java/io/trino/plugin/geospatial/TestBingTileFunctions.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -454,7 +454,7 @@ public void testBingTilePolygon()

assertThat(assertions.function("ST_AsText", "ST_Centroid(bing_tile_polygon(bing_tile('123030123010121')))"))
.hasType(VARCHAR)
.isEqualTo("POINT (60.0018310442288 30.121372968273892)");
.isEqualTo("POINT (60.0018310546875 30.121372973521975)");

// Check bottom right corner of a stack of tiles at different zoom levels
assertThat(assertions.function("ST_AsText", "apply(bing_tile_polygon(bing_tile(1, 1, 1)), g -> ST_Point(ST_XMax(g), ST_YMin(g)))"))
Expand Down
44 changes: 44 additions & 0 deletions plugin/trino-geospatial/src/test/java/io/trino/plugin/geospatial/TestGeoFunctions.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,10 +14,12 @@
package io.trino.plugin.geospatial;

import com.esri.core.geometry.Point;
import com.esri.core.geometry.ogc.OGCGeometry;
import com.esri.core.geometry.ogc.OGCPoint;
import com.google.common.collect.ImmutableList;
import io.trino.geospatial.KdbTreeUtils;
import io.trino.geospatial.Rectangle;
import io.trino.geospatial.serde.GeometrySerde;
import io.trino.spi.block.Block;
import io.trino.spi.block.BlockBuilder;
import io.trino.spi.type.ArrayType;
Expand All @@ -35,13 +37,15 @@

import static com.google.common.collect.ImmutableList.toImmutableList;
import static io.trino.geospatial.KdbTree.buildKdbTree;
import static io.trino.plugin.geospatial.GeoFunctions.stCentroid;
import static io.trino.plugin.geospatial.GeometryType.GEOMETRY;
import static io.trino.spi.type.BooleanType.BOOLEAN;
import static io.trino.spi.type.DoubleType.DOUBLE;
import static io.trino.spi.type.IntegerType.INTEGER;
import static io.trino.spi.type.VarcharType.VARCHAR;
import static io.trino.testing.assertions.TrinoExceptionAssert.assertTrinoExceptionThrownBy;
import static org.assertj.core.api.Assertions.assertThat;
import static org.junit.jupiter.api.Assertions.assertEquals;
import static org.junit.jupiter.api.TestInstance.Lifecycle.PER_CLASS;
import static org.junit.jupiter.api.parallel.ExecutionMode.CONCURRENT;

Expand Down Expand Up @@ -259,6 +263,19 @@ public void testSTBuffer()

assertTrinoExceptionThrownBy(assertions.function("ST_Buffer", "ST_Point(0, 0)", "nan()")::evaluate)
.hasMessage("distance is NaN");

// For small polygons, there was a bug in ESRI that throw an NPE. This
// was fixed (https://github.com/Esri/geometry-api-java/pull/243) to
// return an empty geometry instead. Ideally, these would return
// something approximately like `ST_Buffer(ST_Centroid(geometry))`.
assertThat(assertions.function("ST_IsEmpty", "ST_Buffer(ST_Buffer(ST_Point(177.50102959662, 64.726807421691), 0.0000000001), 0.00005)"))
.hasType(BOOLEAN)
.isEqualTo(true);

assertThat(assertions.function("ST_IsEmpty", "ST_Buffer(ST_GeometryFromText(" +
"'POLYGON ((177.0 64.0, 177.0000000001 64.0, 177.0000000001 64.0000000001, 177.0 64.0000000001, 177.0 64.0))'), 0.01)"))
.hasType(BOOLEAN)
.isEqualTo(true);
}

@Test
Expand Down Expand Up @@ -299,6 +316,33 @@ public void testSTCentroid()
assertThat(assertions.function("ST_AsText", "ST_Centroid(ST_GeometryFromText('POLYGON ((0 0, 0 5, 5 5, 5 0, 0 0), (1 1, 1 2, 2 2, 2 1, 1 1))'))"))
.hasType(VARCHAR)
.isEqualTo("POINT (2.5416666666666665 2.5416666666666665)");

assertApproximateCentroid("MULTIPOLYGON (((4.903234300000006 52.08474289999999, 4.903234265193165 52.084742934806826, 4.903234299999999 52.08474289999999, 4.903234300000006 52.08474289999999)))", new Point(4.9032343, 52.0847429), 1e-7);

// Numerical stability tests
assertApproximateCentroid(
"MULTIPOLYGON (((153.492818 -28.13729, 153.492821 -28.137291, 153.492816 -28.137289, 153.492818 -28.13729)))",
new Point(153.49282, -28.13729), 1e-5);
assertApproximateCentroid(
"MULTIPOLYGON (((153.112475 -28.360526, 153.1124759 -28.360527, 153.1124759 -28.360526, 153.112475 -28.360526)))",
new Point(153.112475, -28.360526), 1e-5);
assertApproximateCentroid(
"POLYGON ((4.903234300000006 52.08474289999999, 4.903234265193165 52.084742934806826, 4.903234299999999 52.08474289999999, 4.903234300000006 52.08474289999999))",
new Point(4.9032343, 52.0847429), 1e-6);
assertApproximateCentroid(
"MULTIPOLYGON (((4.903234300000006 52.08474289999999, 4.903234265193165 52.084742934806826, 4.903234299999999 52.08474289999999, 4.903234300000006 52.08474289999999)))",
new Point(4.9032343, 52.0847429), 1e-6);
assertApproximateCentroid(
"POLYGON ((-81.0387349 29.20822, -81.039974 29.210597, -81.0410331 29.2101579, -81.0404758 29.2090879, -81.0404618 29.2090609, -81.040433 29.209005, -81.0404269 29.208993, -81.0404161 29.2089729, -81.0398001 29.20779, -81.0387349 29.20822), (-81.0404229 29.208986, -81.04042 29.2089809, -81.0404269 29.208993, -81.0404229 29.208986))",
new Point(-81.039885, 29.209191), 1e-6);
}

private void assertApproximateCentroid(String wkt, Point expectedCentroid, double epsilon)
{
OGCPoint actualCentroid = (OGCPoint) GeometrySerde.deserialize(
stCentroid(GeometrySerde.serialize(OGCGeometry.fromText(wkt))));
assertEquals(actualCentroid.X(), expectedCentroid.getX(), epsilon);
assertEquals(actualCentroid.Y(), expectedCentroid.getY(), epsilon);
}

@Test
Expand Down

0 comments on commit 4155e9f

Please sign in to comment.