Apache
/
lucene
mirror of https://github.com/apache/lucene.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

LUCENE-7157 

Add and use clockwise/counterclockwise construction of geo polygons

All code in place but isn't working

Add GeoConcavePolygon, and make more efficient intersection code for GeoConvexPolygon

Finish code, but it fails tests due to looping

Fix some problems

Fix typo

Get one test to run at least

Get it to work

Turn off debugging

Get rid of unused code

Fix beasting failure

Boundary edges have to use sided planes that are inverted in order to constrain the intersections right.

Add test

Use a more rigorous membership structure for bounding edges.

Fix bounding issue

Remove dependence on hash order

Get everything working again

Get everything working

Remove System.out
This commit is contained in:
Karl Wright 2016-04-01 07:07:05 -04:00
parent ef6a0d0013
commit 1040fcbf42
6 changed files with 1525 additions and 150 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
lucene/spatial3d/src/java/org/apache/lucene/spatial3d/Geo3DPoint.java
View File

@ -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);
}

393
lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoConcavePolygon.java Normal file
View File

@ -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<GeoPoint> 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<GeoPolygon> 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<SidedPlane, Membership> 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<GeoPoint> 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<GeoPoint> pointList, final List<GeoPolygon> 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<GeoPoint> 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<GeoPoint> pointList,
final List<GeoPolygon> 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<GeoPolygon> 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) + "}";
}
}

175
lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoConvexPolygon.java
View File

@ -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<GeoPoint> 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<GeoPolygon> 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<SidedPlane, Membership> 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<GeoPoint> 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<GeoPoint> pointList, final List<GeoPolygon> 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<GeoPoint> pointList, final BitSet internalEdgeFlags,
final boolean returnEdgeInternal) {
public GeoConvexPolygon(final PlanetModel planetModel,
final List<GeoPoint> 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<GeoPoint> pointList,
final List<GeoPolygon> 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<GeoPolygon> 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) + "}";
}
}

903
lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoPolygonFactory.java
View File

@ -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<GeoPoint> pointList, final int convexPointIndex) {
public static GeoPolygon makeGeoPolygon(final PlanetModel planetModel,
final List<GeoPoint> 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<GeoPoint> pointList,
final int convexPointIndex,
final List<GeoPolygon> 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<GeoPoint> 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<GeoPoint> pointList,
final List<GeoPolygon> 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<GeoPoint> pointList,
final List<GeoPolygon> 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<GeoPoint> 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<GeoPoint> 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<GeoPoint> pointsList,
final int startPointIndex,
final int endPointIndex,
final SidedPlane startingEdge,
final boolean isInternalEdge,
final List<GeoPolygon> 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<GeoPoint> recursionList = new ArrayList<GeoPoint>();
final List<GeoPoint> currentList = new ArrayList<GeoPoint>();
final BitSet internalEdgeList = new BitSet();
final List<SidedPlane> currentPlanes = new ArrayList<SidedPlane>();
// 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<GeoPolygon> 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<GeoPoint> points = new ArrayList<GeoPoint>(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<GeoPolygon> 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<Edge> 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<Edge> 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<Edge> 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<GeoPoint> points = new ArrayList<GeoPoint>(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<Edge> edges = new ArrayList<Edge>(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<SidedPlane> currentPlanes) {
for (SidedPlane p : currentPlanes) {
if (!p.isWithin(newPoint))
protected static boolean isWithin(final GeoPoint point, final Set<Edge> 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<Edge> 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<Edge> {
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<Edge> edges = new HashSet<>();
protected final Map<Edge, Edge> previousEdges = new HashMap<>();
protected final Map<Edge, Edge> nextEdges = new HashMap<>();
public EdgeBuffer(final List<GeoPoint> 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<Edge> 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<Edge> iterator() {
return new EdgeBufferIterator(this);
}
public Edge pickOne() {
return oneEdge;
}
}
}

20
lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/SidedPlane.java
View File

@ -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.
*

176
lucene/spatial3d/src/test/org/apache/lucene/spatial3d/geom/GeoPolygonTest.java
View File

@ -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<GeoPoint> points;
// Points go counterclockwise, so
points = new ArrayList<GeoPoint>();
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<GeoPoint>();
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<GeoPoint> points;
@ -162,4 +199,141 @@ public class GeoPolygonTest {
assertEquals(0.1, b.getMaxLatitude(), 0.000001);
}
@Test
public void testPolygonBoundsCase1() {
GeoPolygon c;
LatLonBounds b;
List<GeoPoint> 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<GeoPoint> 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<GeoPoint> 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<GeoPoint> 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<GeoPoint> 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);
}
}