From 2c79710618e7b6d8b7ecff5dc17cfe7abbd59434 Mon Sep 17 00:00:00 2001 From: Karl Wright Date: Mon, 25 Sep 2017 02:20:16 -0400 Subject: [PATCH] LUCENE-7970: Add an exact version of a circle, which uses Vicenty logic to construct out of planes. --- .../spatial3d/geom/GeoCircleFactory.java | 16 + .../lucene/spatial3d/geom/GeoExactCircle.java | 403 ++++++++++++++++++ .../spatial3d/geom/StandardObjects.java | 1 + .../lucene/spatial3d/geom/GeoCircleTest.java | 28 ++ 4 files changed, 448 insertions(+) create mode 100644 lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoExactCircle.java diff --git a/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoCircleFactory.java b/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoCircleFactory.java index 292790f4a0a..7b56e587134 100644 --- a/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoCircleFactory.java +++ b/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoCircleFactory.java @@ -40,4 +40,20 @@ public class GeoCircleFactory { return new GeoStandardCircle(planetModel, latitude, longitude, radius); } + /** Create an exact GeoCircle given specified bounds and desired accuracy. + * @param planetModel is the planet model. + * @param latitude is the center latitude. + * @param longitude is the center longitude. + * @param radius is the radius angle. + * @param accuracy is the maximum difference between the circle approximation and the real circle, as computed using + * the Vincenty formula. + * @return a GeoCircle corresponding to what was specified. + */ + public static GeoCircle makeExactGeoCircle(final PlanetModel planetModel, final double latitude, final double longitude, final double radius, final double accuracy) { + if (radius < Vector.MINIMUM_ANGULAR_RESOLUTION) { + return new GeoDegeneratePoint(planetModel, latitude, longitude); + } + return new GeoExactCircle(planetModel, latitude, longitude, radius, accuracy); + } + } diff --git a/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoExactCircle.java b/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoExactCircle.java new file mode 100644 index 00000000000..7d96cfa0ffe --- /dev/null +++ b/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoExactCircle.java @@ -0,0 +1,403 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.lucene.spatial3d.geom; + +import java.util.Map; +import java.util.List; +import java.util.HashMap; +import java.util.ArrayList; + +import java.io.IOException; +import java.io.InputStream; +import java.io.OutputStream; + +/** + * Circular area with a center and radius. + * + * @lucene.experimental + */ +class GeoExactCircle extends GeoBaseCircle { + /** Center of circle */ + protected final GeoPoint center; + /** Cutoff angle of circle (not quite the same thing as radius) */ + protected final double cutoffAngle; + /** Actual accuracy */ + protected final double actualAccuracy; + /** Planes describing the circle */ + protected final List circlePlanes; + /** Bounds for the planes */ + protected final Map eitherBounds; + /** A point that is on the world and on the circle plane */ + protected final GeoPoint[] edgePoints; + /** The set of notable points for each edge */ + protected final List notableEdgePoints; + /** Notable points for a circle -- there aren't any */ + protected static final GeoPoint[] circlePoints = new GeoPoint[0]; + + /** Constructor. + *@param planetModel is the planet model. + *@param lat is the center latitude. + *@param lon is the center longitude. + *@param cutoffAngle is the cutoff angle for the circle. + *@param accuracy is the allowed error value. + */ + public GeoExactCircle(final PlanetModel planetModel, final double lat, final double lon, final double cutoffAngle, final double accuracy) { + super(planetModel); + if (lat < -Math.PI * 0.5 || lat > Math.PI * 0.5) + throw new IllegalArgumentException("Latitude out of bounds"); + if (lon < -Math.PI || lon > Math.PI) + throw new IllegalArgumentException("Longitude out of bounds"); + if (cutoffAngle < 0.0 || cutoffAngle > Math.PI) + throw new IllegalArgumentException("Cutoff angle out of bounds"); + if (cutoffAngle < Vector.MINIMUM_RESOLUTION) + throw new IllegalArgumentException("Cutoff angle cannot be effectively zero"); + + this.center = new GeoPoint(planetModel, lat, lon); + this.cutoffAngle = cutoffAngle; + + if (accuracy < Vector.MINIMUM_RESOLUTION) { + actualAccuracy = Vector.MINIMUM_RESOLUTION; + } else { + actualAccuracy = accuracy; + } + + if (Math.abs(cutoffAngle - Math.PI) < Vector.MINIMUM_RESOLUTION) { + // Circle is the whole world + this.circlePlanes = null; + this.eitherBounds = null; + this.edgePoints = new GeoPoint[0]; + this.notableEdgePoints = null; + } else if (planetModel.c == planetModel.ab) { + // Sphere + this.eitherBounds = null; + this.circlePlanes = new ArrayList<>(1); + // Compute two points on the circle, with the right angle from the center. We'll use these + // to obtain the perpendicular plane to the circle. + double upperLat = lat + cutoffAngle; + double upperLon = lon; + if (upperLat > Math.PI * 0.5) { + upperLon += Math.PI; + if (upperLon > Math.PI) + upperLon -= 2.0 * Math.PI; + upperLat = Math.PI - upperLat; + } + double lowerLat = lat - cutoffAngle; + double lowerLon = lon; + if (lowerLat < -Math.PI * 0.5) { + lowerLon += Math.PI; + if (lowerLon > Math.PI) + lowerLon -= 2.0 * Math.PI; + lowerLat = -Math.PI - lowerLat; + } + final GeoPoint upperPoint = new GeoPoint(planetModel, upperLat, upperLon); + final GeoPoint lowerPoint = new GeoPoint(planetModel, lowerLat, lowerLon); + // Construct normal plane + final Plane normalPlane = Plane.constructNormalizedZPlane(upperPoint, lowerPoint, center); + // Construct a sided plane that goes through the two points and whose normal is in the normalPlane. + final SidedPlane circlePlane = SidedPlane.constructNormalizedPerpendicularSidedPlane(center, normalPlane, upperPoint, lowerPoint); + if (circlePlane == null) + throw new IllegalArgumentException("Couldn't construct circle plane, probably too small? Cutoff angle = "+cutoffAngle+"; upperPoint = "+upperPoint+"; lowerPoint = "+lowerPoint); + circlePlanes.add(circlePlane); + final GeoPoint recomputedIntersectionPoint = circlePlane.getSampleIntersectionPoint(planetModel, normalPlane); + if (recomputedIntersectionPoint == null) + throw new IllegalArgumentException("Couldn't construct intersection point, probably circle too small? Plane = "+circlePlane); + this.edgePoints = new GeoPoint[]{recomputedIntersectionPoint}; + this.notableEdgePoints = null; + } else { + this.circlePlanes = new ArrayList<>(); + this.notableEdgePoints = new ArrayList<>(); + this.eitherBounds = new HashMap<>(); + // We construct approximation planes until we have a low enough error estimate + final List slices = new ArrayList<>(100); + // Construct four cardinal points, and then we'll build the first two planes + final GeoPoint northPoint = planetModel.surfacePointOnBearing(center, cutoffAngle, Math.PI * 0.5); + final GeoPoint southPoint = planetModel.surfacePointOnBearing(center, cutoffAngle, Math.PI * 1.5); + final GeoPoint eastPoint = planetModel.surfacePointOnBearing(center, cutoffAngle, 0.0); + final GeoPoint westPoint = planetModel.surfacePointOnBearing(center, cutoffAngle, Math.PI); + + this.edgePoints = new GeoPoint[]{northPoint}; + + if (planetModel.c > planetModel.ab) { + // z can be greater than x or y, so ellipse is longer in height than width + slices.add(new ApproximationSlice(center, eastPoint, 0.0, westPoint, Math.PI, northPoint, Math.PI * 0.5)); + slices.add(new ApproximationSlice(center, westPoint, Math.PI, eastPoint, 0.0, southPoint, Math.PI * 1.5)); + } else { + // z will be less than x or y, so ellipse is shorter than it is tall + slices.add(new ApproximationSlice(center, northPoint, Math.PI * 0.5, southPoint, Math.PI * 1.5, eastPoint, 0.0)); + slices.add(new ApproximationSlice(center, southPoint, Math.PI * 1.5, northPoint, Math.PI * 0.5, westPoint, Math.PI)); + } + + // Now, iterate over slices until we have converted all of them into safe SidedPlanes. + while (slices.size() > 0) { + // Peel off a slice from the back + final ApproximationSlice thisSlice = slices.remove(slices.size()-1); + // Assess it to see if it is OK as it is, or needs to be split. + // To do this, we need to look at the part of the circle that will have the greatest error. + // We will need to compute bearing points for these. + final double interpPoint1Bearing = (thisSlice.point1Bearing + thisSlice.middlePointBearing) * 0.5; + final GeoPoint interpPoint1 = planetModel.surfacePointOnBearing(center, cutoffAngle, interpPoint1Bearing); + final double interpPoint2Bearing = (thisSlice.point2Bearing + thisSlice.middlePointBearing) * 0.5; + final GeoPoint interpPoint2 = planetModel.surfacePointOnBearing(center, cutoffAngle, interpPoint2Bearing); + // Is this point on the plane? (that is, is the approximation good enough?) + if (Math.abs(thisSlice.plane.evaluate(interpPoint1)) < actualAccuracy && Math.abs(thisSlice.plane.evaluate(interpPoint2)) < actualAccuracy) { + // Good enough; add it to the list of planes + circlePlanes.add(thisSlice.plane); + notableEdgePoints.add(new GeoPoint[]{thisSlice.endPoint1, thisSlice.endPoint2}); + } else { + // Split the plane into two, and add it back to the end + slices.add(new ApproximationSlice(center, + thisSlice.endPoint1, thisSlice.point1Bearing, + thisSlice.middlePoint, thisSlice.middlePointBearing, + interpPoint1, interpPoint1Bearing)); + slices.add(new ApproximationSlice(center, + thisSlice.middlePoint, thisSlice.middlePointBearing, + thisSlice.endPoint2, thisSlice.point2Bearing, + interpPoint2, interpPoint2Bearing)); + } + } + + //System.out.println("Number of planes needed: "+circlePlanes.size()); + + // Compute bounds + for (int i = 0; i < circlePlanes.size(); i++) { + final SidedPlane thisPlane = circlePlanes.get(i); + final SidedPlane previousPlane = (i == 0)?circlePlanes.get(circlePlanes.size()-1):circlePlanes.get(i-1); + final SidedPlane nextPlane = (i == circlePlanes.size()-1)?circlePlanes.get(0):circlePlanes.get(i+1); + eitherBounds.put(thisPlane, new EitherBound(previousPlane, nextPlane)); + } + } + } + + + /** + * Constructor for deserialization. + * @param planetModel is the planet model. + * @param inputStream is the input stream. + */ + public GeoExactCircle(final PlanetModel planetModel, final InputStream inputStream) throws IOException { + this(planetModel, + SerializableObject.readDouble(inputStream), + SerializableObject.readDouble(inputStream), + SerializableObject.readDouble(inputStream), + SerializableObject.readDouble(inputStream)); + } + + @Override + public void write(final OutputStream outputStream) throws IOException { + SerializableObject.writeDouble(outputStream, center.getLatitude()); + SerializableObject.writeDouble(outputStream, center.getLongitude()); + SerializableObject.writeDouble(outputStream, cutoffAngle); + SerializableObject.writeDouble(outputStream, actualAccuracy); + } + + @Override + public double getRadius() { + return cutoffAngle; + } + + @Override + public GeoPoint getCenter() { + return center; + } + + @Override + protected double distance(final DistanceStyle distanceStyle, final double x, final double y, final double z) { + return distanceStyle.computeDistance(this.center, x, y, z); + } + + @Override + protected void distanceBounds(final Bounds bounds, final DistanceStyle distanceStyle, final double distanceValue) { + // TBD: Compute actual bounds based on distance + getBounds(bounds); + } + + @Override + protected double outsideDistance(final DistanceStyle distanceStyle, final double x, final double y, final double z) { + if (circlePlanes == null) { + return 0.0; + } + if (circlePlanes.size() == 1) { + return distanceStyle.computeDistance(planetModel, circlePlanes.get(0), x, y, z); + } + double outsideDistance = Double.POSITIVE_INFINITY; + for (final SidedPlane plane : circlePlanes) { + final double distance = distanceStyle.computeDistance(planetModel, plane, x, y, z, eitherBounds.get(plane)); + if (distance < outsideDistance) { + outsideDistance = distance; + } + } + return outsideDistance; + } + + @Override + public boolean isWithin(final double x, final double y, final double z) { + if (circlePlanes == null) { + return true; + } + for (final Membership plane : circlePlanes) { + if (!plane.isWithin(x, y, z)) { + return false; + } + } + return true; + } + + @Override + public GeoPoint[] getEdgePoints() { + return edgePoints; + } + + @Override + public boolean intersects(final Plane p, final GeoPoint[] notablePoints, final Membership... bounds) { + if (circlePlanes == null) { + return false; + } + if (circlePlanes.size() == 1) { + return circlePlanes.get(0).intersects(planetModel, p, notablePoints, circlePoints, bounds); + } + for (int edgeIndex = 0; edgeIndex < circlePlanes.size(); edgeIndex++) { + final SidedPlane edge = circlePlanes.get(edgeIndex); + final GeoPoint[] points = notableEdgePoints.get(edgeIndex); + if (edge.intersects(planetModel, p, notablePoints, points, bounds, eitherBounds.get(edge))) { + return true; + } + } + return false; + } + + @Override + public boolean intersects(GeoShape geoShape) { + if (circlePlanes == null) { + return false; + } + if (circlePlanes.size() == 1) { + return geoShape.intersects(circlePlanes.get(0), circlePoints); + } + for (int edgeIndex = 0; edgeIndex < circlePlanes.size(); edgeIndex++) { + final SidedPlane edge = circlePlanes.get(edgeIndex); + final GeoPoint[] points = notableEdgePoints.get(edgeIndex); + if (geoShape.intersects(edge, points, eitherBounds.get(edge))) { + return true; + } + } + return false; + } + + + @Override + public void getBounds(Bounds bounds) { + super.getBounds(bounds); + if (circlePlanes == null) { + return; + } + bounds.addPoint(center); + if (circlePlanes.size() == 1) { + bounds.addPlane(planetModel, circlePlanes.get(0)); + return; + } + // Add bounds for all circle planes + for (final SidedPlane plane : circlePlanes) { + bounds.addPlane(planetModel, plane, eitherBounds.get(plane)); + // We don't bother to compute the intersection bounds since, unless the planet model is pathological, we expect planes to be intersecting at shallow + // angles. + } + } + + @Override + public boolean equals(Object o) { + if (!(o instanceof GeoExactCircle)) + return false; + GeoExactCircle other = (GeoExactCircle) o; + return super.equals(other) && other.center.equals(center) && other.cutoffAngle == cutoffAngle && other.actualAccuracy == actualAccuracy; + } + + @Override + public int hashCode() { + int result = super.hashCode(); + result = 31 * result + center.hashCode(); + long temp = Double.doubleToLongBits(cutoffAngle); + result = 31 * result + (int) (temp ^ (temp >>> 32)); + temp = Double.doubleToLongBits(actualAccuracy); + result = 31 * result + (int) (temp ^ (temp >>> 32)); + return result; + } + + @Override + public String toString() { + return "GeoExactCircle: {planetmodel=" + planetModel+", center=" + center + ", radius=" + cutoffAngle + "(" + cutoffAngle * 180.0 / Math.PI + "), accuracy=" + actualAccuracy + "}"; + } + + /** A membership implementation representing edges that must apply. + */ + protected static class EitherBound implements Membership { + + protected final SidedPlane sideBound1; + protected final SidedPlane sideBound2; + + /** Constructor. + * @param sideBound1 is the first side bound. + * @param sideBound2 is the second side bound. + */ + public EitherBound(final SidedPlane sideBound1, final SidedPlane sideBound2) { + this.sideBound1 = sideBound1; + this.sideBound2 = sideBound2; + } + + @Override + public boolean isWithin(final Vector v) { + return sideBound1.isWithin(v) && sideBound2.isWithin(v); + } + + @Override + public boolean isWithin(final double x, final double y, final double z) { + return sideBound1.isWithin(x,y,z) && sideBound2.isWithin(x,y,z); + } + + @Override + public String toString() { + return "(" + sideBound1 + "," + sideBound2 + ")"; + } + } + + /** A temporary description of a section of circle. + */ + protected static class ApproximationSlice { + public final SidedPlane plane; + public final GeoPoint endPoint1; + public final double point1Bearing; + public final GeoPoint endPoint2; + public final double point2Bearing; + public final GeoPoint middlePoint; + public final double middlePointBearing; + + public ApproximationSlice(final GeoPoint center, + final GeoPoint endPoint1, final double point1Bearing, + final GeoPoint endPoint2, final double point2Bearing, + final GeoPoint middlePoint, final double middlePointBearing) { + this.endPoint1 = endPoint1; + this.point1Bearing = point1Bearing; + this.endPoint2 = endPoint2; + this.point2Bearing = point2Bearing; + this.middlePoint = middlePoint; + this.middlePointBearing = middlePointBearing; + // Construct the plane going through the three given points + this.plane = SidedPlane.constructNormalizedThreePointSidedPlane(center, endPoint1, endPoint2, middlePoint); + } + + } + +} diff --git a/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/StandardObjects.java b/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/StandardObjects.java index bad4651ddbc..bcebf43262c 100644 --- a/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/StandardObjects.java +++ b/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/StandardObjects.java @@ -74,6 +74,7 @@ class StandardObjects { classRegsitry.put(StandardXYZSolid.class, 34); classRegsitry.put(PlanetModel.class, 35); classRegsitry.put(GeoDegeneratePath.class, 36); + classRegsitry.put(GeoExactCircle.class, 37); for (Class clazz : classRegsitry.keySet()){ codeRegsitry.put(classRegsitry.get(clazz), clazz); diff --git a/lucene/spatial3d/src/test/org/apache/lucene/spatial3d/geom/GeoCircleTest.java b/lucene/spatial3d/src/test/org/apache/lucene/spatial3d/geom/GeoCircleTest.java index 84d1b77c6c8..21f8d6aacd4 100755 --- a/lucene/spatial3d/src/test/org/apache/lucene/spatial3d/geom/GeoCircleTest.java +++ b/lucene/spatial3d/src/test/org/apache/lucene/spatial3d/geom/GeoCircleTest.java @@ -22,6 +22,34 @@ import org.junit.Test; public class GeoCircleTest extends LuceneTestCase { + @Test + public void testExactCircle() { + GeoCircle c; + GeoPoint gp; + + // Construct a variety of circles to see how many actual planes are involved + c = new GeoExactCircle(PlanetModel.WGS84, 0.0, 0.0, 0.1, 1e-6); + gp = new GeoPoint(PlanetModel.WGS84, 0.0, 0.2); + assertTrue(!c.isWithin(gp)); + gp = new GeoPoint(PlanetModel.WGS84, 0.0, 0.0); + assertTrue(c.isWithin(gp)); + + c = new GeoExactCircle(PlanetModel.WGS84, 0.1, 0.0, 0.1, 1e-6); + + c = new GeoExactCircle(PlanetModel.WGS84, 0.2, 0.0, 0.1, 1e-6); + + c = new GeoExactCircle(PlanetModel.WGS84, 0.3, 0.0, 0.1, 1e-6); + + c = new GeoExactCircle(PlanetModel.WGS84, 0.4, 0.0, 0.1, 1e-6); + + c = new GeoExactCircle(PlanetModel.WGS84, Math.PI * 0.5, 0.0, 0.1, 1e-6); + gp = new GeoPoint(PlanetModel.WGS84, Math.PI * 0.5 - 0.2, 0.0); + assertTrue(!c.isWithin(gp)); + gp = new GeoPoint(PlanetModel.WGS84, Math.PI * 0.5, 0.0); + assertTrue(c.isWithin(gp)); + + } + @Test public void testCircleDistance() { GeoCircle c;