+ * Construction takes {@code O(n log n)} time for sorting and tree construction. + * {@link #relate relate()} are {@code O(n)}, but for most + * practical lines and polygons are much faster than brute force. + * @lucene.internal + */ +public abstract class EdgeTree { + /** minimum latitude of this geometry's bounding box area */ + public final double minLat; + /** maximum latitude of this geometry's bounding box area */ + public final double maxLat; + /** minimum longitude of this geometry's bounding box area */ + public final double minLon; + /** maximum longitude of this geometry's bounding box area */ + public final double maxLon; + + // each component is a node in an augmented 2d kd-tree: we alternate splitting between latitude/longitude, + // and pull up max values for both dimensions to each parent node (regardless of split). + + /** maximum latitude of this component or any of its children */ + protected double maxY; + /** maximum longitude of this component or any of its children */ + protected double maxX; + /** which dimension was this node split on */ + // TODO: its implicit based on level, but boolean keeps code simple + protected boolean splitX; + + // child components, or null + protected EdgeTree left; + protected EdgeTree right; + + /** root node of edge tree */ + protected final Edge tree; + + protected EdgeTree(final double minLat, final double maxLat, final double minLon, final double maxLon, double[] lats, double[] lons) { + this.minLat = minLat; + this.maxLat = maxLat; + this.minLon = minLon; + this.maxLon = maxLon; + this.maxY = maxLat; + this.maxX = maxLon; + + // create interval tree of edges + this.tree = createTree(lats, lons); + } + + /** Returns relation to the provided triangle */ + public Relation relateTriangle(double ax, double ay, double bx, double by, double cx, double cy) { + // compute bounding box of triangle + double minLat = StrictMath.min(StrictMath.min(ay, by), cy); + double minLon = StrictMath.min(StrictMath.min(ax, bx), cx); + double maxLat = StrictMath.max(StrictMath.max(ay, by), cy); + double maxLon = StrictMath.max(StrictMath.max(ax, bx), cx); + if (minLat <= maxY && minLon <= maxX) { + Relation relation = internalComponentRelateTriangle(ax, ay, bx, by, cx, cy); + if (relation != Relation.CELL_OUTSIDE_QUERY) { + return relation; + } + if (left != null) { + relation = left.relateTriangle(ax, ay, bx, by, cx, cy); + if (relation != Relation.CELL_OUTSIDE_QUERY) { + return relation; + } + } + if (right != null && ((splitX == false && maxLat >= this.minLat) || (splitX && maxLon >= this.minLon))) { + relation = right.relateTriangle(ax, ay, bx, by, cx, cy); + if (relation != Relation.CELL_OUTSIDE_QUERY) { + return relation; + } + } + } + return Relation.CELL_OUTSIDE_QUERY; + } + + /** Returns relation to the provided rectangle */ + public Relation relate(double minLat, double maxLat, double minLon, double maxLon) { + if (minLat <= maxY && minLon <= maxX) { + Relation relation = internalComponentRelate(minLat, maxLat, minLon, maxLon); + if (relation != Relation.CELL_OUTSIDE_QUERY) { + return relation; + } + if (left != null) { + relation = left.relate(minLat, maxLat, minLon, maxLon); + if (relation != Relation.CELL_OUTSIDE_QUERY) { + return relation; + } + } + if (right != null && ((splitX == false && maxLat >= this.minLat) || (splitX && maxLon >= this.minLon))) { + relation = right.relate(minLat, maxLat, minLon, maxLon); + if (relation != Relation.CELL_OUTSIDE_QUERY) { + return relation; + } + } + } + return Relation.CELL_OUTSIDE_QUERY; + } + + protected Relation componentRelate(double minLat, double maxLat, double minLon, double maxLon) { + return null; + } + protected Relation componentRelateTriangle(double ax, double ay, double bx, double by, double cx, double cy) { + return null; + } + + private Relation internalComponentRelateTriangle(double ax, double ay, double bx, double by, double cx, double cy) { + // compute bounding box of triangle + double minLat = StrictMath.min(StrictMath.min(ay, by), cy); + double minLon = StrictMath.min(StrictMath.min(ax, bx), cx); + double maxLat = StrictMath.max(StrictMath.max(ay, by), cy); + double maxLon = StrictMath.max(StrictMath.max(ax, bx), cx); + if (maxLon < this.minLon || minLon > this.maxLon || maxLat < this.minLat || minLat > this.maxLat) { + return Relation.CELL_OUTSIDE_QUERY; + } + + Relation shapeRelation = componentRelateTriangle(ax, ay, bx, by, cx, cy); + if (shapeRelation != null) { + return shapeRelation; + } + + // we cross + if (tree.crossesTriangle(ax, ay, bx, by, cx, cy)) { + return Relation.CELL_CROSSES_QUERY; + } + return Relation.CELL_OUTSIDE_QUERY; + } + + + /** Returns relation to the provided rectangle for this component */ + protected Relation internalComponentRelate(double minLat, double maxLat, double minLon, double maxLon) { + // if the bounding boxes are disjoint then the shape does not cross + if (maxLon < this.minLon || minLon > this.maxLon || maxLat < this.minLat || minLat > this.maxLat) { + return Relation.CELL_OUTSIDE_QUERY; + } + // if the rectangle fully encloses us, we cross. + if (minLat <= this.minLat && maxLat >= this.maxLat && minLon <= this.minLon && maxLon >= this.maxLon) { + return Relation.CELL_CROSSES_QUERY; + } + + Relation shapeRelation = componentRelate(minLat, maxLat, minLon, maxLon); + if (shapeRelation != null) { + return shapeRelation; + } + + // we cross + if (tree.crosses(minLat, maxLat, minLon, maxLon)) { + return Relation.CELL_CROSSES_QUERY; + } + + return Relation.CELL_OUTSIDE_QUERY; + } + + /** Creates tree from sorted components (with range low and high inclusive) */ + protected static EdgeTree createTree(EdgeTree components[], int low, int high, boolean splitX) { + if (low > high) { + return null; + } + final int mid = (low + high) >>> 1; + if (low < high) { + Comparator comparator; + if (splitX) { + comparator = (left, right) -> { + int ret = Double.compare(left.minLon, right.minLon); + if (ret == 0) { + ret = Double.compare(left.maxX, right.maxX); + } + return ret; + }; + } else { + comparator = (left, right) -> { + int ret = Double.compare(left.minLat, right.minLat); + if (ret == 0) { + ret = Double.compare(left.maxY, right.maxY); + } + return ret; + }; + } + ArrayUtil.select(components, low, high + 1, mid, comparator); + } + // add midpoint + EdgeTree newNode = components[mid]; + newNode.splitX = splitX; + // add children + newNode.left = createTree(components, low, mid - 1, !splitX); + newNode.right = createTree(components, mid + 1, high, !splitX); + // pull up max values to this node + if (newNode.left != null) { + newNode.maxX = Math.max(newNode.maxX, newNode.left.maxX); + newNode.maxY = Math.max(newNode.maxY, newNode.left.maxY); + } + if (newNode.right != null) { + newNode.maxX = Math.max(newNode.maxX, newNode.right.maxX); + newNode.maxY = Math.max(newNode.maxY, newNode.right.maxY); + } + return newNode; + } + + /** + * Internal tree node: represents geometry edge from lat1,lon1 to lat2,lon2. + * The sort value is {@code low}, which is the minimum latitude of the edge. + * {@code max} stores the maximum latitude of this edge or any children. + */ + static class Edge { + // lat-lon pair (in original order) of the two vertices + final double lat1, lat2; + final double lon1, lon2; + /** min of this edge */ + final double low; + /** max latitude of this edge or any children */ + double max; + + /** left child edge, or null */ + Edge left; + /** right child edge, or null */ + Edge right; + + Edge(double lat1, double lon1, double lat2, double lon2, double low, double max) { + this.lat1 = lat1; + this.lon1 = lon1; + this.lat2 = lat2; + this.lon2 = lon2; + this.low = low; + this.max = max; + } + + /** Returns true if the triangle crosses any edge in this edge subtree */ + boolean crossesTriangle(double ax, double ay, double bx, double by, double cx, double cy) { + // compute bounding box of triangle + double minLat = StrictMath.min(StrictMath.min(ay, by), cy); + double minLon = StrictMath.min(StrictMath.min(ax, bx), cx); + double maxLat = StrictMath.max(StrictMath.max(ay, by), cy); + double maxLon = StrictMath.max(StrictMath.max(ax, bx), cx); + + if (minLat <= max) { + double dy = lat1; + double ey = lat2; + double dx = lon1; + double ex = lon2; + + // optimization: see if the rectangle is outside of the "bounding box" of the polyline at all + // if not, don't waste our time trying more complicated stuff + boolean outside = (dy < minLat && ey < minLat) || + (dy > maxLat && ey > maxLat) || + (dx < minLon && ex < minLon) || + (dx > maxLon && ex > maxLon); + + if (outside == false) { + // does triangle's first edge intersect polyline? + // ax, ay -> bx, by + if (lineCrossesLine(ax, ay, bx, by, dx, dy, ex, ey)) { + return true; + } + + // does triangle's second edge intersect polyline? + // bx, by -> cx, cy + if (lineCrossesLine(bx, by, cx, cy, dx, dy, ex, ey)) { + return true; + } + + // does triangle's third edge intersect polyline? + // cx, cy -> ax, ay + if (lineCrossesLine(cx, cy, ax, ay, dx, dy, ex, ey)) { + return true; + } + } + + if (left != null) { + if (left.crossesTriangle(ax, ay, bx, by, cx, cy)) { + return true; + } + } + + if (right != null && maxLat >= low) { + if (right.crossesTriangle(ax, ay, bx, by, cx, cy)) { + return true; + } + } + } + return false; + } + + /** Returns true if the box crosses any edge in this edge subtree */ + boolean crosses(double minLat, double maxLat, double minLon, double maxLon) { + // we just have to cross one edge to answer the question, so we descend the tree and return when we do. + if (minLat <= max) { + // we compute line intersections of every polygon edge with every box line. + // if we find one, return true. + // for each box line (AB): + // for each poly line (CD): + // intersects = orient(C,D,A) * orient(C,D,B) <= 0 && orient(A,B,C) * orient(A,B,D) <= 0 + double cy = lat1; + double dy = lat2; + double cx = lon1; + double dx = lon2; + + // optimization: see if the rectangle is outside of the "bounding box" of the polyline at all + // if not, don't waste our time trying more complicated stuff + boolean outside = (cy < minLat && dy < minLat) || + (cy > maxLat && dy > maxLat) || + (cx < minLon && dx < minLon) || + (cx > maxLon && dx > maxLon); + // optimization: see if either end of the line segment is contained by the rectangle + if (Rectangle.containsPoint(cy, cx, minLat, maxLat, minLon, maxLon) + || Rectangle.containsPoint(dy, dx, minLat, maxLat, minLon, maxLon)) { + return true; + } + + if (outside == false) { + // does box's top edge intersect polyline? + // ax = minLon, bx = maxLon, ay = maxLat, by = maxLat + if (orient(cx, cy, dx, dy, minLon, maxLat) * orient(cx, cy, dx, dy, maxLon, maxLat) <= 0 && + orient(minLon, maxLat, maxLon, maxLat, cx, cy) * orient(minLon, maxLat, maxLon, maxLat, dx, dy) <= 0) { + return true; + } + + // does box's right edge intersect polyline? + // ax = maxLon, bx = maxLon, ay = maxLat, by = minLat + if (orient(cx, cy, dx, dy, maxLon, maxLat) * orient(cx, cy, dx, dy, maxLon, minLat) <= 0 && + orient(maxLon, maxLat, maxLon, minLat, cx, cy) * orient(maxLon, maxLat, maxLon, minLat, dx, dy) <= 0) { + return true; + } + + // does box's bottom edge intersect polyline? + // ax = maxLon, bx = minLon, ay = minLat, by = minLat + if (orient(cx, cy, dx, dy, maxLon, minLat) * orient(cx, cy, dx, dy, minLon, minLat) <= 0 && + orient(maxLon, minLat, minLon, minLat, cx, cy) * orient(maxLon, minLat, minLon, minLat, dx, dy) <= 0) { + return true; + } + + // does box's left edge intersect polyline? + // ax = minLon, bx = minLon, ay = minLat, by = maxLat + if (orient(cx, cy, dx, dy, minLon, minLat) * orient(cx, cy, dx, dy, minLon, maxLat) <= 0 && + orient(minLon, minLat, minLon, maxLat, cx, cy) * orient(minLon, minLat, minLon, maxLat, dx, dy) <= 0) { + return true; + } + } + + if (left != null) { + if (left.crosses(minLat, maxLat, minLon, maxLon)) { + return true; + } + } + + if (right != null && maxLat >= low) { + if (right.crosses(minLat, maxLat, minLon, maxLon)) { + return true; + } + } + } + return false; + } + } + + /** + * Creates an edge interval tree from a set of geometry vertices. + * @return root node of the tree. + */ + private static Edge createTree(double[] lats, double[] lons) { + Edge edges[] = new Edge[lats.length - 1]; + for (int i = 1; i < lats.length; i++) { + double lat1 = lats[i-1]; + double lon1 = lons[i-1]; + double lat2 = lats[i]; + double lon2 = lons[i]; + edges[i - 1] = new Edge(lat1, lon1, lat2, lon2, Math.min(lat1, lat2), Math.max(lat1, lat2)); + } + // sort the edges then build a balanced tree from them + Arrays.sort(edges, (left, right) -> { + int ret = Double.compare(left.low, right.low); + if (ret == 0) { + ret = Double.compare(left.max, right.max); + } + return ret; + }); + return createTree(edges, 0, edges.length - 1); + } + + /** Creates tree from sorted edges (with range low and high inclusive) */ + private static Edge createTree(Edge edges[], int low, int high) { + if (low > high) { + return null; + } + // add midpoint + int mid = (low + high) >>> 1; + Edge newNode = edges[mid]; + // add children + newNode.left = createTree(edges, low, mid - 1); + newNode.right = createTree(edges, mid + 1, high); + // pull up max values to this node + if (newNode.left != null) { + newNode.max = Math.max(newNode.max, newNode.left.max); + } + if (newNode.right != null) { + newNode.max = Math.max(newNode.max, newNode.right.max); + } + return newNode; + } +} http://git-wip-us.apache.org/repos/asf/lucene-solr/blob/0cbefe8b/lucene/core/src/java/org/apache/lucene/geo/GeoUtils.java ---------------------------------------------------------------------- diff --git a/lucene/core/src/java/org/apache/lucene/geo/GeoUtils.java b/lucene/core/src/java/org/apache/lucene/geo/GeoUtils.java index 468de93..0c73032 100644 --- a/lucene/core/src/java/org/apache/lucene/geo/GeoUtils.java +++ b/lucene/core/src/java/org/apache/lucene/geo/GeoUtils.java @@ -194,6 +194,20 @@ public final class GeoUtils { } } + /** uses orient method to compute whether two line segments cross */ + public static boolean lineCrossesLine(double a1x, double a1y, double b1x, double b1y, double a2x, double a2y, double b2x, double b2y) { + // shortcut: either "line" is actually a point + if ((a1x == b1x && a1y == b1y) || (a2x == b2x && a2y == b2y)) { + return false; + } + + if (orient(a2x, a2y, b2x, b2y, a1x, a1y) * orient(a2x, a2y, b2x, b2y, b1x, b1y) <= 0 && + orient(a1x, a1y, b1x, b1y, a2x, a2y) * orient(a1x, a1y, b1x, b1y, b2x, b2y) <= 0) { + return true; + } + return false; + } + /** * used to define the orientation of 3 points * -1 = Clockwise http://git-wip-us.apache.org/repos/asf/lucene-solr/blob/0cbefe8b/lucene/core/src/java/org/apache/lucene/geo/Polygon2D.java ---------------------------------------------------------------------- diff --git a/lucene/core/src/java/org/apache/lucene/geo/Polygon2D.java b/lucene/core/src/java/org/apache/lucene/geo/Polygon2D.java index 64a3784..fee23d0 100644 --- a/lucene/core/src/java/org/apache/lucene/geo/Polygon2D.java +++ b/lucene/core/src/java/org/apache/lucene/geo/Polygon2D.java @@ -16,72 +16,29 @@ */ package org.apache.lucene.geo; -import java.util.Arrays; -import java.util.Comparator; - import org.apache.lucene.index.PointValues.Relation; -import org.apache.lucene.util.ArrayUtil; - -import static org.apache.lucene.geo.GeoUtils.orient; /** * 2D polygon implementation represented as a balanced interval tree of edges. *

- * Construction takes {@code O(n log n)} time for sorting and tree construction. - * {@link #contains contains()} and {@link #relate relate()} are {@code O(n)}, but for most - * practical polygons are much faster than brute force. - *

* Loosely based on the algorithm described in * http://www-ma2.upc.es/geoc/Schirra-pointPolygon.pdf. * @lucene.internal */ // Both Polygon.contains() and Polygon.crossesSlowly() loop all edges, and first check that the edge is within a range. -// we just organize the edges to do the same computations on the same subset of edges more efficiently. -public final class Polygon2D { - /** minimum latitude of this polygon's bounding box area */ - public final double minLat; - /** maximum latitude of this polygon's bounding box area */ - public final double maxLat; - /** minimum longitude of this polygon's bounding box area */ - public final double minLon; - /** maximum longitude of this polygon's bounding box area */ - public final double maxLon; - +// we just organize the edges to do the same computations on the same subset of edges more efficiently. +public final class Polygon2D extends EdgeTree { // each component/hole is a node in an augmented 2d kd-tree: we alternate splitting between latitude/longitude, // and pull up max values for both dimensions to each parent node (regardless of split). - - /** maximum latitude of this component or any of its children */ - private double maxY; - /** maximum longitude of this component or any of its children */ - private double maxX; - /** which dimension was this node split on */ - // TODO: its implicit based on level, but boolean keeps code simple - private boolean splitX; - - // child components, or null - private Polygon2D left; - private Polygon2D right; - /** tree of holes, or null */ private final Polygon2D holes; - - /** root node of edge tree */ - private final Edge tree; private Polygon2D(Polygon polygon, Polygon2D holes) { + super(polygon.minLat, polygon.maxLat, polygon.minLon, polygon.maxLon, polygon.getPolyLats(), polygon.getPolyLons()); this.holes = holes; - this.minLat = polygon.minLat; - this.maxLat = polygon.maxLat; - this.minLon = polygon.minLon; - this.maxLon = polygon.maxLon; - this.maxY = maxLat; - this.maxX = maxLon; - - // create interval tree of edges - this.tree = createTree(polygon.getPolyLats(), polygon.getPolyLons()); } - /** + /** * Returns true if the point is contained within this polygon. *

+ * See + * https://www.ecse.rpi.edu/~wrf/Research/Short_Notes/pnpoly.html for more information. */ - static final class Edge { - // lat-lon pair (in original order) of the two vertices - final double lat1, lat2; - final double lon1, lon2; - /** min of this edge */ - final double low; - /** max latitude of this edge or any children */ - double max; - - /** left child edge, or null */ - Edge left; - /** right child edge, or null */ - Edge right; - - Edge(double lat1, double lon1, double lat2, double lon2, double low, double max) { - this.lat1 = lat1; - this.lon1 = lon1; - this.lat2 = lat2; - this.lon2 = lon2; - this.low = low; - this.max = max; - } - - /** - * Returns true if the point crosses this edge subtree an odd number of times - *

+ * Note: + *

+ *
• {@code QueryRelation.WITHIN} queries are not yet supported
• + *
• Dateline crossing is not yet supported
• + *
+ *

+ * todo: + *

+ *
• Add distance support for buffered queries
• + *
+ *

The field must be indexed using + * {@link org.apache.lucene.document.LatLonShape#createIndexableFields} added per document. + * + * @lucene.experimental + **/ +final class LatLonShapeLineQuery extends LatLonShapeQuery { + final Line[] lines; + final private Line2D line2D; + + public LatLonShapeLineQuery(String field, QueryRelation queryRelation, Line... lines) { + super(field, queryRelation); + /** line queries do not support within relations, only intersects and disjoint */ + if (queryRelation == QueryRelation.WITHIN) { + throw new IllegalArgumentException("LatLonShapeLineQuery does not support " + QueryRelation.WITHIN + " queries"); + } + + if (lines == null) { + throw new IllegalArgumentException("lines must not be null"); + } + if (lines.length == 0) { + throw new IllegalArgumentException("lines must not be empty"); + } + for (int i = 0; i < lines.length; ++i) { + if (lines[i] == null) { + throw new IllegalArgumentException("line[" + i + "] must not be null"); + } else if (lines[i].minLon > lines[i].maxLon) { + throw new IllegalArgumentException("LatLonShapeLineQuery does not currently support querying across dateline."); + } + } + this.lines = lines.clone(); + this.line2D = Line2D.create(lines); + } + + @Override + protected Relation relateRangeBBoxToQuery(int minXOffset, int minYOffset, byte[] minTriangle, + int maxXOffset, int maxYOffset, byte[] maxTriangle) { + double minLat = GeoEncodingUtils.decodeLatitude(LatLonShape.decodeTriangleBoxVal(minTriangle, minYOffset)); + double minLon = GeoEncodingUtils.decodeLongitude(LatLonShape.decodeTriangleBoxVal(minTriangle, minXOffset)); + double maxLat = GeoEncodingUtils.decodeLatitude(LatLonShape.decodeTriangleBoxVal(maxTriangle, maxYOffset)); + double maxLon = GeoEncodingUtils.decodeLongitude(LatLonShape.decodeTriangleBoxVal(maxTriangle, maxXOffset)); + + // check internal node against query + return line2D.relate(minLat, maxLat, minLon, maxLon); + } + + @Override + protected boolean queryMatches(byte[] t) { + long a = NumericUtils.sortableBytesToLong(t, 4 * LatLonShape.BYTES); + long b = NumericUtils.sortableBytesToLong(t, 5 * LatLonShape.BYTES); + long c = NumericUtils.sortableBytesToLong(t, 6 * LatLonShape.BYTES); + + int aX = (int)((a >>> 32) & 0x00000000FFFFFFFFL); + int bX = (int)((b >>> 32) & 0x00000000FFFFFFFFL); + int cX = (int)((c >>> 32) & 0x00000000FFFFFFFFL); + int aY = (int)(a & 0x00000000FFFFFFFFL); + int bY = (int)(b & 0x00000000FFFFFFFFL); + int cY = (int)(c & 0x00000000FFFFFFFFL); + + double alat = GeoEncodingUtils.decodeLatitude(aY); + double alon = GeoEncodingUtils.decodeLongitude(aX); + double blat = GeoEncodingUtils.decodeLatitude(bY); + double blon = GeoEncodingUtils.decodeLongitude(bX); + double clat = GeoEncodingUtils.decodeLatitude(cY); + double clon = GeoEncodingUtils.decodeLongitude(cX); + + if (queryRelation == LatLonShape.QueryRelation.WITHIN) { + return line2D.relateTriangle(alon, alat, blon, blat, clon, clat) == Relation.CELL_INSIDE_QUERY; + } + // INTERSECTS + return line2D.relateTriangle(alon, alat, blon, blat, clon, clat) != Relation.CELL_OUTSIDE_QUERY; + } + + @Override + public String toString(String field) { + final StringBuilder sb = new StringBuilder(); + sb.append(getClass().getSimpleName()); + sb.append(':'); + if (this.field.equals(field) == false) { + sb.append(" field="); + sb.append(this.field); + sb.append(':'); + } + sb.append("Line(" + lines[0].toGeoJSON() + ")"); + return sb.toString(); + } + + @Override + protected boolean equalsTo(Object o) { + return super.equalsTo(o) && Arrays.equals(lines, ((LatLonShapeLineQuery)o).lines); + } + + @Override + public int hashCode() { + int hash = super.hashCode(); + hash = 31 * hash + Arrays.hashCode(lines); + return hash; + } +} http://git-wip-us.apache.org/repos/asf/lucene-solr/blob/0cbefe8b/lucene/sandbox/src/java/org/apache/lucene/document/LatLonShapePolygonQuery.java ---------------------------------------------------------------------- diff --git a/lucene/sandbox/src/java/org/apache/lucene/document/LatLonShapePolygonQuery.java b/lucene/sandbox/src/java/org/apache/lucene/document/LatLonShapePolygonQuery.java index a587112..2b342a8 100644 --- a/lucene/sandbox/src/java/org/apache/lucene/document/LatLonShapePolygonQuery.java +++ b/lucene/sandbox/src/java/org/apache/lucene/document/LatLonShapePolygonQuery.java @@ -29,7 +29,7 @@ import org.apache.lucene.util.NumericUtils; * Finds all previously indexed shapes that intersect the specified arbitrary. * *