diff --git a/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/Geo3DPoint.java b/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/Geo3DPoint.java index 45e17b72f3b..f0f202049db 100644 --- a/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/Geo3DPoint.java +++ b/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/Geo3DPoint.java @@ -22,6 +22,7 @@ import java.util.ArrayList; import org.apache.lucene.document.Field; import org.apache.lucene.document.FieldType; import org.apache.lucene.index.PointValues; +import org.apache.lucene.spatial3d.geom.Vector; import org.apache.lucene.spatial3d.geom.GeoPoint; import org.apache.lucene.spatial3d.geom.GeoShape; import org.apache.lucene.spatial3d.geom.PlanetModel; @@ -149,11 +150,8 @@ public final class Geo3DPoint extends Field { checkLongitude(longitude); polyPoints.add(new GeoPoint(PlanetModel.WGS84, fromDegrees(latitude), fromDegrees(longitude))); } - // We don't know what the sense of the polygon is without providing the index of one vertex we know to be convex. - // Since that doesn't fit with the "super-simple API" requirements, we just use the index of the first one, and people have to just - // know to do it that way. - final int convexPointIndex = 0; - final GeoShape shape = GeoPolygonFactory.makeGeoPolygon(PlanetModel.WGS84, polyPoints, convexPointIndex); + // We use the polygon constructor that looks at point order. + final GeoShape shape = GeoPolygonFactory.makeGeoPolygon(PlanetModel.WGS84, polyPoints, null); return newShapeQuery(field, shape); } diff --git a/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoConcavePolygon.java b/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoConcavePolygon.java new file mode 100644 index 00000000000..3250a80b6ec --- /dev/null +++ b/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoConcavePolygon.java @@ -0,0 +1,393 @@ +/* + * 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.ArrayList; +import java.util.BitSet; +import java.util.List; +import java.util.HashMap; +import java.util.Map; + +/** + * GeoConcavePolygon objects are generic building blocks of more complex structures. + * The only restrictions on these objects are: (1) they must be concave; (2) they must have + * a maximum extent larger than PI. Violating either one of these limits will + * cause the logic to fail. + * + * @lucene.experimental + */ +public class GeoConcavePolygon extends GeoBasePolygon { + /** The list of polygon points */ + protected final List points; + /** A bitset describing, for each edge, whether it is internal or not */ + protected final BitSet isInternalEdges; + /** The list of holes. If a point is in the hole, it is *not* in the polygon */ + protected final List holes; + + /** A list of edges */ + protected SidedPlane[] edges = null; + /** A list of inverted edges */ + protected SidedPlane[] invertedEdges = null; + /** The set of notable points for each edge */ + protected GeoPoint[][] notableEdgePoints = null; + /** A point which is on the boundary of the polygon */ + protected GeoPoint[] edgePoints = null; + /** Tracking the maximum distance we go at any one time, so to be sure it's legal */ + protected double fullDistance = 0.0; + /** Set to true when the polygon is complete */ + protected boolean isDone = false; + /** A bounds object for each sided plane */ + protected Map eitherBounds = null; + + /** + * Create a concave polygon from a list of points. The first point must be on the + * external edge. + *@param planetModel is the planet model. + *@param pointList is the list of points to create the polygon from. + */ + public GeoConcavePolygon(final PlanetModel planetModel, final List pointList) { + this(planetModel, pointList, null); + } + + /** + * Create a concave polygon from a list of points. The first point must be on the + * external edge. + *@param planetModel is the planet model. + *@param pointList is the list of points to create the polygon from. + *@param holes is the list of GeoPolygon objects that describe holes in the concave polygon. Null == no holes. + */ + public GeoConcavePolygon(final PlanetModel planetModel, final List pointList, final List holes) { + super(planetModel); + this.points = pointList; + this.holes = holes; + this.isInternalEdges = new BitSet(); + done(false); + } + + /** + * Create a concave polygon from a list of points, keeping track of which boundaries + * are internal. This is used when creating a polygon as a building block for another shape. + *@param planetModel is the planet model. + *@param pointList is the set of points to create the polygon from. + *@param internalEdgeFlags is a bitset describing whether each edge is internal or not. + *@param returnEdgeInternal is true when the final return edge is an internal one. + */ + public GeoConcavePolygon(final PlanetModel planetModel, + final List pointList, + final BitSet internalEdgeFlags, + final boolean returnEdgeInternal) { + this(planetModel, pointList, null, internalEdgeFlags, returnEdgeInternal); + } + + /** + * Create a concave polygon from a list of points, keeping track of which boundaries + * are internal. This is used when creating a polygon as a building block for another shape. + *@param planetModel is the planet model. + *@param pointList is the set of points to create the polygon from. + *@param holes is the list of GeoPolygon objects that describe holes in the concave polygon. Null == no holes. + *@param internalEdgeFlags is a bitset describing whether each edge is internal or not. + *@param returnEdgeInternal is true when the final return edge is an internal one. + */ + public GeoConcavePolygon(final PlanetModel planetModel, + final List pointList, + final List holes, + final BitSet internalEdgeFlags, + final boolean returnEdgeInternal) { + super(planetModel); + this.points = pointList; + this.holes = holes; + this.isInternalEdges = internalEdgeFlags; + done(returnEdgeInternal); + } + + /** + * Create a concave polygon, with a starting latitude and longitude. + * Accepts only values in the following ranges: lat: {@code -PI/2 -> PI/2}, lon: {@code -PI -> PI} + *@param planetModel is the planet model. + *@param startLatitude is the latitude of the first point. + *@param startLongitude is the longitude of the first point. + */ + public GeoConcavePolygon(final PlanetModel planetModel, + final double startLatitude, + final double startLongitude) { + this(planetModel, startLatitude, startLongitude, null); + } + + /** + * Create a concave polygon, with a starting latitude and longitude. + * Accepts only values in the following ranges: lat: {@code -PI/2 -> PI/2}, lon: {@code -PI -> PI} + *@param planetModel is the planet model. + *@param startLatitude is the latitude of the first point. + *@param startLongitude is the longitude of the first point. + *@param holes is the list of GeoPolygon objects that describe holes in the concave polygon. Null == no holes. + */ + public GeoConcavePolygon(final PlanetModel planetModel, + final double startLatitude, + final double startLongitude, + final List holes) { + super(planetModel); + points = new ArrayList<>(); + this.holes = holes; + isInternalEdges = new BitSet(); + points.add(new GeoPoint(planetModel, startLatitude, startLongitude)); + } + + /** + * Add a point to the polygon. + * Accepts only values in the following ranges: lat: {@code -PI/2 -> PI/2}, lon: {@code -PI -> PI} + * + * @param latitude is the latitude of the next point. + * @param longitude is the longitude of the next point. + * @param isInternalEdge is true if the edge just added with this point should be considered "internal", and not + * intersected as part of the intersects() operation. + */ + public void addPoint(final double latitude, final double longitude, final boolean isInternalEdge) { + if (isDone) + throw new IllegalStateException("Can't call addPoint() if done() already called"); + if (isInternalEdge) + isInternalEdges.set(points.size() - 1); + points.add(new GeoPoint(planetModel, latitude, longitude)); + } + + /** + * Finish the polygon, by connecting the last added point with the starting point. + *@param isInternalReturnEdge is true if the return edge (back to start) is an internal one. + */ + public void done(final boolean isInternalReturnEdge) { + if (isDone) + throw new IllegalStateException("Can't call done() more than once"); + // If fewer than 3 points, can't do it. + if (points.size() < 3) + throw new IllegalArgumentException("Polygon needs at least three points."); + + if (isInternalReturnEdge) + isInternalEdges.set(points.size() - 1); + + isDone = true; + + // Time to construct the planes. If the polygon is truly concave then any adjacent point + // to a segment can provide an exterior measurement. Note: We build the true planes + // here and use the logic to return what *isn't* inside all of them. + edges = new SidedPlane[points.size()]; + invertedEdges = new SidedPlane[points.size()]; + notableEdgePoints = new GeoPoint[points.size()][]; + + for (int i = 0; i < points.size(); i++) { + final GeoPoint start = points.get(i); + final GeoPoint end = points.get(legalIndex(i + 1)); + final double distance = start.arcDistance(end); + if (distance > fullDistance) + fullDistance = distance; + final GeoPoint check = points.get(legalIndex(i + 2)); + // Here note the flip of the sense of the sided plane!! + final SidedPlane sp = new SidedPlane(check, false, start, end); + //System.out.println("Created edge "+sp+" using start="+start+" end="+end+" check="+check); + edges[i] = sp; + invertedEdges[i] = new SidedPlane(sp); + notableEdgePoints[i] = new GeoPoint[]{start, end}; + } + // In order to naively confirm that the polygon is concave, I would need to + // check every edge, and verify that every point (other than the edge endpoints) + // is within the edge's sided plane. This is an order n^2 operation. That's still + // not wrong, though, because everything else about polygons has a similar cost. + for (int edgeIndex = 0; edgeIndex < edges.length; edgeIndex++) { + final SidedPlane edge = edges[edgeIndex]; + for (int pointIndex = 0; pointIndex < points.size(); pointIndex++) { + if (pointIndex != edgeIndex && pointIndex != legalIndex(edgeIndex + 1)) { + if (edge.isWithin(points.get(pointIndex))) + throw new IllegalArgumentException("Polygon is not concave: Point " + points.get(pointIndex) + " Edge " + edge); + } + } + } + + // For each edge, create a bounds object. + eitherBounds = new HashMap<>(edges.length); + for (int edgeIndex = 0; edgeIndex < edges.length; edgeIndex++) { + eitherBounds.put(edges[edgeIndex], new EitherBound(invertedEdges[edgeIndex])); + } + + // Pick an edge point arbitrarily + edgePoints = new GeoPoint[]{points.get(0)}; + } + + /** Compute a legal point index from a possibly illegal one, that may have wrapped. + *@param index is the index. + *@return the normalized index. + */ + protected int legalIndex(int index) { + while (index >= points.size()) + index -= points.size(); + return index; + } + + @Override + public boolean isWithin(final double x, final double y, final double z) { + // If present within *any* plane, then it is a member, except where there are holes. + boolean isMember = false; + for (final SidedPlane edge : edges) { + if (edge.isWithin(x, y, z)) { + isMember = true; + break; + } + } + if (isMember == false) { + return false; + } + if (holes != null) { + for (final GeoPolygon polygon : holes) { + if (polygon.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) { + // The bounding planes are inverted and complementary. For intersection computation, we + // cannot use them as bounds. They are independent hemispheres. + for (int edgeIndex = 0; edgeIndex < edges.length; edgeIndex++) { + final SidedPlane edge = edges[edgeIndex]; + final GeoPoint[] points = this.notableEdgePoints[edgeIndex]; + if (!isInternalEdges.get(edgeIndex)) { + if (edge.intersects(planetModel, p, notablePoints, points, bounds, eitherBounds.get(edge))) { + //System.err.println(" intersects!"); + return true; + } + } + } + if (holes != null) { + // Each hole needs to be looked at for intersection too, since a shape can be entirely within the hole + for (final GeoPolygon hole : holes) { + if (hole.intersects(p, notablePoints, bounds)) { + return true; + } + } + } + //System.err.println(" no intersection"); + return false; + } + + /** A membership implementation representing polygon edges that all must apply. + */ + protected class EitherBound implements Membership { + + protected final SidedPlane exception; + + /** Constructor. + */ + public EitherBound(final SidedPlane exception) { + this.exception = exception; + } + + @Override + public boolean isWithin(final Vector v) { + for (final SidedPlane edge : invertedEdges) { + if (edge != exception && !edge.isWithin(v)) { + return false; + } + } + return true; + } + + @Override + public boolean isWithin(final double x, final double y, final double z) { + for (final SidedPlane edge : invertedEdges) { + if (edge != exception && !edge.isWithin(x, y, z)) { + return false; + } + } + return true; + } + } + + @Override + public void getBounds(Bounds bounds) { + super.getBounds(bounds); + bounds.isWide(); + + // Add all the points + for (final GeoPoint point : points) { + bounds.addPoint(point); + } + + // Add planes with membership. + for (final SidedPlane edge : edges) { + bounds.addPlane(planetModel, edge, eitherBounds.get(edge)); + } + } + + @Override + protected double outsideDistance(final DistanceStyle distanceStyle, final double x, final double y, final double z) { + double minimumDistance = Double.MAX_VALUE; + for (final GeoPoint edgePoint : points) { + final double newDist = distanceStyle.computeDistance(edgePoint, x,y,z); + if (newDist < minimumDistance) { + minimumDistance = newDist; + } + } + for (final SidedPlane edgePlane : edges) { + final double newDist = distanceStyle.computeDistance(planetModel, edgePlane, x, y, z, eitherBounds.get(edgePlane)); + if (newDist < minimumDistance) { + minimumDistance = newDist; + } + } + return minimumDistance; + } + + @Override + public boolean equals(Object o) { + if (!(o instanceof GeoConcavePolygon)) + return false; + GeoConcavePolygon other = (GeoConcavePolygon) o; + if (!super.equals(other)) + return false; + if (!other.isInternalEdges.equals(isInternalEdges)) + return false; + if (other.holes != null || holes != null) { + if (other.holes == null || holes == null) { + return false; + } + if (!other.holes.equals(holes)) { + return false; + } + } + return (other.points.equals(points)); + } + + @Override + public int hashCode() { + int result = super.hashCode(); + result = 31 * result + points.hashCode(); + if (holes != null) { + result = 31 * result + holes.hashCode(); + } + return result; + } + + @Override + public String toString() { + return "GeoConcavePolygon: {planetmodel=" + planetModel + ", points=" + points + ", internalEdges=" + isInternalEdges + ((holes== null)?"":", holes=" + holes) + "}"; + } +} + diff --git a/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoConvexPolygon.java b/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoConvexPolygon.java index fb024b6efee..749971f2f8f 100755 --- a/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoConvexPolygon.java +++ b/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoConvexPolygon.java @@ -19,6 +19,8 @@ package org.apache.lucene.spatial3d.geom; import java.util.ArrayList; import java.util.BitSet; import java.util.List; +import java.util.HashMap; +import java.util.Map; /** * GeoConvexPolygon objects are generic building blocks of more complex structures. @@ -33,6 +35,8 @@ public class GeoConvexPolygon extends GeoBasePolygon { protected final List points; /** A bitset describing, for each edge, whether it is internal or not */ protected final BitSet isInternalEdges; + /** The list of holes. If a point is in the hole, it is *not* in the polygon */ + protected final List holes; /** A list of edges */ protected SidedPlane[] edges = null; @@ -44,6 +48,8 @@ public class GeoConvexPolygon extends GeoBasePolygon { protected double fullDistance = 0.0; /** Set to true when the polygon is complete */ protected boolean isDone = false; + /** A bounds object for each sided plane */ + protected Map eitherBounds = null; /** * Create a convex polygon from a list of points. The first point must be on the @@ -52,8 +58,20 @@ public class GeoConvexPolygon extends GeoBasePolygon { *@param pointList is the list of points to create the polygon from. */ public GeoConvexPolygon(final PlanetModel planetModel, final List pointList) { + this(planetModel, pointList, null); + } + + /** + * Create a convex polygon from a list of points. The first point must be on the + * external edge. + *@param planetModel is the planet model. + *@param pointList is the list of points to create the polygon from. + *@param holes is the list of GeoPolygon objects that describe holes in the complex polygon. Null == no holes. + */ + public GeoConvexPolygon(final PlanetModel planetModel, final List pointList, final List holes) { super(planetModel); this.points = pointList; + this.holes = holes; this.isInternalEdges = new BitSet(); done(false); } @@ -66,10 +84,30 @@ public class GeoConvexPolygon extends GeoBasePolygon { *@param internalEdgeFlags is a bitset describing whether each edge is internal or not. *@param returnEdgeInternal is true when the final return edge is an internal one. */ - public GeoConvexPolygon(final PlanetModel planetModel, final List pointList, final BitSet internalEdgeFlags, - final boolean returnEdgeInternal) { + public GeoConvexPolygon(final PlanetModel planetModel, + final List pointList, + final BitSet internalEdgeFlags, + final boolean returnEdgeInternal) { + this(planetModel, pointList, null, internalEdgeFlags, returnEdgeInternal); + } + + /** + * Create a convex polygon from a list of points, keeping track of which boundaries + * are internal. This is used when creating a polygon as a building block for another shape. + *@param planetModel is the planet model. + *@param pointList is the set of points to create the polygon from. + *@param holes is the list of GeoPolygon objects that describe holes in the complex polygon. Null == no holes. + *@param internalEdgeFlags is a bitset describing whether each edge is internal or not. + *@param returnEdgeInternal is true when the final return edge is an internal one. + */ + public GeoConvexPolygon(final PlanetModel planetModel, + final List pointList, + final List holes, + final BitSet internalEdgeFlags, + final boolean returnEdgeInternal) { super(planetModel); this.points = pointList; + this.holes = holes; this.isInternalEdges = internalEdgeFlags; done(returnEdgeInternal); } @@ -81,9 +119,27 @@ public class GeoConvexPolygon extends GeoBasePolygon { *@param startLatitude is the latitude of the first point. *@param startLongitude is the longitude of the first point. */ - public GeoConvexPolygon(final PlanetModel planetModel, final double startLatitude, final double startLongitude) { + public GeoConvexPolygon(final PlanetModel planetModel, + final double startLatitude, + final double startLongitude) { + this(planetModel, startLatitude, startLongitude, null); + } + + /** + * Create a convex polygon, with a starting latitude and longitude. + * Accepts only values in the following ranges: lat: {@code -PI/2 -> PI/2}, lon: {@code -PI -> PI} + *@param planetModel is the planet model. + *@param startLatitude is the latitude of the first point. + *@param startLongitude is the longitude of the first point. + *@param holes is the list of GeoPolygon objects that describe holes in the complex polygon. Null == no holes. + */ + public GeoConvexPolygon(final PlanetModel planetModel, + final double startLatitude, + final double startLongitude, + final List holes) { super(planetModel); points = new ArrayList<>(); + this.holes = holes; isInternalEdges = new BitSet(); points.add(new GeoPoint(planetModel, startLatitude, startLongitude)); } @@ -138,12 +194,6 @@ public class GeoConvexPolygon extends GeoBasePolygon { edges[i] = sp; notableEdgePoints[i] = new GeoPoint[]{start, end}; } - createCenterPoint(); - } - - /** Compute a reasonable center point. - */ - protected void createCenterPoint() { // In order to naively confirm that the polygon is convex, I would need to // check every edge, and verify that every point (other than the edge endpoints) // is within the edge's sided plane. This is an order n^2 operation. That's still @@ -157,6 +207,14 @@ public class GeoConvexPolygon extends GeoBasePolygon { } } } + + // For each edge, create a bounds object. + eitherBounds = new HashMap<>(edges.length); + for (final SidedPlane edge : edges) { + eitherBounds.put(edge, new EitherBound(edge)); + } + + // Pick an edge point arbitrarily edgePoints = new GeoPoint[]{points.get(0)}; } @@ -176,6 +234,13 @@ public class GeoConvexPolygon extends GeoBasePolygon { if (!edge.isWithin(x, y, z)) return false; } + if (holes != null) { + for (final GeoPolygon polygon : holes) { + if (polygon.isWithin(x, y, z)) { + return false; + } + } + } return true; } @@ -191,26 +256,58 @@ public class GeoConvexPolygon extends GeoBasePolygon { final SidedPlane edge = edges[edgeIndex]; final GeoPoint[] points = this.notableEdgePoints[edgeIndex]; if (!isInternalEdges.get(edgeIndex)) { - //System.err.println(" non-internal edge "+edge); - // Edges flagged as 'internal only' are excluded from the matching - // Construct boundaries - final Membership[] membershipBounds = new Membership[edges.length - 1]; - int count = 0; - for (int otherIndex = 0; otherIndex < edges.length; otherIndex++) { - if (otherIndex != edgeIndex) { - membershipBounds[count++] = edges[otherIndex]; - } - } - if (edge.intersects(planetModel, p, notablePoints, points, bounds, membershipBounds)) { + if (edge.intersects(planetModel, p, notablePoints, points, bounds, eitherBounds.get(edge))) { //System.err.println(" intersects!"); return true; } } } + if (holes != null) { + // Each hole needs to be looked at for intersection too, since a shape can be entirely within the hole + for (final GeoPolygon hole : holes) { + if (hole.intersects(p, notablePoints, bounds)) { + return true; + } + } + } //System.err.println(" no intersection"); return false; } + /** A membership implementation representing polygon edges that all must apply. + */ + protected class EitherBound implements Membership { + + protected final SidedPlane exception; + + /** Constructor. + */ + public EitherBound(final SidedPlane exception) { + this.exception = exception; + } + + @Override + public boolean isWithin(final Vector v) { + for (final SidedPlane edge : edges) { + if (edge != exception && !edge.isWithin(v)) { + return false; + } + } + return true; + } + + @Override + public boolean isWithin(final double x, final double y, final double z) { + for (final SidedPlane edge : edges) { + if (edge != exception && !edge.isWithin(x, y, z)) { + return false; + } + } + return true; + } + } + + @Override public void getBounds(Bounds bounds) { super.getBounds(bounds); @@ -221,17 +318,8 @@ public class GeoConvexPolygon extends GeoBasePolygon { } // Add planes with membership. - for (int edgeIndex = 0; edgeIndex < edges.length; edgeIndex++) { - final SidedPlane edge = edges[edgeIndex]; - // Construct boundaries - final Membership[] membershipBounds = new Membership[edges.length - 1]; - int count = 0; - for (int otherIndex = 0; otherIndex < edges.length; otherIndex++) { - if (otherIndex != edgeIndex) { - membershipBounds[count++] = edges[otherIndex]; - } - } - bounds.addPlane(planetModel, edge, membershipBounds); + for (final SidedPlane edge : edges) { + bounds.addPlane(planetModel, edge, eitherBounds.get(edge)); } } @@ -244,16 +332,8 @@ public class GeoConvexPolygon extends GeoBasePolygon { minimumDistance = newDist; } } - for (int edgeIndex = 0; edgeIndex < edges.length; edgeIndex++) { - final Plane edgePlane = edges[edgeIndex]; - final Membership[] membershipBounds = new Membership[edges.length - 1]; - int count = 0; - for (int otherIndex = 0; otherIndex < edges.length; otherIndex++) { - if (otherIndex != edgeIndex) { - membershipBounds[count++] = edges[otherIndex]; - } - } - final double newDist = distanceStyle.computeDistance(planetModel, edgePlane, x, y, z, membershipBounds); + for (final SidedPlane edgePlane : edges) { + final double newDist = distanceStyle.computeDistance(planetModel, edgePlane, x, y, z, eitherBounds.get(edgePlane)); if (newDist < minimumDistance) { minimumDistance = newDist; } @@ -270,6 +350,14 @@ public class GeoConvexPolygon extends GeoBasePolygon { return false; if (!other.isInternalEdges.equals(isInternalEdges)) return false; + if (other.holes != null || holes != null) { + if (other.holes == null || holes == null) { + return false; + } + if (!other.holes.equals(holes)) { + return false; + } + } return (other.points.equals(points)); } @@ -277,12 +365,15 @@ public class GeoConvexPolygon extends GeoBasePolygon { public int hashCode() { int result = super.hashCode(); result = 31 * result + points.hashCode(); + if (holes != null) { + result = 31 * result + holes.hashCode(); + } return result; } @Override public String toString() { - return "GeoConvexPolygon: {planetmodel=" + planetModel + ", points=" + points + "}"; + return "GeoConvexPolygon: {planetmodel=" + planetModel + ", points=" + points + ", internalEdges="+isInternalEdges+((holes== null)?"":", holes=" + holes) + "}"; } } diff --git a/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoPolygonFactory.java b/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoPolygonFactory.java index 8ee4290df3f..d8040de64e8 100755 --- a/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoPolygonFactory.java +++ b/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoPolygonFactory.java @@ -19,6 +19,12 @@ package org.apache.lucene.spatial3d.geom; import java.util.ArrayList; import java.util.BitSet; import java.util.List; +import java.util.Random; +import java.util.Iterator; +import java.util.Set; +import java.util.HashSet; +import java.util.Map; +import java.util.HashMap; /** * Class which constructs a GeoMembershipShape representing an arbitrary polygon. @@ -37,138 +43,620 @@ public class GeoPolygonFactory { * its neighbors determines inside/outside for the entire polygon. * @return a GeoPolygon corresponding to what was specified. */ - public static GeoPolygon makeGeoPolygon(final PlanetModel planetModel, final List pointList, final int convexPointIndex) { + public static GeoPolygon makeGeoPolygon(final PlanetModel planetModel, + final List pointList, + final int convexPointIndex) { + return makeGeoPolygon(planetModel, pointList, convexPointIndex, null); + } + + /** + * Create a GeoMembershipShape of the right kind given the specified bounds. + * + * @param pointList is a list of the GeoPoints to build an arbitrary polygon out of. + * @param convexPointIndex is the index of a single convex point whose conformation with + * its neighbors determines inside/outside for the entire polygon. + * @param holes is a list of polygons representing "holes" in the outside polygon. Null == none. + * @return a GeoPolygon corresponding to what was specified. + */ + public static GeoPolygon makeGeoPolygon(final PlanetModel planetModel, + final List pointList, + final int convexPointIndex, + final List holes) { // The basic operation uses a set of points, two points determining one particular edge, and a sided plane // describing membership. return buildPolygonShape(planetModel, pointList, convexPointIndex, getLegalIndex(convexPointIndex + 1, pointList.size()), new SidedPlane(pointList.get(getLegalIndex(convexPointIndex - 1, pointList.size())), pointList.get(convexPointIndex), pointList.get(getLegalIndex(convexPointIndex + 1, pointList.size()))), - false); + false, + holes, + null); } - /** Build a GeoMembershipShape given points, starting edge, and whether starting edge is internal or not. + /** Create a GeoPolygon using the specified points and holes, using order to determine + * siding of the polygon. Much like ESRI, this method uses clockwise to indicate the space + * on the same side of the shape as being inside, and counter-clockwise to indicate the + * space on the opposite side as being inside. + * @param pointList is a list of the GeoPoints to build an arbitrary polygon out of. If points go + * clockwise from a given pole, then that pole should be within the polygon. If points go + * counter-clockwise, then that pole should be outside the polygon. + * @return a GeoPolygon corresponding to what was specified. + */ + public static GeoPolygon makeGeoPolygon(final PlanetModel planetModel, + final List pointList) { + return makeGeoPolygon(planetModel, pointList, null); + } + + /** Create a GeoPolygon using the specified points and holes, using order to determine + * siding of the polygon. Much like ESRI, this method uses clockwise to indicate the space + * on the same side of the shape as being inside, and counter-clockwise to indicate the + * space on the opposite side as being inside. + * @param pointList is a list of the GeoPoints to build an arbitrary polygon out of. If points go + * clockwise from a given pole, then that pole should be within the polygon. If points go + * counter-clockwise, then that pole should be outside the polygon. + * @param holes is a list of polygons representing "holes" in the outside polygon. Null == none. + * @return a GeoPolygon corresponding to what was specified. + */ + public static GeoPolygon makeGeoPolygon(final PlanetModel planetModel, + final List pointList, + final List holes) { + // Create a random number generator. Effectively this furnishes us with a repeatable sequence + // of points to use for poles. + final Random generator = new Random(1234); + while (true) { + // Pick the next random pole + final double poleLat = generator.nextDouble() * Math.PI - Math.PI * 0.5; + final double poleLon = generator.nextDouble() * Math.PI * 2.0 - Math.PI; + final GeoPoint pole = new GeoPoint(planetModel, poleLat, poleLon); + // Is it inside or outside? + final Boolean isPoleInside = isInsidePolygon(pole, pointList); + if (isPoleInside != null) { + // Legal pole + return makeGeoPolygon(planetModel, pointList, holes, pole, isPoleInside); + } + // If pole choice was illegal, try another one + } + } + + /** + * Create a GeoPolygon using the specified points and holes and a test point. + * + * @param pointList is a list of the GeoPoints to build an arbitrary polygon out of. + * @param holes is a list of polygons representing "holes" in the outside polygon. Null == none. + * @param testPoint is a test point that is either known to be within the polygon area, or not. + * @param testPointInside is true if the test point is within the area, false otherwise. + * @return a GeoPolygon corresponding to what was specified. + */ + public static GeoPolygon makeGeoPolygon(final PlanetModel planetModel, + final List pointList, + final List holes, + final GeoPoint testPoint, + final boolean testPointInside) { + // We will be trying twice to find the right GeoPolygon, using alternate siding choices for the first polygon + // side. While this looks like it might be 2x as expensive as it could be, there's really no other choice I can + // find. + final SidedPlane initialPlane = new SidedPlane(testPoint, pointList.get(0), pointList.get(1)); + // We don't know if this is the correct siding choice. We will only know as we build the complex polygon. + // So we need to be prepared to try both possibilities. + final GeoPolygon trial = buildPolygonShape(planetModel, pointList, 0, 1, initialPlane, false, holes, testPoint); + if (trial == null) { + // The testPoint was within the shape. Was that intended? + if (testPointInside) { + // Yes: build it for real + return buildPolygonShape(planetModel, pointList, 0, 1, initialPlane, false, holes, null); + } + // No: do the complement and return that. + return buildPolygonShape(planetModel, pointList, 0, 1, new SidedPlane(initialPlane), false, holes, null); + } else { + // The testPoint was outside the shape. Was that intended? + if (!testPointInside) { + // Yes: return what we just built + return trial; + } + // No: return the complement + return buildPolygonShape(planetModel, pointList, 0, 1, new SidedPlane(initialPlane), false, holes, null); + } + } + + /** For a specified point and a list of poly points, determine based on point order whether the + * point should be considered in or out of the polygon. + * @param point is the point to check. + * @param polyPoints is the list of points comprising the polygon. + * @return null if the point is illegal, otherwise false if the point is inside and true if the point is outside + * of the polygon. + */ + protected static Boolean isInsidePolygon(final GeoPoint point, final List polyPoints) { + // First, compute sine and cosine of pole point latitude and longitude + final double norm = 1.0 / point.magnitude(); + final double xyDenom = Math.sqrt(point.x * point.x + point.y * point.y); + final double sinLatitude = point.z * norm; + final double cosLatitude = xyDenom * norm; + final double sinLongitude; + final double cosLongitude; + if (Math.abs(xyDenom) < Vector.MINIMUM_RESOLUTION) { + sinLongitude = 0.0; + cosLongitude = 1.0; + } else { + final double xyNorm = 1.0 / xyDenom; + sinLongitude = point.y * xyNorm; + cosLongitude = point.x * xyNorm; + } + + // Now, compute the incremental arc distance around the points of the polygon + double arcDistance = 0.0; + Double prevAngle = null; + for (final GeoPoint polyPoint : polyPoints) { + final Double angle = computeAngle(polyPoint, sinLatitude, cosLatitude, sinLongitude, cosLongitude); + if (angle == null) { + return null; + } + //System.out.println("Computed angle: "+angle); + if (prevAngle != null) { + // Figure out delta between prevAngle and current angle, and add it to arcDistance + double angleDelta = angle - prevAngle; + if (angleDelta < -Math.PI) { + angleDelta += Math.PI * 2.0; + } + if (angleDelta > Math.PI) { + angleDelta -= Math.PI * 2.0; + } + if (Math.abs(angleDelta - Math.PI) < Vector.MINIMUM_RESOLUTION) { + return null; + } + //System.out.println(" angle delta = "+angleDelta); + arcDistance += angleDelta; + } + prevAngle = angle; + } + if (prevAngle != null) { + final Double lastAngle = computeAngle(polyPoints.get(0), sinLatitude, cosLatitude, sinLongitude, cosLongitude); + if (lastAngle == null) { + return null; + } + //System.out.println("Computed last angle: "+lastAngle); + // Figure out delta and add it + double angleDelta = lastAngle - prevAngle; + if (angleDelta < -Math.PI) { + angleDelta += Math.PI * 2.0; + } + if (angleDelta > Math.PI) { + angleDelta -= Math.PI * 2.0; + } + if (Math.abs(angleDelta - Math.PI) < Vector.MINIMUM_RESOLUTION) { + return null; + } + //System.out.println(" angle delta = "+angleDelta); + arcDistance += angleDelta; + } + // Clockwise == inside == negative + //System.out.println("Arcdistance = "+arcDistance); + if (Math.abs(arcDistance) < Vector.MINIMUM_RESOLUTION) { + // No idea what direction, so try another pole. + return null; + } + return arcDistance > 0.0; + } + + protected static Double computeAngle(final GeoPoint point, + final double sinLatitude, + final double cosLatitude, + final double sinLongitude, + final double cosLongitude) { + // Coordinate rotation formula: + // x1 = x0 cos T - y0 sin T + // y1 = x0 sin T + y0 cos T + // We need to rotate the point in question into the coordinate frame specified by + // the lat and lon trig functions. + // To do this we need to do two rotations on it. First rotation is in x/y. Second rotation is in x/z. + // So: + // x1 = x0 cos az - y0 sin az + // y1 = x0 sin az + y0 cos az + // z1 = z0 + // x2 = x1 cos al - z1 sin al + // y2 = y1 + // z2 = x1 sin al + z1 cos al + + final double x1 = point.x * cosLongitude - point.y * sinLongitude; + final double y1 = point.x * sinLongitude + point.y * cosLongitude; + final double z1 = point.z; + //final double x2 = x1 * cosLatitude - z1 * sinLatitude; + final double y2 = y1; + final double z2 = x1 * sinLatitude + z1 * cosLatitude; + + // Now we should be looking down the X axis; the original point has rotated coordinates (N, 0, 0). + // So we can just compute the angle using y2 and z2. (If Math.sqrt(y2*y2 + z2 * z2) is 0.0, then the point is on the pole and we need another one). + if (Math.sqrt(y2*y2 + z2*z2) < Vector.MINIMUM_RESOLUTION) { + return null; + } + + return Math.atan2(z2, y2); + } + + /** Build a GeoPolygon out of one concave part and multiple convex parts given points, starting edge, and whether starting edge is internal or not. * @param pointsList is a list of the GeoPoints to build an arbitrary polygon out of. - * @param startPointIndex is one of the points constituting the starting edge. - * @param endPointIndex is another of the points constituting the starting edge. + * @param startPointIndex is the first of the points, constituting the starting edge. * @param startingEdge is the plane describing the starting edge. * @param isInternalEdge is true if the specified edge is an internal one. - * @return a GeoMembershipShape corresponding to what was specified. + * @param holes is the list of holes in the polygon, or null if none. + * @param testPoint is an (optional) test point, which will be used to determine if we are generating + * a shape with the proper sidedness. It is passed in only when the test point is supposed to be outside + * of the generated polygon. In this case, if the generated polygon is found to contain the point, the + * method exits early with a null return value. + * This only makes sense in the context of evaluating both possible choices and using logic to determine + * which result to use. If the test point is supposed to be within the shape, then it must be outside of the + * complement shape. If the test point is supposed to be outside the shape, then it must be outside of the + * original shape. Either way, we can figure out the right thing to use. + * @return a GeoMembershipShape corresponding to what was specified, or null if what was specified + * was inconsistent with what we generated. Specifically, if we specify an exterior point that is + * found in the interior of the shape we create here we return null, which is a signal that we chose + * our initial plane sidedness backwards. */ - public static GeoPolygon buildPolygonShape(final PlanetModel planetModel, final List pointsList, final int startPointIndex, final int endPointIndex, final SidedPlane startingEdge, final boolean isInternalEdge) { - // Algorithm as follows: - // Start with sided edge. Go through all points in some order. For each new point, determine if the point is within all edges considered so far. - // If not, put it into a list of points for recursion. If it is within, add new edge and keep going. - // Once we detect a point that is within, if there are points put aside for recursion, then call recursively. + public static GeoPolygon buildPolygonShape( + final PlanetModel planetModel, + final List pointsList, + final int startPointIndex, + final int endPointIndex, + final SidedPlane startingEdge, + final boolean isInternalEdge, + final List holes, + final GeoPoint testPoint) { - // Current composite. This is what we'll actually be returning. + // It could be the case that we need a concave polygon. So we need to try and look for that case + // as part of the general code for constructing complex polygons. + + // Note that there can be only one concave polygon. + + // The code here must keep track of two lists of sided planes. The first list contains the planes consistent with + // a concave polygon. This list will grow and shrink. The second list is built starting at the current edge that + // was last consistent with the concave polygon, and contains all edges consistent with a convex polygon. + // When that sequence of edges is done, then an internal edge is created and the identified points are converted to a + // convex polygon. That internal edge is used to extend the list of edges in the concave polygon edge list. + + // The edge buffer. + final EdgeBuffer edgeBuffer = new EdgeBuffer(pointsList, startPointIndex, endPointIndex, startingEdge, isInternalEdge); + + // Current composite. This is what we'll actually be returning. This will have a number of convex polygons, and + // maybe a single concave one too. final GeoCompositePolygon rval = new GeoCompositePolygon(); - final List recursionList = new ArrayList(); - final List currentList = new ArrayList(); - final BitSet internalEdgeList = new BitSet(); - final List currentPlanes = new ArrayList(); + // Starting state: + // The stopping point + Edge stoppingPoint = edgeBuffer.pickOne(); + // The current edge + Edge currentEdge = stoppingPoint; + + // Progressively look for convex sections. If we find one, we emit it and replace it. + // Keep going until we have been around once and nothing needed to change, and then + // do the concave polygon, if necessary. + while (true) { - // Initialize the current list and current planes - currentList.add(pointsList.get(startPointIndex)); - currentList.add(pointsList.get(endPointIndex)); - internalEdgeList.set(currentPlanes.size(), isInternalEdge); - currentPlanes.add(startingEdge); + if (currentEdge == null) { + // We're done! + break; + } + + // Find convexity around the current edge, if any + final Boolean foundIt = findConvexPolygon(planetModel, currentEdge, rval, edgeBuffer, holes, testPoint); + if (foundIt == null) { + return null; + } + + if (foundIt) { + // New start point + stoppingPoint = edgeBuffer.pickOne(); + currentEdge = stoppingPoint; + // back around + continue; + } + + // Otherwise, go on to the next + currentEdge = edgeBuffer.getNext(currentEdge); + if (currentEdge == stoppingPoint) { + break; + } + } + + // If there's anything left in the edge buffer, convert to concave polygon. + if (makeConcavePolygon(planetModel, rval, edgeBuffer, holes, testPoint) == false) { + return null; + } + + return rval; + } + + /** Look for a concave polygon in the remainder of the edgebuffer. + * By this point, if there are any edges in the edgebuffer, they represent a concave polygon. + * @param planetModel is the planet model. + * @param rval is the composite polygon we're building. + * @param edgeBuffer is the edge buffer. + * @param holes is the optional list of holes. + * @param testPoint is the optional test point. + * @return true unless the testPoint caused failure. + */ + protected static boolean makeConcavePolygon(final PlanetModel planetModel, + final GeoCompositePolygon rval, + final EdgeBuffer edgeBuffer, + final List holes, + final GeoPoint testPoint) { + if (edgeBuffer.size() == 0) { + return true; + } + + // If there are less than three edges, something got messed up somehow. Don't know how this + // can happen but check. + if (edgeBuffer.size() < 3) { + throw new IllegalStateException("Ending edge buffer had only "+edgeBuffer.size()+" edges"); + } + + // Create the list of points + final List points = new ArrayList(edgeBuffer.size()); + final BitSet internalEdges = new BitSet(edgeBuffer.size()-1); - // Now, scan all remaining points, in order. We'll use an index and just add to it. - for (int i = 0; i < pointsList.size() - 2; i++) { - GeoPoint newPoint = pointsList.get(getLegalIndex(i + endPointIndex + 1, pointsList.size())); - if (isWithin(newPoint, currentPlanes)) { - // Construct a sided plane based on the last two points, and the previous point - SidedPlane newBoundary = new SidedPlane(currentList.get(currentList.size() - 2), newPoint, currentList.get(currentList.size() - 1)); - // Construct a sided plane based on the return trip - SidedPlane returnBoundary = new SidedPlane(currentList.get(currentList.size() - 1), currentList.get(0), newPoint); - // Verify that none of the points beyond the new point in the list are inside the polygon we'd - // be creating if we stopped making the current polygon right now. - boolean pointInside = false; - for (int j = i + 1; j < pointsList.size() - 2; j++) { - GeoPoint checkPoint = pointsList.get(getLegalIndex(j + endPointIndex + 1, pointsList.size())); - boolean isInside = true; - if (isInside && !newBoundary.isWithin(checkPoint)) - isInside = false; - if (isInside && !returnBoundary.isWithin(checkPoint)) - isInside = false; - if (isInside) { - for (SidedPlane plane : currentPlanes) { - if (!plane.isWithin(checkPoint)) { - isInside = false; + Edge edge = edgeBuffer.pickOne(); + boolean isInternal = false; + for (int i = 0; i < edgeBuffer.size(); i++) { + points.add(edge.startPoint); + if (i < edgeBuffer.size() - 1) { + internalEdges.set(i, edge.isInternal); + } else { + isInternal = edge.isInternal; + } + edge = edgeBuffer.getNext(edge); + } + + if (testPoint != null && holes != null && holes.size() > 0) { + // No holes, for test + final GeoPolygon testPolygon = new GeoConcavePolygon(planetModel, points, null, internalEdges, isInternal); + if (testPolygon.isWithin(testPoint)) { + return false; + } + } + + final GeoPolygon realPolygon = new GeoConcavePolygon(planetModel, points, holes, internalEdges, isInternal); + if (testPoint != null && (holes == null || holes.size() == 0)) { + if (realPolygon.isWithin(testPoint)) { + return false; + } + } + + rval.addShape(realPolygon); + return true; + } + + /** Look for a convex polygon at the specified edge. If we find it, create one and adjust the edge buffer. + * @param planetModel is the planet model. + * @param currentEdge is the current edge to use starting the search. + * @param rval is the composite polygon to build. + * @param edgeBuffer is the edge buffer. + * @param holes is the optional list of holes. + * @param testPoint is the optional test point. + * @return null if the testPoint is within any polygon detected, otherwise true if a convex polygon was created. + */ + protected static Boolean findConvexPolygon(final PlanetModel planetModel, + final Edge currentEdge, + final GeoCompositePolygon rval, + final EdgeBuffer edgeBuffer, + final List holes, + final GeoPoint testPoint) { + + //System.out.println("Looking at edge "+currentEdge+" with startpoint "+currentEdge.startPoint+" endpoint "+currentEdge.endPoint); + + // Initialize the structure. + // We don't keep track of order here; we just care about membership. + // The only exception is the head and tail pointers. + final Set includedEdges = new HashSet<>(); + includedEdges.add(currentEdge); + Edge firstEdge = currentEdge; + Edge lastEdge = currentEdge; + + // First, walk towards the end until we need to stop + while (true) { + if (firstEdge.startPoint == lastEdge.endPoint) { + break; + } + final Edge newLastEdge = edgeBuffer.getNext(lastEdge); + if (isWithin(newLastEdge.endPoint, includedEdges)) { + //System.out.println(" maybe can extend to next edge"); + // Found a candidate for extension. But do some other checks first. Basically, we need to know if we construct a polygon + // here will overlap with other remaining points? + final SidedPlane returnBoundary; + if (firstEdge.startPoint != newLastEdge.endPoint) { + returnBoundary = new SidedPlane(firstEdge.endPoint, firstEdge.startPoint, newLastEdge.endPoint); + } else { + returnBoundary = null; + } + // The complete set of sided planes for the tentative new polygon include the ones in includedEdges, plus the one from newLastEdge, + // plus the new tentative return boundary. We have to make sure there are no points from elsewhere within the tentative convex polygon. + boolean foundPointInside = false; + final Iterator edgeIterator = edgeBuffer.iterator(); + while (edgeIterator.hasNext()) { + final Edge edge = edgeIterator.next(); + if (!includedEdges.contains(edge) && edge != newLastEdge) { + // This edge has a point to check + if (edge.startPoint != newLastEdge.endPoint) { + // look at edge.startPoint + if (isWithin(edge.startPoint, includedEdges, newLastEdge, returnBoundary)) { + //System.out.println(" nope; point within found: "+edge.startPoint); + foundPointInside = true; + break; + } + } + if (edge.endPoint != firstEdge.startPoint) { + // look at edge.endPoint + if (isWithin(edge.endPoint, includedEdges, newLastEdge, returnBoundary)) { + //System.out.println(" nope; point within found: "+edge.endPoint); + foundPointInside = true; break; } } } - if (isInside) { - pointInside = true; - break; - } } - if (!pointInside) { - // Any excluded points? - boolean isInternalBoundary = recursionList.size() > 0; - if (isInternalBoundary) { - // Handle exclusion - recursionList.add(newPoint); - recursionList.add(currentList.get(currentList.size() - 1)); - if (recursionList.size() == pointsList.size()) { - // We are trying to recurse with a list the same size as the one we started with. - // Clearly, the polygon cannot be constructed - throw new IllegalArgumentException("Polygon is illegal; cannot be decomposed into convex parts"); - } - // We want the other side for the recursion - SidedPlane otherSideNewBoundary = new SidedPlane(newBoundary); - rval.addShape(buildPolygonShape(planetModel, recursionList, recursionList.size() - 2, recursionList.size() - 1, otherSideNewBoundary, true)); - recursionList.clear(); - } - currentList.add(newPoint); - internalEdgeList.set(currentPlanes.size(), isInternalBoundary); - currentPlanes.add(newBoundary); + + if (!foundPointInside) { + //System.out.println(" extending!"); + // Extend the polygon by the new last edge + includedEdges.add(newLastEdge); + lastEdge = newLastEdge; + // continue extending in this direction + continue; + } + } + // We can't extend any more in this direction, so break from the loop. + break; + } + + // Now, walk towards the beginning until we need to stop + while (true) { + if (firstEdge.startPoint == lastEdge.endPoint) { + break; + } + final Edge newFirstEdge = edgeBuffer.getPrevious(firstEdge); + if (isWithin(newFirstEdge.startPoint, includedEdges)) { + //System.out.println(" maybe can extend to previous edge"); + // Found a candidate for extension. But do some other checks first. Basically, we need to know if we construct a polygon + // here will overlap with other remaining points? + final SidedPlane returnBoundary; + if (newFirstEdge.startPoint != lastEdge.endPoint) { + returnBoundary = new SidedPlane(lastEdge.startPoint, lastEdge.endPoint, newFirstEdge.startPoint); } else { - recursionList.add(newPoint); + returnBoundary = null; + } + // The complete set of sided planes for the tentative new polygon include the ones in includedEdges, plus the one from newLastEdge, + // plus the new tentative return boundary. We have to make sure there are no points from elsewhere within the tentative convex polygon. + boolean foundPointInside = false; + final Iterator edgeIterator = edgeBuffer.iterator(); + while (edgeIterator.hasNext()) { + final Edge edge = edgeIterator.next(); + if (!includedEdges.contains(edge) && edge != newFirstEdge) { + // This edge has a point to check + if (edge.startPoint != lastEdge.endPoint) { + // look at edge.startPoint + if (isWithin(edge.startPoint, includedEdges, newFirstEdge, returnBoundary)) { + //System.out.println(" nope; point within found: "+edge.startPoint); + foundPointInside = true; + break; + } + } + if (edge.endPoint != newFirstEdge.startPoint) { + // look at edge.endPoint + if (isWithin(edge.endPoint, includedEdges, newFirstEdge, returnBoundary)) { + //System.out.println(" nope; point within found: "+edge.endPoint); + foundPointInside = true; + break; + } + } + } + } + + if (!foundPointInside) { + //System.out.println(" extending!"); + // Extend the polygon by the new last edge + includedEdges.add(newFirstEdge); + firstEdge = newFirstEdge; + // continue extending in this direction + continue; } - } else { - recursionList.add(newPoint); } + // We can't extend any more in this direction, so break from the loop. + break; } - boolean returnEdgeInternalBoundary = recursionList.size() > 0; - if (returnEdgeInternalBoundary) { - // The last step back to the start point had a recursion, so take care of that before we complete our work - recursionList.add(currentList.get(0)); - recursionList.add(currentList.get(currentList.size() - 1)); - if (recursionList.size() == pointsList.size()) { - // We are trying to recurse with a list the same size as the one we started with. - // Clearly, the polygon cannot be constructed - throw new IllegalArgumentException("Polygon is illegal; cannot be decomposed into convex parts"); - } - // Construct a sided plane based on these two points, and the previous point - SidedPlane newBoundary = new SidedPlane(currentList.get(currentList.size() - 2), currentList.get(0), currentList.get(currentList.size() - 1)); - // We want the other side for the recursion - SidedPlane otherSideNewBoundary = new SidedPlane(newBoundary); - rval.addShape(buildPolygonShape(planetModel, recursionList, recursionList.size() - 2, recursionList.size() - 1, otherSideNewBoundary, true)); - recursionList.clear(); + // Ok, figure out what we've accumulated. If it is enough for a polygon, build it. + if (includedEdges.size() < 2) { + //System.out.println("Done edge "+currentEdge+": no poly found"); + return false; } - // Now, add in the current shape. - rval.addShape(new GeoConvexPolygon(planetModel, currentList, internalEdgeList, returnEdgeInternalBoundary)); - //System.out.println("Done creating polygon"); - return rval; + + // It's enough to build a convex polygon + //System.out.println("Edge "+currentEdge+": Found complex poly"); + + // Create the point list and edge list, starting with the first edge and going to the last. The return edge will be between + // the start point of the first edge and the end point of the last edge. If the first edge start point is the same as the last edge end point, + // it's a degenerate case and we want to just clean out the edge buffer entirely. + + final List points = new ArrayList(includedEdges.size()); + final BitSet internalEdges = new BitSet(includedEdges.size()-1); + final boolean returnIsInternal; + + if (firstEdge.startPoint == lastEdge.endPoint) { + // Degenerate case!! There is no return edge -- or rather, we already have it. + Edge edge = firstEdge; + points.add(edge.startPoint); + int i = 0; + while (true) { + if (edge == lastEdge) { + break; + } + points.add(edge.endPoint); + internalEdges.set(i++, edge.isInternal); + edge = edgeBuffer.getNext(edge); + } + returnIsInternal = lastEdge.isInternal; + edgeBuffer.clear(); + } else { + // Build the return edge (internal, of course) + final SidedPlane returnSidedPlane = new SidedPlane(firstEdge.endPoint, false, firstEdge.startPoint, lastEdge.endPoint); + final Edge returnEdge = new Edge(firstEdge.startPoint, lastEdge.endPoint, returnSidedPlane, true); + + // Build point list and edge list + final List edges = new ArrayList(includedEdges.size()); + returnIsInternal = true; + + Edge edge = firstEdge; + points.add(edge.startPoint); + int i = 0; + while (true) { + points.add(edge.endPoint); + internalEdges.set(i++, edge.isInternal); + edges.add(edge); + if (edge == lastEdge) { + break; + } + edge = edgeBuffer.getNext(edge); + } + + // Modify the edge buffer + edgeBuffer.replace(edges, returnEdge); + } + + // Now, construct the polygon + if (testPoint != null && holes != null && holes.size() > 0) { + // No holes, for test + final GeoPolygon testPolygon = new GeoConvexPolygon(planetModel, points, null, internalEdges, returnIsInternal); + if (testPolygon.isWithin(testPoint)) { + return null; + } + } + + final GeoPolygon realPolygon = new GeoConvexPolygon(planetModel, points, holes, internalEdges, returnIsInternal); + if (testPoint != null && (holes == null || holes.size() == 0)) { + if (realPolygon.isWithin(testPoint)) { + return null; + } + } + + rval.addShape(realPolygon); + return true; } - - /** Check if a point is within a described list of planes. - *@param newPoint is the point. - *@param currentPlanes is the list of planes. - *@return true if within. - */ - protected static boolean isWithin(final GeoPoint newPoint, final List currentPlanes) { - for (SidedPlane p : currentPlanes) { - if (!p.isWithin(newPoint)) + + protected static boolean isWithin(final GeoPoint point, final Set edgeSet, final Edge extension, final SidedPlane returnBoundary) { + if (!extension.plane.isWithin(point)) { + return false; + } + if (returnBoundary != null && !returnBoundary.isWithin(point)) { + return false; + } + return isWithin(point, edgeSet); + } + + protected static boolean isWithin(final GeoPoint point, final Set edgeSet) { + for (final Edge edge : edgeSet) { + if (!edge.plane.isWithin(point)) { return false; + } } return true; } - + /** Convert raw point index into valid array position. *@param index is the array index. *@param size is the array size. @@ -184,4 +672,219 @@ public class GeoPolygonFactory { return index; } + /** Class representing a single (unused) edge. + */ + protected static class Edge { + public final SidedPlane plane; + public final GeoPoint startPoint; + public final GeoPoint endPoint; + public final boolean isInternal; + + public Edge(final GeoPoint startPoint, final GeoPoint endPoint, final SidedPlane plane, final boolean isInternal) { + this.startPoint = startPoint; + this.endPoint = endPoint; + this.plane = plane; + this.isInternal = isInternal; + } + + @Override + public int hashCode() { + return System.identityHashCode(this); + } + + @Override + public boolean equals(final Object o) { + return o == this; + } + } + + /** Class representing an iterator over an EdgeBuffer. + */ + protected static class EdgeBufferIterator implements Iterator { + protected final EdgeBuffer edgeBuffer; + protected final Edge firstEdge; + protected Edge currentEdge; + + public EdgeBufferIterator(final EdgeBuffer edgeBuffer) { + this.edgeBuffer = edgeBuffer; + this.currentEdge = edgeBuffer.pickOne(); + this.firstEdge = currentEdge; + } + + @Override + public boolean hasNext() { + return currentEdge != null; + } + + @Override + public Edge next() { + final Edge rval = currentEdge; + if (currentEdge != null) { + currentEdge = edgeBuffer.getNext(currentEdge); + if (currentEdge == firstEdge) { + currentEdge = null; + } + } + return rval; + } + + @Override + public void remove() { + throw new RuntimeException("Unsupported operation"); + } + } + + /** Class representing a pool of unused edges, all linked together by vertices. + */ + protected static class EdgeBuffer { + protected Edge oneEdge; + protected final Set edges = new HashSet<>(); + protected final Map previousEdges = new HashMap<>(); + protected final Map nextEdges = new HashMap<>(); + + public EdgeBuffer(final List pointList, final int startPlaneStartIndex, final int startPlaneEndIndex, final SidedPlane startPlane, final boolean startPlaneIsInternal) { + /* + System.out.println("Initial points:"); + for (final GeoPoint p : pointList) { + System.out.println(" "+p); + } + System.out.println("For start plane, the following points are in/out:"); + for (final GeoPoint p: pointList) { + System.out.println(" "+p+" is: "+(startPlane.isWithin(p)?"in":"out")); + } + */ + + final Edge startEdge = new Edge(pointList.get(startPlaneStartIndex), pointList.get(startPlaneEndIndex), startPlane, startPlaneIsInternal); + // Fill in the EdgeBuffer by walking around creating more stuff + Edge currentEdge = startEdge; + int startIndex = startPlaneStartIndex; + int endIndex = startPlaneEndIndex; + boolean isInternal = startPlaneIsInternal; + while (true) { + // Compute the next edge + startIndex = endIndex; + endIndex++; + if (endIndex >= pointList.size()) { + endIndex -= pointList.size(); + } + // Get the next point + final GeoPoint newPoint = pointList.get(endIndex); + // Build the new edge + final boolean isNewPointWithin = currentEdge.plane.isWithin(newPoint); + final SidedPlane newPlane = new SidedPlane(currentEdge.startPoint, isNewPointWithin, pointList.get(startIndex), newPoint); + /* + System.out.println("For next plane, the following points are in/out:"); + for (final GeoPoint p: pointList) { + System.out.println(" "+p+" is: "+(newPlane.isWithin(p)?"in":"out")); + } + */ + final Edge newEdge = new Edge(pointList.get(startIndex), pointList.get(endIndex), newPlane, false); + + // Link it in + previousEdges.put(newEdge, currentEdge); + nextEdges.put(currentEdge, newEdge); + edges.add(newEdge); + currentEdge = newEdge; + + if (currentEdge.endPoint == startEdge.startPoint) { + // We finish here. Link the current edge to the start edge, and exit + previousEdges.put(startEdge, currentEdge); + nextEdges.put(currentEdge, startEdge); + edges.add(startEdge); + break; + } + } + + oneEdge = startEdge; + + // Verify the structure. + //verify(); + } + + /* + protected void verify() { + if (edges.size() != previousEdges.size() || edges.size() != nextEdges.size()) { + throw new IllegalStateException("broken structure"); + } + // Confirm each edge + for (final Edge e : edges) { + final Edge previousEdge = getPrevious(e); + final Edge nextEdge = getNext(e); + if (e.endPoint != nextEdge.startPoint) { + throw new IllegalStateException("broken structure"); + } + if (e.startPoint != previousEdge.endPoint) { + throw new IllegalStateException("broken structure"); + } + if (getNext(previousEdge) != e) { + throw new IllegalStateException("broken structure"); + } + if (getPrevious(nextEdge) != e) { + throw new IllegalStateException("broken structure"); + } + } + if (oneEdge != null && !edges.contains(oneEdge)) { + throw new IllegalStateException("broken structure"); + } + if (oneEdge == null && edges.size() > 0) { + throw new IllegalStateException("broken structure"); + } + } + */ + + public Edge getPrevious(final Edge currentEdge) { + return previousEdges.get(currentEdge); + } + + public Edge getNext(final Edge currentEdge) { + return nextEdges.get(currentEdge); + } + + public void replace(final List removeList, final Edge newEdge) { + /* + System.out.println("Replacing: "); + for (final Edge e : removeList) { + System.out.println(" "+e.startPoint+"-->"+e.endPoint); + } + System.out.println("...with: "+newEdge.startPoint+"-->"+newEdge.endPoint); + */ + final Edge previous = previousEdges.get(removeList.get(0)); + final Edge next = nextEdges.get(removeList.get(removeList.size()-1)); + edges.add(newEdge); + previousEdges.put(newEdge, previous); + nextEdges.put(previous, newEdge); + previousEdges.put(next, newEdge); + nextEdges.put(newEdge, next); + for (final Edge edge : removeList) { + if (edge == oneEdge) { + oneEdge = newEdge; + } + edges.remove(edge); + previousEdges.remove(edge); + nextEdges.remove(edge); + } + //verify(); + } + + public void clear() { + edges.clear(); + previousEdges.clear(); + nextEdges.clear(); + oneEdge = null; + } + + public int size() { + return edges.size(); + } + + public Iterator iterator() { + return new EdgeBufferIterator(this); + } + + public Edge pickOne() { + return oneEdge; + } + + } + } diff --git a/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/SidedPlane.java b/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/SidedPlane.java index e080bc04f40..0c19a9eadf4 100755 --- a/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/SidedPlane.java +++ b/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/SidedPlane.java @@ -31,7 +31,7 @@ public class SidedPlane extends Plane implements Membership { * * @param sidedPlane is the existing plane. */ - public SidedPlane(SidedPlane sidedPlane) { + public SidedPlane(final SidedPlane sidedPlane) { super(sidedPlane, sidedPlane.D); this.sigNum = -sidedPlane.sigNum; } @@ -44,13 +44,29 @@ public class SidedPlane extends Plane implements Membership { * @param A is the first in-plane point * @param B is the second in-plane point */ - public SidedPlane(Vector p, Vector A, Vector B) { + public SidedPlane(final Vector p, final Vector A, final Vector B) { super(A, B); sigNum = Math.signum(evaluate(p)); if (sigNum == 0.0) throw new IllegalArgumentException("Cannot determine sidedness because check point is on plane."); } + /** + * Construct a sided plane from a pair of vectors describing points, and including + * origin, plus a point p which describes the side. + * + * @param p point to evaluate + * @param onSide is true if the point is on the correct side of the plane, false otherwise. + * @param A is the first in-plane point + * @param B is the second in-plane point + */ + public SidedPlane(final Vector p, final boolean onSide, final Vector A, final Vector B) { + super(A, B); + sigNum = onSide?Math.signum(evaluate(p)):-Math.signum(evaluate(p)); + if (sigNum == 0.0) + throw new IllegalArgumentException("Cannot determine sidedness because check point is on plane."); + } + /** * Construct a sided plane from a point and a Z coordinate. * diff --git a/lucene/spatial3d/src/test/org/apache/lucene/spatial3d/geom/GeoPolygonTest.java b/lucene/spatial3d/src/test/org/apache/lucene/spatial3d/geom/GeoPolygonTest.java index d9b220de925..9f2bea41b48 100755 --- a/lucene/spatial3d/src/test/org/apache/lucene/spatial3d/geom/GeoPolygonTest.java +++ b/lucene/spatial3d/src/test/org/apache/lucene/spatial3d/geom/GeoPolygonTest.java @@ -18,8 +18,10 @@ package org.apache.lucene.spatial3d.geom; import java.util.ArrayList; import java.util.List; +import java.util.BitSet; import org.junit.Test; +import org.junit.Ignore; import static org.junit.Assert.assertEquals; import static org.junit.Assert.assertFalse; @@ -27,10 +29,45 @@ import static org.junit.Assert.assertTrue; public class GeoPolygonTest { + @Test + public void testPolygonClockwise() { + GeoPolygon c; + GeoPoint gp; + List points; + + // Points go counterclockwise, so + points = new ArrayList(); + points.add(new GeoPoint(PlanetModel.SPHERE, -0.1, -0.5)); + points.add(new GeoPoint(PlanetModel.SPHERE, 0.0, -0.6)); + points.add(new GeoPoint(PlanetModel.SPHERE, 0.1, -0.5)); + points.add(new GeoPoint(PlanetModel.SPHERE, 0.0, -0.4)); + + c = GeoPolygonFactory.makeGeoPolygon(PlanetModel.SPHERE, points); + //System.out.println(c); + + // Middle point should NOT be within!! + gp = new GeoPoint(PlanetModel.SPHERE, 0.0, -0.5); + assertTrue(!c.isWithin(gp)); + + // Now, go clockwise + points = new ArrayList(); + points.add(new GeoPoint(PlanetModel.SPHERE, 0.0, -0.4)); + points.add(new GeoPoint(PlanetModel.SPHERE, 0.1, -0.5)); + points.add(new GeoPoint(PlanetModel.SPHERE, 0.0, -0.6)); + points.add(new GeoPoint(PlanetModel.SPHERE, -0.1, -0.5)); + + c = GeoPolygonFactory.makeGeoPolygon(PlanetModel.SPHERE, points); + //System.out.println(c); + + // Middle point should be within!! + gp = new GeoPoint(PlanetModel.SPHERE, 0.0, -0.5); + assertTrue(c.isWithin(gp)); + + } @Test public void testPolygonPointWithin() { - GeoMembershipShape c; + GeoPolygon c; GeoPoint gp; List points; @@ -162,4 +199,141 @@ public class GeoPolygonTest { assertEquals(0.1, b.getMaxLatitude(), 0.000001); } + @Test + public void testPolygonBoundsCase1() { + GeoPolygon c; + LatLonBounds b; + List points; + XYZBounds xyzb; + GeoPoint point1; + GeoPoint point2; + GeoArea area; + + // Build the polygon + points = new ArrayList<>(); + points.add(new GeoPoint(PlanetModel.WGS84, 0.7769776943105245, -2.157536559188766)); + points.add(new GeoPoint(PlanetModel.WGS84, -0.9796549195552824, -0.25078026625235256)); + points.add(new GeoPoint(PlanetModel.WGS84, 0.17644522781457245, 2.4225312555674967)); + points.add(new GeoPoint(PlanetModel.WGS84, -1.4459804612164617, -1.2970934639728127)); + c = GeoPolygonFactory.makeGeoPolygon(PlanetModel.WGS84, points, 3); + // GeoCompositeMembershipShape: {[GeoConvexPolygon: {planetmodel=PlanetModel.WGS84, points= + // [[lat=0.17644522781457245, lon=2.4225312555674967], + // [lat=-1.4459804612164617, lon=-1.2970934639728127], + // [lat=0.7769776943105245, lon=-2.157536559188766]]}, + // GeoConcavePolygon: {planetmodel=PlanetModel.WGS84, points= + // [[lat=-0.9796549195552824, lon=-0.25078026625235256], + // [lat=0.17644522781457245, lon=2.4225312555674967], + // [lat=0.7769776943105245, lon=-2.157536559188766]]}]} + point1 = new GeoPoint(PlanetModel.WGS84, -1.2013743680763862, 0.48458963747230094); + point2 = new GeoPoint(0.3189285805649921, 0.16790264636909197, -0.9308557496413026); + + assertTrue(c.isWithin(point1)); + assertTrue(c.isWithin(point2)); + + // Now try bounds + xyzb = new XYZBounds(); + c.getBounds(xyzb); + area = GeoAreaFactory.makeGeoArea(PlanetModel.WGS84, + xyzb.getMinimumX(), xyzb.getMaximumX(), xyzb.getMinimumY(), xyzb.getMaximumY(), xyzb.getMinimumZ(), xyzb.getMaximumZ()); + + assertTrue(area.isWithin(point1)); + assertTrue(area.isWithin(point2)); + } + + @Test + public void testGeoPolygonBoundsCase2() { + // [junit4] 1> TEST: iter=23 shape=GeoCompositeMembershipShape: {[GeoConvexPolygon: {planetmodel=PlanetModel(ab=0.7563871189161702 c=1.2436128810838298), points= + // [[lat=0.014071770744627236, lon=0.011030818292803128], + // [lat=0.006772117088906782, lon=-0.0012531892445234592], + // [lat=0.0022201615609504792, lon=0.005941293187389326]]}, GeoConcavePolygon: {planetmodel=PlanetModel(ab=0.7563871189161702 c=1.2436128810838298), points= + // [[lat=-0.005507100238396111, lon=-0.008487706131259667], + // [lat=0.014071770744627236, lon=0.011030818292803128], + // [lat=0.0022201615609504792, lon=0.005941293187389326]]}]} + + PlanetModel pm = new PlanetModel(0.7563871189161702, 1.2436128810838298); + // Build the polygon + GeoCompositeMembershipShape c = new GeoCompositeMembershipShape(); + List points1 = new ArrayList<>(); + points1.add(new GeoPoint(pm, 0.014071770744627236, 0.011030818292803128)); + points1.add(new GeoPoint(pm, 0.006772117088906782, -0.0012531892445234592)); + points1.add(new GeoPoint(pm, 0.0022201615609504792, 0.005941293187389326)); + BitSet p1bits = new BitSet(); + c.addShape(new GeoConvexPolygon(pm, points1, p1bits, true)); + List points2 = new ArrayList<>(); + points2.add(new GeoPoint(pm, -0.005507100238396111, -0.008487706131259667)); + points2.add(new GeoPoint(pm, 0.014071770744627236, 0.011030818292803128)); + points2.add(new GeoPoint(pm, 0.0022201615609504792, 0.005941293187389326)); + BitSet p2bits = new BitSet(); + p2bits.set(1, true); + c.addShape(new GeoConcavePolygon(pm, points2, p2bits, false)); + //System.out.println(c); + + // [junit4] 1> point=[lat=0.003540694517552105, lon=-9.99517927901697E-4] + // [junit4] 1> quantized=[X=0.7563849869428783, Y=-7.560204674780763E-4, Z=0.0026781405884151086] + GeoPoint point = new GeoPoint(pm, 0.003540694517552105, -9.99517927901697E-4); + GeoPoint pointQuantized = new GeoPoint(0.7563849869428783, -7.560204674780763E-4, 0.0026781405884151086); + + // Now try bounds + XYZBounds xyzb = new XYZBounds(); + c.getBounds(xyzb); + GeoArea area = GeoAreaFactory.makeGeoArea(pm, + xyzb.getMinimumX(), xyzb.getMaximumX(), xyzb.getMinimumY(), xyzb.getMaximumY(), xyzb.getMinimumZ(), xyzb.getMaximumZ()); + + assertTrue(c.isWithin(point)); + assertTrue(c.isWithin(pointQuantized)); + // This fails!! + assertTrue(area.isWithin(point)); + assertTrue(area.isWithin(pointQuantized)); + } + + @Test + public void testGeoConcaveRelationshipCase1() { + /* + [junit4] 1> doc=906 matched but should not + [junit4] 1> point=[lat=-0.9825762558001477, lon=2.4832136904725273] + [junit4] 1> quantized=[X=-0.4505446160475436, Y=0.34850109186970535, Z=-0.8539966368663765] + +doc=906 added here: + + [junit4] 1> cycle: cell=107836 parentCellID=107835 x: -1147288468 TO -742350917, y: -1609508490 TO 1609508490, z: -2147483647 TO 2147483647, splits: 3 queue.size()=1 + [junit4] 1> minx=-0.6107484000858642 maxx=-0.39518364125756916 miny=-0.8568069517709872 maxy=0.8568069517709872 minz=-1.1431930485939341 maxz=1.1431930485939341 + [junit4] 1> GeoArea.CONTAINS: now addAll + +shape: + [junit4] 1> TEST: iter=18 shape=GeoCompositeMembershipShape: {[GeoConvexPolygon: { + planetmodel=PlanetModel(ab=0.8568069516722363 c=1.1431930483277637), points= + [[lat=1.1577814487635816, lon=1.6283601832010004], + [lat=0.6664570999069251, lon=2.0855825542851574], + [lat=-0.23953537010974632, lon=1.8498724094352876]]}, GeoConcavePolygon: {planetmodel=PlanetModel(ab=0.8568069516722363 c=1.1431930483277637), points= + [[lat=1.1577814487635816, lon=1.6283601832010004], + [lat=-0.23953537010974632, lon=1.8498724094352876], + [lat=-1.1766904875978805, lon=-2.1346828411344436]]}]} + */ + PlanetModel pm = new PlanetModel(0.8568069516722363, 1.1431930483277637); + // Build the polygon + GeoCompositeMembershipShape c = new GeoCompositeMembershipShape(); + List points1 = new ArrayList<>(); + points1.add(new GeoPoint(pm, 1.1577814487635816, 1.6283601832010004)); + points1.add(new GeoPoint(pm, 0.6664570999069251, 2.0855825542851574)); + points1.add(new GeoPoint(pm, -0.23953537010974632, 1.8498724094352876)); + BitSet p1bits = new BitSet(); + c.addShape(new GeoConvexPolygon(pm, points1, p1bits, true)); + List points2 = new ArrayList<>(); + points2.add(new GeoPoint(pm, 1.1577814487635816, 1.6283601832010004)); + points2.add(new GeoPoint(pm, -0.23953537010974632, 1.8498724094352876)); + points2.add(new GeoPoint(pm, -1.1766904875978805, -2.1346828411344436)); + BitSet p2bits = new BitSet(); + p2bits.set(1, true); + c.addShape(new GeoConcavePolygon(pm, points2, p2bits, false)); + //System.out.println(c); + + GeoPoint point = new GeoPoint(pm, -0.9825762558001477, 2.4832136904725273); + GeoPoint quantizedPoint = new GeoPoint(-0.4505446160475436, 0.34850109186970535, -0.8539966368663765); + + GeoArea xyzSolid = GeoAreaFactory.makeGeoArea(pm, + -0.6107484000858642, -0.39518364125756916, -0.8568069517709872, 0.8568069517709872, -1.1431930485939341, 1.1431930485939341); + //System.out.println("relationship = "+xyzSolid.getRelationship(c)); + assertTrue(xyzSolid.getRelationship(c) == GeoArea.OVERLAPS); + } + }