diff --git a/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoComplexPolygon.java b/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoComplexPolygon.java index df89f55a7fb..9ffbcc720dc 100644 --- a/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoComplexPolygon.java +++ b/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoComplexPolygon.java @@ -120,8 +120,27 @@ class GeoComplexPolygon extends GeoBasePolygon { @Override public boolean intersects(final Plane p, final GeoPoint[] notablePoints, final Membership... bounds) { - // MHL - return false; + // Create the intersector + final EdgeIterator intersector = new IntersectorEdgeIterator(p, notablePoints, bounds); + // First, compute the bounds for the the plane + final XYZBounds xyzBounds = new XYZBounds(); + p.recordBounds(xyzBounds); + // Figure out which tree likely works best + final double xDelta = xyzBounds.getMaximumX() - xyzBounds.getMinimumX(); + final double yDelta = xyzBounds.getMaximumY() - xyzBounds.getMinimumY(); + final double zDelta = xyzBounds.getMaximumZ() - xyzBounds.getMinimumZ(); + // Select the smallest range + if (xDelta <= yDelta && xDelta <= zDelta) { + // Drill down in x + return !xtree.traverse(intersector, xyzBounds.getMinimumX(), xyzBounds.getMaximumX()); + } else if (yDelta <= xDelta && yDelta <= zDelta) { + // Drill down in y + return !ytree.traverse(intersector, xyzBounds.getMinimumY(), xyzBounds.getMaximumY()); + } else if (zDelta <= xDelta && zDelta <= yDelta) { + // Drill down in z + return !ztree.traverse(intersector, xyzBounds.getMinimumZ(), xyzBounds.getMaximumZ()); + } + return true; } @@ -153,6 +172,7 @@ class GeoComplexPolygon extends GeoBasePolygon { private static class Edge { public final GeoPoint startPoint; public final GeoPoint endPoint; + public final GeoPoint[] notablePoints; public final SidedPlane startPlane; public final SidedPlane endPlane; public final Plane plane; @@ -163,6 +183,7 @@ class GeoComplexPolygon extends GeoBasePolygon { public Edge(final PlanetModel pm, final GeoPoint startPoint, final GeoPoint endPoint) { this.startPoint = startPoint; this.endPoint = endPoint; + this.notablePoints = new GeoPoint[]{startPoint, endPoint}; this.plane = new Plane(startPoint, endPoint); this.startPlane = new SidedPlane(endPoint, plane, startPoint); this.endPlane = new SidedPlane(startPoint, plane, endPoint); @@ -193,10 +214,11 @@ class GeoComplexPolygon extends GeoBasePolygon { /** * Compare an edge. * @param edge is the edge to compare. - * @param value is the value to compare. + * @param minValue is the minimum value to compare (bottom of the range) + * @param maxValue is the maximum value to compare (top of the range) * @return -1 if "less" than this one, 0 if overlaps, or 1 if "greater". */ - public int compare(final Edge edge, final double value); + public int compare(final Edge edge, final double minValue, final double maxValue); } @@ -263,10 +285,10 @@ class GeoComplexPolygon extends GeoBasePolygon { } } - public boolean traverse(final EdgeIterator edgeIterator, final TraverseComparator edgeComparator, final double value) { + public boolean traverse(final EdgeIterator edgeIterator, final TraverseComparator edgeComparator, final double minValue, final double maxValue) { Node currentNode = this; while (currentNode != null) { - final int result = edgeComparator.compare(currentNode.edge, value); + final int result = edgeComparator.compare(currentNode.edge, minValue, maxValue); if (result < 0) { currentNode = lesser; } else if (result > 0) { @@ -298,11 +320,11 @@ class GeoComplexPolygon extends GeoBasePolygon { } } - public boolean traverse(final EdgeIterator edgeIterator, final double value) { + public boolean traverse(final EdgeIterator edgeIterator, final double minValue, final double maxValue) { if (rootNode == null) { return true; } - return rootNode.traverse(edgeIterator, this, value); + return rootNode.traverse(edgeIterator, this, minValue, maxValue); } @Override @@ -316,10 +338,10 @@ class GeoComplexPolygon extends GeoBasePolygon { } @Override - public int compare(final Edge edge, final double value) { - if (edge.planeBounds.getMinimumZ() > value) { + public int compare(final Edge edge, final double minValue, final double maxValue) { + if (edge.planeBounds.getMinimumZ() > maxValue) { return -1; - } else if (edge.planeBounds.getMaximumZ() < value) { + } else if (edge.planeBounds.getMaximumZ() < minValue) { return 1; } return 0; @@ -343,11 +365,11 @@ class GeoComplexPolygon extends GeoBasePolygon { } } - public boolean traverse(final EdgeIterator edgeIterator, final double value) { + public boolean traverse(final EdgeIterator edgeIterator, final double minValue, final double maxValue) { if (rootNode == null) { return true; } - return rootNode.traverse(edgeIterator, this, value); + return rootNode.traverse(edgeIterator, this, minValue, maxValue); } @Override @@ -361,10 +383,10 @@ class GeoComplexPolygon extends GeoBasePolygon { } @Override - public int compare(final Edge edge, final double value) { - if (edge.planeBounds.getMinimumY() > value) { + public int compare(final Edge edge, final double minValue, final double maxValue) { + if (edge.planeBounds.getMinimumY() > maxValue) { return -1; - } else if (edge.planeBounds.getMaximumY() < value) { + } else if (edge.planeBounds.getMaximumY() < minValue) { return 1; } return 0; @@ -388,11 +410,11 @@ class GeoComplexPolygon extends GeoBasePolygon { } } - public boolean traverse(final EdgeIterator edgeIterator, final double value) { + public boolean traverse(final EdgeIterator edgeIterator, final double minValue, final double maxValue) { if (rootNode == null) { return true; } - return rootNode.traverse(edgeIterator, this, value); + return rootNode.traverse(edgeIterator, this, minValue, maxValue); } @Override @@ -406,10 +428,10 @@ class GeoComplexPolygon extends GeoBasePolygon { } @Override - public int compare(final Edge edge, final double value) { - if (edge.planeBounds.getMinimumX() > value) { + public int compare(final Edge edge, final double minValue, final double maxValue) { + if (edge.planeBounds.getMinimumX() > maxValue) { return -1; - } else if (edge.planeBounds.getMaximumX() < value) { + } else if (edge.planeBounds.getMaximumX() < minValue) { return 1; } return 0; @@ -417,6 +439,27 @@ class GeoComplexPolygon extends GeoBasePolygon { } + /** Assess whether edge intersects the provided plane plus bounds. + */ + private class IntersectorEdgeIterator implements EdgeIterator { + + private final Plane plane; + private final GeoPoint[] notablePoints; + private final Membership[] bounds; + + public IntersectorEdgeIterator(final Plane plane, final GeoPoint[] notablePoints, final Membership... bounds) { + this.plane = plane; + this notablePoints = notablePoints; + this.bounds = bounds; + } + + @Override + public boolean matches(final Edge edge) { + return !plane.intersects(planetModel, edge.plane, notablePoints, edge.notablePoints, bounds, edge.startPlane, edge.endPlane); + } + + } + @Override public boolean equals(Object o) { // MHL @@ -426,6 +469,7 @@ class GeoComplexPolygon extends GeoBasePolygon { @Override public int hashCode() { // MHL + return 0; } @Override diff --git a/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/Plane.java b/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/Plane.java index f0df49d7b27..a205ca915df 100755 --- a/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/Plane.java +++ b/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/Plane.java @@ -614,11 +614,11 @@ public class Plane extends Vector { } /** - * Public version of findIntersections. + * Find the intersection points between two planes, given a set of bounds. * @param planetModel is the planet model. * @param q is the plane to intersect with. * @param bounds are the bounds to consider to determine legal intersection points. - * @return the set of legal intersection points. + * @return the set of legal intersection points, or null if the planes are numerically identical. */ public GeoPoint[] findIntersections(final PlanetModel planetModel, final Plane q, final Membership... bounds) { if (isNumericallyIdentical(q)) { @@ -626,6 +626,23 @@ public class Plane extends Vector { } return findIntersections(planetModel, q, bounds, NO_BOUNDS); } + + /** + * Find the points between two planes, where one plane crosses the other, given a set of bounds. + * Crossing is not just intersection; the planes cannot touch at just one point on the ellipsoid, + * but must cross at two. + * + * @param planetModel is the planet model. + * @param q is the plane to intersect with. + * @param bounds are the bounds to consider to determine legal intersection points. + * @return the set of legal crossing points, or null if the planes are numerically identical. + */ + public GeoPoint[] findCrossings(final PlanetModel planetModel, final Plane q, final Membership... bounds) { + if (isNumericallyIdentical(q)) { + return null; + } + return findCrossings(planetModel, q, bounds, NO_BOUNDS); + } /** * Find the intersection points between two planes, given a set of bounds. @@ -755,6 +772,124 @@ public class Plane extends Vector { } } + /** + * Find the points between two planes, where one plane crosses the other, given a set of bounds. + * Crossing is not just intersection; the planes cannot touch at just one point on the ellipsoid, + * but must cross at two. + * + * @param planetModel is the planet model to use in finding points. + * @param q is the plane to intersect with. + * @param bounds is the set of bounds. + * @param moreBounds is another set of bounds. + * @return the intersection point(s) on the ellipsoid, if there are any. + */ + protected GeoPoint[] findCrosses(final PlanetModel planetModel, final Plane q, final Membership[] bounds, final Membership[] moreBounds) { + // This code in this method is very similar to findIntersections(), but eliminates the cases where + // crossings are detected. + // Unnormalized, unchecked... + final Vector lineVector = new Vector(y * q.z - z * q.y, z * q.x - x * q.z, x * q.y - y * q.x); + if (Math.abs(lineVector.x) < MINIMUM_RESOLUTION && Math.abs(lineVector.y) < MINIMUM_RESOLUTION && Math.abs(lineVector.z) < MINIMUM_RESOLUTION) { + // Degenerate case: parallel planes + return NO_POINTS; + } + + // The line will have the equation: A t + A0 = x, B t + B0 = y, C t + C0 = z. + // We have A, B, and C. In order to come up with A0, B0, and C0, we need to find a point that is on both planes. + // To do this, we find the largest vector value (either x, y, or z), and look for a point that solves both plane equations + // simultaneous. For example, let's say that the vector is (0.5,0.5,1), and the two plane equations are: + // 0.7 x + 0.3 y + 0.1 z + 0.0 = 0 + // and + // 0.9 x - 0.1 y + 0.2 z + 4.0 = 0 + // Then we'd pick z = 0, so the equations to solve for x and y would be: + // 0.7 x + 0.3y = 0.0 + // 0.9 x - 0.1y = -4.0 + // ... which can readily be solved using standard linear algebra. Generally: + // Q0 x + R0 y = S0 + // Q1 x + R1 y = S1 + // ... can be solved by Cramer's rule: + // x = det(S0 R0 / S1 R1) / det(Q0 R0 / Q1 R1) + // y = det(Q0 S0 / Q1 S1) / det(Q0 R0 / Q1 R1) + // ... where det( a b / c d ) = ad - bc, so: + // x = (S0 * R1 - R0 * S1) / (Q0 * R1 - R0 * Q1) + // y = (Q0 * S1 - S0 * Q1) / (Q0 * R1 - R0 * Q1) + double x0; + double y0; + double z0; + // We try to maximize the determinant in the denominator + final double denomYZ = this.y * q.z - this.z * q.y; + final double denomXZ = this.x * q.z - this.z * q.x; + final double denomXY = this.x * q.y - this.y * q.x; + if (Math.abs(denomYZ) >= Math.abs(denomXZ) && Math.abs(denomYZ) >= Math.abs(denomXY)) { + // X is the biggest, so our point will have x0 = 0.0 + if (Math.abs(denomYZ) < MINIMUM_RESOLUTION_SQUARED) { + return NO_POINTS; + } + final double denom = 1.0 / denomYZ; + x0 = 0.0; + y0 = (-this.D * q.z - this.z * -q.D) * denom; + z0 = (this.y * -q.D + this.D * q.y) * denom; + } else if (Math.abs(denomXZ) >= Math.abs(denomXY) && Math.abs(denomXZ) >= Math.abs(denomYZ)) { + // Y is the biggest, so y0 = 0.0 + if (Math.abs(denomXZ) < MINIMUM_RESOLUTION_SQUARED) { + return NO_POINTS; + } + final double denom = 1.0 / denomXZ; + x0 = (-this.D * q.z - this.z * -q.D) * denom; + y0 = 0.0; + z0 = (this.x * -q.D + this.D * q.x) * denom; + } else { + // Z is the biggest, so Z0 = 0.0 + if (Math.abs(denomXY) < MINIMUM_RESOLUTION_SQUARED) { + return NO_POINTS; + } + final double denom = 1.0 / denomXY; + x0 = (-this.D * q.y - this.y * -q.D) * denom; + y0 = (this.x * -q.D + this.D * q.x) * denom; + z0 = 0.0; + } + + // Once an intersecting line is determined, the next step is to intersect that line with the ellipsoid, which + // will yield zero, one, or two points. + // The ellipsoid equation: 1,0 = x^2/a^2 + y^2/b^2 + z^2/c^2 + // 1.0 = (At+A0)^2/a^2 + (Bt+B0)^2/b^2 + (Ct+C0)^2/c^2 + // A^2 t^2 / a^2 + 2AA0t / a^2 + A0^2 / a^2 + B^2 t^2 / b^2 + 2BB0t / b^2 + B0^2 / b^2 + C^2 t^2 / c^2 + 2CC0t / c^2 + C0^2 / c^2 - 1,0 = 0.0 + // [A^2 / a^2 + B^2 / b^2 + C^2 / c^2] t^2 + [2AA0 / a^2 + 2BB0 / b^2 + 2CC0 / c^2] t + [A0^2 / a^2 + B0^2 / b^2 + C0^2 / c^2 - 1,0] = 0.0 + // Use the quadratic formula to determine t values and candidate point(s) + final double A = lineVector.x * lineVector.x * planetModel.inverseAbSquared + + lineVector.y * lineVector.y * planetModel.inverseAbSquared + + lineVector.z * lineVector.z * planetModel.inverseCSquared; + final double B = 2.0 * (lineVector.x * x0 * planetModel.inverseAbSquared + lineVector.y * y0 * planetModel.inverseAbSquared + lineVector.z * z0 * planetModel.inverseCSquared); + final double C = x0 * x0 * planetModel.inverseAbSquared + y0 * y0 * planetModel.inverseAbSquared + z0 * z0 * planetModel.inverseCSquared - 1.0; + + final double BsquaredMinus = B * B - 4.0 * A * C; + if (Math.abs(BsquaredMinus) < MINIMUM_RESOLUTION_SQUARED) { + // One point of intersection: cannot be a crossing. + return NO_POINTS; + } else if (BsquaredMinus > 0.0) { + final double inverse2A = 1.0 / (2.0 * A); + // Two solutions + final double sqrtTerm = Math.sqrt(BsquaredMinus); + final double t1 = (-B + sqrtTerm) * inverse2A; + final double t2 = (-B - sqrtTerm) * inverse2A; + GeoPoint point1 = new GeoPoint(lineVector.x * t1 + x0, lineVector.y * t1 + y0, lineVector.z * t1 + z0); + GeoPoint point2 = new GeoPoint(lineVector.x * t2 + x0, lineVector.y * t2 + y0, lineVector.z * t2 + z0); + //verifyPoint(planetModel, point1, q); + //verifyPoint(planetModel, point2, q); + //System.err.println(" "+point1+" and "+point2); + if (point1.isWithin(bounds, moreBounds)) { + if (point2.isWithin(bounds, moreBounds)) + return new GeoPoint[]{point1, point2}; + return new GeoPoint[]{point1}; + } + if (point2.isWithin(bounds, moreBounds)) + return new GeoPoint[]{point2}; + return NO_POINTS; + } else { + // No solutions. + return NO_POINTS; + } + } + /* protected void verifyPoint(final PlanetModel planetModel, final GeoPoint point, final Plane q) { if (!evaluateIsZero(point))