Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Improve numerical stability for polygon centroids #13148

Closed
wants to merge 1 commit into from
Closed

Improve numerical stability for polygon centroids #13148

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 11 additions & 0 deletions presto-geospatial/src/main/java/com/facebook/presto/plugin/geospatial/GeoFunctions.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@
import com.esri.core.geometry.ogc.OGCPoint;
import com.esri.core.geometry.ogc.OGCPolygon;
import com.facebook.presto.geospatial.GeometryType;
import com.facebook.presto.geospatial.GeometryUtils;
import com.facebook.presto.geospatial.KdbTree;
import com.facebook.presto.geospatial.Rectangle;
import com.facebook.presto.geospatial.serde.GeometrySerde;
Expand Down Expand Up @@ -104,6 +105,7 @@
import static java.lang.Double.isInfinite;
import static java.lang.Double.isNaN;
import static java.lang.Math.PI;
import static java.lang.Math.abs;
import static java.lang.Math.atan2;
import static java.lang.Math.cos;
import static java.lang.Math.sin;
Expand Down Expand Up @@ -1453,6 +1455,15 @@ private static Point getPolygonSansHolesCentroid(Polygon polygon)
ySum += (current.getY() + next.getY()) * ladder;
signedArea += ladder / 2;
}
if (abs(signedArea) < Double.MIN_NORMAL) {
// We can't trust the calculation, so default to bounding box.
// We've already handled the empty case, so the Envelope is non-null.
OGCGeometry ogcGeom = OGCGeometry.createFromEsriGeometry(polygon, null);
Envelope bounds = GeometryUtils.getEnvelope(ogcGeom);
return new Point(
(bounds.getXMax() + bounds.getXMin()) / 2,
(bounds.getYMax() + bounds.getYMin()) / 2);
}
return new Point(xSum / (signedArea * 6), ySum / (signedArea * 6));
}

Expand Down
53 changes: 38 additions & 15 deletions presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestGeoFunctions.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,14 +14,17 @@
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;
import com.facebook.presto.spi.type.ArrayType;
import com.google.common.collect.ImmutableList;
import io.airlift.slice.Slice;
import org.testng.annotations.BeforeClass;
import org.testng.annotations.Test;

Expand All @@ -30,6 +33,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 All @@ -44,6 +48,17 @@
public class TestGeoFunctions
extends AbstractTestFunctions
{
private static String makeKdbTreeJson()
{
ImmutableList.Builder<Rectangle> rectangles = ImmutableList.builder();
for (double x = 0; x < 10; x += 1) {
for (double y = 0; y < 5; y += 1) {
rectangles.add(new Rectangle(x, y, x + 1, y + 2));
}
}
return KdbTreeUtils.toJson(buildKdbTree(10, new Rectangle(0, 0, 9, 4), rectangles.build()));
}

@BeforeClass
protected void registerFunctions()
{
Expand Down Expand Up @@ -87,17 +102,6 @@ public void testSpatialPartitions()
assertSpatialPartitions(kdbTreeJson, "MULTIPOINT (2 6, 3 7)", 1.2, ImmutableList.of());
}

private static String makeKdbTreeJson()
{
ImmutableList.Builder<Rectangle> rectangles = ImmutableList.builder();
for (double x = 0; x < 10; x += 1) {
for (double y = 0; y < 5; y += 1) {
rectangles.add(new Rectangle(x, y, x + 1, y + 2));
}
}
return KdbTreeUtils.toJson(buildKdbTree(10, new Rectangle(0, 0, 9, 4), rectangles.build()));
}

private void assertSpatialPartitions(String kdbTreeJson, String wkt, List<Integer> expectedPartitions)
{
assertFunction(format("spatial_partitions(cast('%s' as KdbTree), ST_GeometryFromText('%s'))", kdbTreeJson, wkt), new ArrayType(INTEGER), expectedPartitions);
Expand Down Expand Up @@ -212,11 +216,30 @@ public void testSTCentroid()
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.");
}

@Test
public void testSTCentroidNumericStability()
{
String wkt1 = "MULTIPOLYGON (((153.492818 -28.13729, 153.492821 -28.137291, 153.492816 -28.137289, 153.492818 -28.13729)))";
assertApproximateCentroid(wkt1, new Point(153.49282, -28.13729), 1e-5);

String wkt2 = "MULTIPOLYGON (((153.112475 -28.360526, 153.1124759 -28.360527, 153.1124759 -28.360526, 153.112475 -28.360526)))";
assertApproximateCentroid(wkt2, new Point(153.112475, -28.360526), 1e-5);
}

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 centroid, double epsilon)
{
OGCGeometry geometry = OGCGeometry.fromText(wkt);
Slice inputSlice = GeometrySerde.serialize(geometry);
OGCPoint calculatedCentroid = (OGCPoint) GeometrySerde.deserialize(stCentroid(inputSlice));
assertEquals(calculatedCentroid.X(), centroid.getX(), epsilon);
assertEquals(calculatedCentroid.Y(), centroid.getY(), epsilon);
}

@Test
public void testSTConvexHull()
{
Expand Down Expand Up @@ -1084,8 +1107,8 @@ private void assertMultiPoint(String expectedWkt, String... pointWkts)
format(
"ST_MultiPoint(array[%s])",
Arrays.stream(pointWkts)
.map(wkt -> wkt == null ? "null" : format("ST_GeometryFromText('%s')", wkt))
.collect(Collectors.joining(","))),
.map(wkt -> wkt == null ? "null" : format("ST_GeometryFromText('%s')", wkt))
.collect(Collectors.joining(","))),
GEOMETRY,
expectedWkt);
}
Expand All @@ -1096,8 +1119,8 @@ private void assertInvalidMultiPoint(String errorMessage, String... pointWkts)
format(
"ST_MultiPoint(array[%s])",
Arrays.stream(pointWkts)
.map(wkt -> wkt == null ? "null" : format("ST_GeometryFromText('%s')", wkt))
.collect(Collectors.joining(","))),
.map(wkt -> wkt == null ? "null" : format("ST_GeometryFromText('%s')", wkt))
.collect(Collectors.joining(","))),
INVALID_FUNCTION_ARGUMENT,
format("Invalid input to ST_MultiPoint: %s", errorMessage));
}
Expand Down