Skip to content

Commit

Permalink
Fix centroids for small/degenerate polygons
Browse files Browse the repository at this point in the history 
In Esri v2.2.3, the centroid calculate includes centroid fixes for
degenerate (near dimension 1) polygons (Esri/geometry-api-java#227).
This replaces our implemenation (which the original Esri implemenation
was based on) with the fixed Esri version.

Fixes prestodb#10629
  • Loading branch information
@jagill
jagill committed Sep 3, 2019
1 parent 1e1579f commit a0fead3
Show file tree
Hide file tree
Showing 3 changed files with 31 additions and 143 deletions.
142 changes: 1 addition & 141 deletions presto-geospatial/src/main/java/com/facebook/presto/plugin/geospatial/GeoFunctions.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,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.facebook.presto.geospatial.GeometryType;
Expand Down Expand Up @@ -384,33 +383,7 @@ public static Slice stCentroid(@SqlType(GEOMETRY_TYPE_NAME) Slice input)
return serialize(createFromEsriGeometry(new Point(), geometry.getEsriSpatialReference()));
}

Point centroid;
try {
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 PrestoException(INVALID_FUNCTION_ARGUMENT, "Unexpected geometry type: " + geometryType);
}
}
catch (RuntimeException e) {
if (e instanceof PrestoException && ((PrestoException) e).getErrorCode() == INVALID_FUNCTION_ARGUMENT.toErrorCode()) {
throw e;
}
throw new PrestoException(INVALID_FUNCTION_ARGUMENT, format("Cannot compute centroid: %s Use ST_IsValid to confirm that input geometry is valid or compute centroid for a bounding box using ST_Envelope.", e.getMessage()), e);
}
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 @@ -1360,119 +1333,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);
}

final 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 ...geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestBingTileFunctions.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -382,7 +382,7 @@ public void testBingTileZoomLevel()
public void testBingTilePolygon()
{
assertFunction("ST_AsText(bing_tile_polygon(bing_tile('123030123010121')))", VARCHAR, "POLYGON ((59.996337890625 30.11662158281937, 60.00732421875 30.11662158281937, 60.00732421875 30.12612436422458, 59.996337890625 30.12612436422458, 59.996337890625 30.11662158281937))");
assertFunction("ST_AsText(ST_Centroid(bing_tile_polygon(bing_tile('123030123010121'))))", VARCHAR, "POINT (60.0018310442288 30.121372968273892)");
assertFunction("ST_AsText(ST_Centroid(bing_tile_polygon(bing_tile('123030123010121'))))", VARCHAR, "POINT (60.0018310546875 30.121372973521975)");

// Check bottom right corner of a stack of tiles at different zoom levels
assertFunction("ST_AsText(apply(bing_tile_polygon(bing_tile(1, 1, 1)), g -> ST_Point(ST_XMax(g), ST_YMin(g))))", VARCHAR, "POINT (180 -85.05112877980659)");
Expand Down
30 changes: 29 additions & 1 deletion presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestGeoFunctions.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,9 +14,11 @@
package com.facebook.presto.plugin.geospatial;

import com.esri.core.geometry.Point;
import com.esri.core.geometry.ogc.OGCGeometry;
import com.esri.core.geometry.ogc.OGCPoint;
import com.facebook.presto.geospatial.KdbTreeUtils;
import com.facebook.presto.geospatial.Rectangle;
import com.facebook.presto.geospatial.serde.GeometrySerde;
import com.facebook.presto.operator.scalar.AbstractTestFunctions;
import com.facebook.presto.spi.block.Block;
import com.facebook.presto.spi.block.BlockBuilder;
Expand All @@ -30,6 +32,7 @@
import java.util.stream.Collectors;

import static com.facebook.presto.geospatial.KdbTree.buildKdbTree;
import static com.facebook.presto.plugin.geospatial.GeoFunctions.stCentroid;
import static com.facebook.presto.plugin.geospatial.GeometryType.GEOMETRY;
import static com.facebook.presto.spi.StandardErrorCode.INVALID_FUNCTION_ARGUMENT;
import static com.facebook.presto.spi.type.BigintType.BIGINT;
Expand Down Expand Up @@ -219,14 +222,39 @@ public void testSTCentroid()
assertCentroid("POLYGON ((0 0, 0 5, 5 5, 5 0, 0 0), (1 1, 1 2, 2 2, 2 1, 1 1))", new Point(2.5416666666666665, 2.5416666666666665));

// invalid geometry
assertInvalidFunction("ST_Centroid(ST_GeometryFromText('MULTIPOLYGON (((4.903234300000006 52.08474289999999, 4.903234265193165 52.084742934806826, 4.903234299999999 52.08474289999999, 4.903234300000006 52.08474289999999)))'))", "Cannot compute centroid: .* Use ST_IsValid to confirm that input geometry is valid or compute centroid for a bounding box using ST_Envelope.");
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 assertCentroid(String wkt, Point centroid)
{
assertFunction(format("ST_AsText(ST_Centroid(ST_GeometryFromText('%s')))", wkt), VARCHAR, new OGCPoint(centroid, null).asText());
}

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
public void testSTConvexHull()
{
Expand Down

0 comments on commit a0fead3

Please sign in to comment.