You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@lucene.apache.org by mi...@apache.org on 2015/09/08 18:29:23 UTC
svn commit: r1701837 - in /lucene/dev/branches/branch_5x: ./ lucene/
lucene/spatial3d/ lucene/spatial3d/src/java/org/apache/lucene/geo3d/
lucene/spatial3d/src/test/org/apache/lucene/bkdtree3d/
lucene/spatial3d/src/test/org/apache/lucene/geo3d/
Author: mikemccand
Date: Tue Sep 8 16:29:22 2015
New Revision: 1701837
URL: http://svn.apache.org/r1701837
Log:
LUCENE-6776: Fix geo3d math to handle randomly squashed planets
Modified:
lucene/dev/branches/branch_5x/ (props changed)
lucene/dev/branches/branch_5x/lucene/ (props changed)
lucene/dev/branches/branch_5x/lucene/CHANGES.txt (contents, props changed)
lucene/dev/branches/branch_5x/lucene/spatial3d/ (props changed)
lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/GeoPath.java
lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/Plane.java
lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/PlanetModel.java
lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/Vector.java
lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYZBounds.java
lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYZSolid.java
lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYdZSolid.java
lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XdYZSolid.java
lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/dXYZSolid.java
lucene/dev/branches/branch_5x/lucene/spatial3d/src/test/org/apache/lucene/bkdtree3d/TestGeo3DPointField.java
lucene/dev/branches/branch_5x/lucene/spatial3d/src/test/org/apache/lucene/geo3d/GeoCircleTest.java
lucene/dev/branches/branch_5x/lucene/spatial3d/src/test/org/apache/lucene/geo3d/GeoPathTest.java
Modified: lucene/dev/branches/branch_5x/lucene/CHANGES.txt
URL: http://svn.apache.org/viewvc/lucene/dev/branches/branch_5x/lucene/CHANGES.txt?rev=1701837&r1=1701836&r2=1701837&view=diff
==============================================================================
--- lucene/dev/branches/branch_5x/lucene/CHANGES.txt (original)
+++ lucene/dev/branches/branch_5x/lucene/CHANGES.txt Tue Sep 8 16:29:22 2015
@@ -71,6 +71,9 @@ Bug Fixes
* LUCENE-6783: Removed side effects from FuzzyLikeThisQuery.rewrite.
(Adrien Grand)
+* LUCENE-6776: Fix geo3d math to handle randomly squashed planet
+ models (Karl Wright via Mike McCandless)
+
Other
* LUCENE-6174: Improve "ant eclipse" to select right JRE for building.
Modified: lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/GeoPath.java
URL: http://svn.apache.org/viewvc/lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/GeoPath.java?rev=1701837&r1=1701836&r2=1701837&view=diff
==============================================================================
--- lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/GeoPath.java (original)
+++ lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/GeoPath.java Tue Sep 8 16:29:22 2015
@@ -475,7 +475,14 @@ public class GeoPath extends GeoBaseDist
public boolean isWithin(final Vector point) {
if (circlePlane == null)
return false;
- return circlePlane.isWithin(point);
+ if (!circlePlane.isWithin(point))
+ return false;
+ for (final Membership m : cutoffPlanes) {
+ if (!m.isWithin(point)) {
+ return false;
+ }
+ }
+ return true;
}
/** Check if point is within this endpoint.
@@ -487,7 +494,14 @@ public class GeoPath extends GeoBaseDist
public boolean isWithin(final double x, final double y, final double z) {
if (circlePlane == null)
return false;
- return circlePlane.isWithin(x, y, z);
+ if (!circlePlane.isWithin(x, y, z))
+ return false;
+ for (final Membership m : cutoffPlanes) {
+ if (!m.isWithin(x,y,z)) {
+ return false;
+ }
+ }
+ return true;
}
/** Compute interior path distance.
Modified: lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/Plane.java
URL: http://svn.apache.org/viewvc/lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/Plane.java?rev=1701837&r1=1701836&r2=1701837&view=diff
==============================================================================
--- lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/Plane.java (original)
+++ lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/Plane.java Tue Sep 8 16:29:22 2015
@@ -551,7 +551,6 @@ public class Plane extends Vector {
final double startMagnitude = Math.sqrt(modifiedStart.x * modifiedStart.x + modifiedStart.y * modifiedStart.y);
double delta;
- double beginAngle;
double newEndAngle = endAngle;
while (newEndAngle < startAngle) {
@@ -560,14 +559,12 @@ public class Plane extends Vector {
if (newEndAngle - startAngle <= Math.PI) {
delta = newEndAngle - startAngle;
- beginAngle = startAngle;
} else {
double newStartAngle = startAngle;
while (newStartAngle < endAngle) {
newStartAngle += Math.PI * 2.0;
}
delta = newStartAngle - endAngle;
- beginAngle = endAngle;
}
final GeoPoint[] returnValues = new GeoPoint[proportions.length];
@@ -789,30 +786,7 @@ public class Plane extends Vector {
final double B = this.y;
final double C = this.z;
- // For the X and Y values, we need a D value, which is the AVERAGE D value
- // for two planes that pass through the two Z points determined here, for the axis in question.
- final GeoPoint[] zPoints;
- if (!boundsInfo.isSmallestMinX(planetModel) || !boundsInfo.isLargestMaxX(planetModel) ||
- !boundsInfo.isSmallestMinY(planetModel) || !boundsInfo.isLargestMaxY(planetModel)) {
- if ((Math.abs(A) >= MINIMUM_RESOLUTION || Math.abs(B) >= MINIMUM_RESOLUTION)) {
- // We need unconstrained values in order to compute D
- zPoints = findIntersections(planetModel, constructNormalizedZPlane(A,B), NO_BOUNDS, NO_BOUNDS);
- if (zPoints.length == 0) {
- // No intersections, so plane never intersects world.
- //System.err.println(" plane never intersects world");
- return;
- }
- //for (final GeoPoint p : zPoints) {
- // System.err.println("zPoint: "+p);
- //}
- } else {
- zPoints = null;
- }
- } else {
- zPoints = null;
- }
-
- // Do Z.
+ // Do Z. This can be done simply because it is symmetrical.
if (!boundsInfo.isSmallestMinZ(planetModel) || !boundsInfo.isLargestMaxZ(planetModel)) {
//System.err.println(" computing Z bound");
// Compute Z bounds for this arc
@@ -839,65 +813,339 @@ public class Plane extends Vector {
boundsInfo.addZValue(points[0]);
}
}
-
+
+ // First, compute common subexpressions
+ final double k = 1.0 / ((x*x + y*y)*planetModel.ab*planetModel.ab + z*z*planetModel.c*planetModel.c);
+ final double abSquared = planetModel.ab * planetModel.ab;
+ final double cSquared = planetModel.c * planetModel.c;
+ final double ASquared = A * A;
+ final double BSquared = B * B;
+ final double CSquared = C * C;
+
+ final double r = 2.0*D*k;
+ final double rSquared = r * r;
+
if (!boundsInfo.isSmallestMinX(planetModel) || !boundsInfo.isLargestMaxX(planetModel)) {
+ // For min/max x, we need to use lagrange multipliers.
+ //
+ // For this, we need grad(F(x,y,z)) = (dF/dx, dF/dy, dF/dz).
+ //
+ // Minimize and maximize f(x,y,z) = x, with respect to g(x,y,z) = Ax + By + Cz - D and h(x,y,z) = x^2/ab^2 + y^2/ab^2 + z^2/c^2 - 1
+ //
+ // grad(f(x,y,z)) = (1,0,0)
+ // grad(g(x,y,z)) = (A,B,C)
+ // grad(h(x,y,z)) = (2x/ab^2,2y/ab^2,2z/c^2)
+ //
+ // Equations we need to simultaneously solve:
+ //
+ // grad(f(x,y,z)) = l * grad(g(x,y,z)) + m * grad(h(x,y,z))
+ // g(x,y,z) = 0
+ // h(x,y,z) = 0
+ //
+ // Equations:
+ // 1 = l*A + m*2x/ab^2
+ // 0 = l*B + m*2y/ab^2
+ // 0 = l*C + m*2z/c^2
+ // Ax + By + Cz + D = 0
+ // x^2/ab^2 + y^2/ab^2 + z^2/c^2 - 1 = 0
+ //
+ // Solve for x,y,z in terms of (l, m):
+ //
+ // x = ((1 - l*A) * ab^2 ) / (2 * m)
+ // y = (-l*B * ab^2) / ( 2 * m)
+ // z = (-l*C * c^2)/ (2 * m)
+ //
+ // Two equations, two unknowns:
+ //
+ // A * (((1 - l*A) * ab^2 ) / (2 * m)) + B * ((-l*B * ab^2) / ( 2 * m)) + C * ((-l*C * c^2)/ (2 * m)) + D = 0
+ //
+ // and
+ //
+ // (((1 - l*A) * ab^2 ) / (2 * m))^2/ab^2 + ((-l*B * ab^2) / ( 2 * m))^2/ab^2 + ((-l*C * c^2)/ (2 * m))^2/c^2 - 1 = 0
+ //
+ // Simple: solve for l and m, then find x from it.
+ //
+ // (a) Use first equation to find l in terms of m.
+ //
+ // A * (((1 - l*A) * ab^2 ) / (2 * m)) + B * ((-l*B * ab^2) / ( 2 * m)) + C * ((-l*C * c^2)/ (2 * m)) + D = 0
+ // A * ((1 - l*A) * ab^2 ) + B * (-l*B * ab^2) + C * (-l*C * c^2) + D * 2 * m = 0
+ // A * ab^2 - l*A^2* ab^2 - B^2 * l * ab^2 - C^2 * l * c^2 + D * 2 * m = 0
+ // - l *(A^2* ab^2 + B^2 * ab^2 + C^2 * c^2) + (A * ab^2 + D * 2 * m) = 0
+ // l = (A * ab^2 + D * 2 * m) / (A^2* ab^2 + B^2 * ab^2 + C^2 * c^2)
+ // l = A * ab^2 / (A^2* ab^2 + B^2 * ab^2 + C^2 * c^2) + m * 2 * D / (A^2* ab^2 + B^2 * ab^2 + C^2 * c^2)
+ //
+ // For convenience:
+ //
+ // k = 1.0 / (A^2* ab^2 + B^2 * ab^2 + C^2 * c^2)
+ //
+ // Then:
+ //
+ // l = A * ab^2 * k + m * 2 * D * k
+ // l = k * (A*ab^2 + m*2*D)
+ //
+ // For further convenience:
+ //
+ // q = A*ab^2*k
+ // r = 2*D*k
+ //
+ // l = (r*m + q)
+ // l^2 = (r^2 * m^2 + 2*r*m*q + q^2)
+ //
+ // (b) Simplify the second equation before substitution
+ //
+ // (((1 - l*A) * ab^2 ) / (2 * m))^2/ab^2 + ((-l*B * ab^2) / ( 2 * m))^2/ab^2 + ((-l*C * c^2)/ (2 * m))^2/c^2 - 1 = 0
+ // ((1 - l*A) * ab^2 )^2/ab^2 + (-l*B * ab^2)^2/ab^2 + (-l*C * c^2)^2/c^2 = 4 * m^2
+ // (1 - l*A)^2 * ab^2 + (-l*B)^2 * ab^2 + (-l*C)^2 * c^2 = 4 * m^2
+ // (1 - 2*l*A + l^2*A^2) * ab^2 + l^2*B^2 * ab^2 + l^2*C^2 * c^2 = 4 * m^2
+ // ab^2 - 2*A*ab^2*l + A^2*ab^2*l^2 + B^2*ab^2*l^2 + C^2*c^2*l^2 - 4*m^2 = 0
+ //
+ // (c) Substitute for l, l^2
+ //
+ // ab^2 - 2*A*ab^2*(r*m + q) + A^2*ab^2*(r^2 * m^2 + 2*r*m*q + q^2) + B^2*ab^2*(r^2 * m^2 + 2*r*m*q + q^2) + C^2*c^2*(r^2 * m^2 + 2*r*m*q + q^2) - 4*m^2 = 0
+ // ab^2 - 2*A*ab^2*r*m - 2*A*ab^2*q + A^2*ab^2*r^2*m^2 + 2*A^2*ab^2*r*q*m +
+ // A^2*ab^2*q^2 + B^2*ab^2*r^2*m^2 + 2*B^2*ab^2*r*q*m + B^2*ab^2*q^2 + C^2*c^2*r^2*m^2 + 2*C^2*c^2*r*q*m + C^2*c^2*q^2 - 4*m^2 = 0
+ //
+ // (d) Group
+ //
+ // m^2 * [A^2*ab^2*r^2 + B^2*ab^2*r^2 + C^2*c^2*r^2 - 4] +
+ // m * [- 2*A*ab^2*r + 2*A^2*ab^2*r*q + 2*B^2*ab^2*r*q + 2*C^2*c^2*r*q] +
+ // [ab^2 - 2*A*ab^2*q + A^2*ab^2*q^2 + B^2*ab^2*q^2 + C^2*c^2*q^2] = 0
+
//System.err.println(" computing X bound");
- if ((Math.abs(B) >= MINIMUM_RESOLUTION || Math.abs(C) >= MINIMUM_RESOLUTION)) {
- // NOT a degenerate case. Compute D.
- //System.err.println(" not degenerate; B="+B+"; C="+C);
- final Plane originPlane = constructNormalizedXPlane(B,C,0.0);
- double DValue = 0.0;
- if (zPoints != null) {
- for (final GeoPoint p : zPoints) {
- final double zValue = originPlane.evaluate(p);
- //System.err.println(" originPlane.evaluate(zpoint)="+zValue);
- DValue += zValue;
- }
- DValue /= (double)zPoints.length;
- }
- final Plane normalizedXPlane = constructNormalizedXPlane(B,C,-DValue);
- final GeoPoint[] points = findIntersections(planetModel, normalizedXPlane, bounds, NO_BOUNDS);
- for (final GeoPoint point : points) {
- assert planetModel.pointOnSurface(point);
- //System.err.println(" Point = "+point+"; this.evaluate(point)="+this.evaluate(point)+"; normalizedXPlane.evaluate(point)="+normalizedXPlane.evaluate(point));
- addPoint(boundsInfo, bounds, point);
+
+ // Useful subexpressions for this bound
+ final double q = A*abSquared*k;
+ final double qSquared = q * q;
+
+ // Quadratic equation
+ final double a = ASquared*abSquared*rSquared + BSquared*abSquared*rSquared + CSquared*cSquared*rSquared - 4.0;
+ final double b = - 2.0*A*abSquared*r + 2.0*ASquared*abSquared*r*q + 2.0*BSquared*abSquared*r*q + 2.0*CSquared*cSquared*r*q;
+ final double c = abSquared - 2.0*A*abSquared*q + ASquared*abSquared*qSquared + BSquared*abSquared*qSquared + CSquared*cSquared*qSquared;
+
+ if (Math.abs(a) >= MINIMUM_RESOLUTION_SQUARED) {
+ final double sqrtTerm = b*b - 4.0*a*c;
+ if (Math.abs(sqrtTerm) < MINIMUM_RESOLUTION_SQUARED) {
+ // One solution
+ final double m = -b / (2.0 * a);
+ final double l = r * m + q;
+ // x = ((1 - l*A) * ab^2 ) / (2 * m)
+ // y = (-l*B * ab^2) / ( 2 * m)
+ // z = (-l*C * c^2)/ (2 * m)
+ final double denom0 = 0.5 / m;
+ final GeoPoint thePoint = new GeoPoint((1.0-l*A) * abSquared * denom0, -l*B * abSquared * denom0, -l*C * cSquared * denom0);
+ //Math is not quite accurate enough for this
+ //assert planetModel.pointOnSurface(thePoint): "Point: "+thePoint+"; Planetmodel="+planetModel+"; A="+A+" B="+B+" C="+C+" D="+D+" planetfcn="+
+ // (thePoint.x*thePoint.x*planetModel.inverseAb*planetModel.inverseAb + thePoint.y*thePoint.y*planetModel.inverseAb*planetModel.inverseAb + thePoint.z*thePoint.z*planetModel.inverseC*planetModel.inverseC);
+ //assert evaluateIsZero(thePoint): "Evaluation of point: "+evaluate(thePoint);
+ addPoint(boundsInfo, bounds, thePoint);
+ } else if (sqrtTerm > 0.0) {
+ // Two solutions
+ final double sqrtResult = Math.sqrt(sqrtTerm);
+ final double commonDenom = 0.5/a;
+ final double m1 = (-b + sqrtResult) * commonDenom;
+ assert Math.abs(a * m1 * m1 + b * m1 + c) < MINIMUM_RESOLUTION;
+ final double m2 = (-b - sqrtResult) * commonDenom;
+ assert Math.abs(a * m2 * m2 + b * m2 + c) < MINIMUM_RESOLUTION;
+ final double l1 = r * m1 + q;
+ final double l2 = r * m2 + q;
+ // x = ((1 - l*A) * ab^2 ) / (2 * m)
+ // y = (-l*B * ab^2) / ( 2 * m)
+ // z = (-l*C * c^2)/ (2 * m)
+ final double denom1 = 0.5 / m1;
+ final double denom2 = 0.5 / m2;
+ final GeoPoint thePoint1 = new GeoPoint((1.0-l1*A) * abSquared * denom1, -l1*B * abSquared * denom1, -l1*C * cSquared * denom1);
+ final GeoPoint thePoint2 = new GeoPoint((1.0-l2*A) * abSquared * denom2, -l2*B * abSquared * denom2, -l2*C * cSquared * denom2);
+ //Math is not quite accurate enough for this
+ //assert planetModel.pointOnSurface(thePoint1): "Point1: "+thePoint1+"; Planetmodel="+planetModel+"; A="+A+" B="+B+" C="+C+" D="+D+" planetfcn="+
+ // (thePoint1.x*thePoint1.x*planetModel.inverseAb*planetModel.inverseAb + thePoint1.y*thePoint1.y*planetModel.inverseAb*planetModel.inverseAb + thePoint1.z*thePoint1.z*planetModel.inverseC*planetModel.inverseC);
+ //assert planetModel.pointOnSurface(thePoint2): "Point1: "+thePoint2+"; Planetmodel="+planetModel+"; A="+A+" B="+B+" C="+C+" D="+D+" planetfcn="+
+ // (thePoint2.x*thePoint2.x*planetModel.inverseAb*planetModel.inverseAb + thePoint2.y*thePoint2.y*planetModel.inverseAb*planetModel.inverseAb + thePoint2.z*thePoint2.z*planetModel.inverseC*planetModel.inverseC);
+ //assert evaluateIsZero(thePoint1): "Evaluation of point1: "+evaluate(thePoint1);
+ //assert evaluateIsZero(thePoint2): "Evaluation of point2: "+evaluate(thePoint2);
+ addPoint(boundsInfo, bounds, thePoint1);
+ addPoint(boundsInfo, bounds, thePoint2);
+ } else {
+ // No solutions
}
+ } else if (Math.abs(b) > MINIMUM_RESOLUTION_SQUARED) {
+ //System.err.println("Not x quadratic");
+ // a = 0, so m = - c / b
+ final double m = -c / b;
+ final double l = r * m + q;
+ // x = ((1 - l*A) * ab^2 ) / (2 * m)
+ // y = (-l*B * ab^2) / ( 2 * m)
+ // z = (-l*C * c^2)/ (2 * m)
+ final double denom0 = 0.5 / m;
+ final GeoPoint thePoint = new GeoPoint((1.0-l*A) * abSquared * denom0, -l*B * abSquared * denom0, -l*C * cSquared * denom0);
+ //Math is not quite accurate enough for this
+ //assert planetModel.pointOnSurface(thePoint): "Point: "+thePoint+"; Planetmodel="+planetModel+"; A="+A+" B="+B+" C="+C+" D="+D+" planetfcn="+
+ // (thePoint.x*thePoint.x*planetModel.inverseAb*planetModel.inverseAb + thePoint.y*thePoint.y*planetModel.inverseAb*planetModel.inverseAb + thePoint.z*thePoint.z*planetModel.inverseC*planetModel.inverseC);
+ //assert evaluateIsZero(thePoint): "Evaluation of point: "+evaluate(thePoint);
+ addPoint(boundsInfo, bounds, thePoint);
} else {
- // Since b==c==0, any plane including the X axis suffices.
- //System.err.println(" Perpendicular to x");
- final GeoPoint[] points = findIntersections(planetModel, normalZPlane, NO_BOUNDS, NO_BOUNDS);
- boundsInfo.addXValue(points[0]);
+ // Something went very wrong; a = b = 0
}
}
// Do Y
if (!boundsInfo.isSmallestMinY(planetModel) || !boundsInfo.isLargestMaxY(planetModel)) {
+ // For min/max x, we need to use lagrange multipliers.
+ //
+ // For this, we need grad(F(x,y,z)) = (dF/dx, dF/dy, dF/dz).
+ //
+ // Minimize and maximize f(x,y,z) = y, with respect to g(x,y,z) = Ax + By + Cz - D and h(x,y,z) = x^2/ab^2 + y^2/ab^2 + z^2/c^2 - 1
+ //
+ // grad(f(x,y,z)) = (0,1,0)
+ // grad(g(x,y,z)) = (A,B,C)
+ // grad(h(x,y,z)) = (2x/ab^2,2y/ab^2,2z/c^2)
+ //
+ // Equations we need to simultaneously solve:
+ //
+ // grad(f(x,y,z)) = l * grad(g(x,y,z)) + m * grad(h(x,y,z))
+ // g(x,y,z) = 0
+ // h(x,y,z) = 0
+ //
+ // Equations:
+ // 0 = l*A + m*2x/ab^2
+ // 1 = l*B + m*2y/ab^2
+ // 0 = l*C + m*2z/c^2
+ // Ax + By + Cz + D = 0
+ // x^2/ab^2 + y^2/ab^2 + z^2/c^2 - 1 = 0
+ //
+ // Solve for x,y,z in terms of (l, m):
+ //
+ // x = (-l*A * ab^2 ) / (2 * m)
+ // y = ((1 - l*B) * ab^2) / ( 2 * m)
+ // z = (-l*C * c^2)/ (2 * m)
+ //
+ // Two equations, two unknowns:
+ //
+ // A * ((-l*A * ab^2 ) / (2 * m)) + B * (((1 - l*B) * ab^2) / ( 2 * m)) + C * ((-l*C * c^2)/ (2 * m)) + D = 0
+ //
+ // and
+ //
+ // ((-l*A * ab^2 ) / (2 * m))^2/ab^2 + (((1 - l*B) * ab^2) / ( 2 * m))^2/ab^2 + ((-l*C * c^2)/ (2 * m))^2/c^2 - 1 = 0
+ //
+ // Simple: solve for l and m, then find y from it.
+ //
+ // (a) Use first equation to find l in terms of m.
+ //
+ // A * ((-l*A * ab^2 ) / (2 * m)) + B * (((1 - l*B) * ab^2) / ( 2 * m)) + C * ((-l*C * c^2)/ (2 * m)) + D = 0
+ // A * (-l*A * ab^2 ) + B * ((1-l*B) * ab^2) + C * (-l*C * c^2) + D * 2 * m = 0
+ // -A^2*l*ab^2 + B*ab^2 - l*B^2*ab^2 - C^2*l*c^2 + D*2*m = 0
+ // - l *(A^2* ab^2 + B^2 * ab^2 + C^2 * c^2) + (B * ab^2 + D * 2 * m) = 0
+ // l = (B * ab^2 + D * 2 * m) / (A^2* ab^2 + B^2 * ab^2 + C^2 * c^2)
+ // l = B * ab^2 / (A^2* ab^2 + B^2 * ab^2 + C^2 * c^2) + m * 2 * D / (A^2* ab^2 + B^2 * ab^2 + C^2 * c^2)
+ //
+ // For convenience:
+ //
+ // k = 1.0 / (A^2* ab^2 + B^2 * ab^2 + C^2 * c^2)
+ //
+ // Then:
+ //
+ // l = B * ab^2 * k + m * 2 * D * k
+ // l = k * (B*ab^2 + m*2*D)
+ //
+ // For further convenience:
+ //
+ // q = B*ab^2*k
+ // r = 2*D*k
+ //
+ // l = (r*m + q)
+ // l^2 = (r^2 * m^2 + 2*r*m*q + q^2)
+ //
+ // (b) Simplify the second equation before substitution
+ //
+ // ((-l*A * ab^2 ) / (2 * m))^2/ab^2 + (((1 - l*B) * ab^2) / ( 2 * m))^2/ab^2 + ((-l*C * c^2)/ (2 * m))^2/c^2 - 1 = 0
+ // (-l*A * ab^2 )^2/ab^2 + ((1 - l*B) * ab^2)^2/ab^2 + (-l*C * c^2)^2/c^2 = 4 * m^2
+ // (-l*A)^2 * ab^2 + (1 - l*B)^2 * ab^2 + (-l*C)^2 * c^2 = 4 * m^2
+ // l^2*A^2 * ab^2 + (1 - 2*l*B + l^2*B^2) * ab^2 + l^2*C^2 * c^2 = 4 * m^2
+ // A^2*ab^2*l^2 + ab^2 - 2*B*ab^2*l + B^2*ab^2*l^2 + C^2*c^2*l^2 - 4*m^2 = 0
+ //
+ // (c) Substitute for l, l^2
+ //
+ // A^2*ab^2*(r^2 * m^2 + 2*r*m*q + q^2) + ab^2 - 2*B*ab^2*(r*m + q) + B^2*ab^2*(r^2 * m^2 + 2*r*m*q + q^2) + C^2*c^2*(r^2 * m^2 + 2*r*m*q + q^2) - 4*m^2 = 0
+ // A^2*ab^2*r^2*m^2 + 2*A^2*ab^2*r*q*m + A^2*ab^2*q^2 + ab^2 - 2*B*ab^2*r*m - 2*B*ab^2*q + B^2*ab^2*r^2*m^2 +
+ // 2*B^2*ab^2*r*q*m + B^2*ab^2*q^2 + C^2*c^2*r^2*m^2 + 2*C^2*c^2*r*q*m + C^2*c^2*q^2 - 4*m^2 = 0
+ //
+ // (d) Group
+ //
+ // m^2 * [A^2*ab^2*r^2 + B^2*ab^2*r^2 + C^2*c^2*r^2 - 4] +
+ // m * [2*A^2*ab^2*r*q - 2*B*ab^2*r + 2*B^2*ab^2*r*q + 2*C^2*c^2*r*q] +
+ // [A^2*ab^2*q^2 + ab^2 - 2*B*ab^2*q + B^2*ab^2*q^2 + C^2*c^2*q^2] = 0
+
//System.err.println(" computing Y bound");
- if ((Math.abs(A) >= MINIMUM_RESOLUTION || Math.abs(C) >= MINIMUM_RESOLUTION)) {
- // NOT a degenerate case. Compute D.
- //System.err.println(" not degenerate");
- final Plane originPlane = constructNormalizedYPlane(A,C,0.0);
- double DValue = 0.0;
- if (zPoints != null) {
- for (final GeoPoint p : zPoints) {
- DValue += originPlane.evaluate(p);
- }
- DValue /= (double)zPoints.length;
- }
- final Plane normalizedYPlane = constructNormalizedYPlane(A,C,-DValue);
- final GeoPoint[] points = findIntersections(planetModel, normalizedYPlane, bounds, NO_BOUNDS);
- for (final GeoPoint point : points) {
- assert planetModel.pointOnSurface(point);
- //System.err.println(" Point = "+point+"; this.evaluate(point)="+this.evaluate(point)+"; normalizedYPlane.evaluate(point)="+normalizedYPlane.evaluate(point));
- addPoint(boundsInfo, bounds, point);
+
+ // Useful subexpressions for this bound
+ final double q = B*abSquared*k;
+ final double qSquared = q * q;
+
+ // Quadratic equation
+ final double a = ASquared*abSquared*rSquared + BSquared*abSquared*rSquared + CSquared*cSquared*rSquared - 4.0;
+ final double b = 2.0*ASquared*abSquared*r*q - 2.0*B*abSquared*r + 2.0*BSquared*abSquared*r*q + 2.0*CSquared*cSquared*r*q;
+ final double c = ASquared*abSquared*qSquared + abSquared - 2.0*B*abSquared*q + BSquared*abSquared*qSquared + CSquared*cSquared*qSquared;
+
+ if (Math.abs(a) >= MINIMUM_RESOLUTION_SQUARED) {
+ final double sqrtTerm = b*b - 4.0*a*c;
+ if (Math.abs(sqrtTerm) < MINIMUM_RESOLUTION_SQUARED) {
+ // One solution
+ final double m = -b / (2.0 * a);
+ final double l = r * m + q;
+ // x = (-l*A * ab^2 ) / (2 * m)
+ // y = ((1.0-l*B) * ab^2) / ( 2 * m)
+ // z = (-l*C * c^2)/ (2 * m)
+ final double denom0 = 0.5 / m;
+ final GeoPoint thePoint = new GeoPoint(-l*A * abSquared * denom0, (1.0-l*B) * abSquared * denom0, -l*C * cSquared * denom0);
+ //Math is not quite accurate enough for this
+ //assert planetModel.pointOnSurface(thePoint): "Point: "+thePoint+"; Planetmodel="+planetModel+"; A="+A+" B="+B+" C="+C+" D="+D+" planetfcn="+
+ // (thePoint1.x*thePoint.x*planetModel.inverseAb*planetModel.inverseAb + thePoint.y*thePoint.y*planetModel.inverseAb*planetModel.inverseAb + thePoint.z*thePoint.z*planetModel.inverseC*planetModel.inverseC);
+ //assert evaluateIsZero(thePoint): "Evaluation of point: "+evaluate(thePoint);
+ addPoint(boundsInfo, bounds, thePoint);
+ } else if (sqrtTerm > 0.0) {
+ // Two solutions
+ final double sqrtResult = Math.sqrt(sqrtTerm);
+ final double commonDenom = 0.5/a;
+ final double m1 = (-b + sqrtResult) * commonDenom;
+ assert Math.abs(a * m1 * m1 + b * m1 + c) < MINIMUM_RESOLUTION;
+ final double m2 = (-b - sqrtResult) * commonDenom;
+ assert Math.abs(a * m2 * m2 + b * m2 + c) < MINIMUM_RESOLUTION;
+ final double l1 = r * m1 + q;
+ final double l2 = r * m2 + q;
+ // x = (-l*A * ab^2 ) / (2 * m)
+ // y = ((1.0-l*B) * ab^2) / ( 2 * m)
+ // z = (-l*C * c^2)/ (2 * m)
+ final double denom1 = 0.5 / m1;
+ final double denom2 = 0.5 / m2;
+ final GeoPoint thePoint1 = new GeoPoint(-l1*A * abSquared * denom1, (1.0-l1*B) * abSquared * denom1, -l1*C * cSquared * denom1);
+ final GeoPoint thePoint2 = new GeoPoint(-l2*A * abSquared * denom2, (1.0-l2*B) * abSquared * denom2, -l2*C * cSquared * denom2);
+ //Math is not quite accurate enough for this
+ //assert planetModel.pointOnSurface(thePoint1): "Point1: "+thePoint1+"; Planetmodel="+planetModel+"; A="+A+" B="+B+" C="+C+" D="+D+" planetfcn="+
+ // (thePoint1.x*thePoint1.x*planetModel.inverseAb*planetModel.inverseAb + thePoint1.y*thePoint1.y*planetModel.inverseAb*planetModel.inverseAb + thePoint1.z*thePoint1.z*planetModel.inverseC*planetModel.inverseC);
+ //assert planetModel.pointOnSurface(thePoint2): "Point2: "+thePoint2+"; Planetmodel="+planetModel+"; A="+A+" B="+B+" C="+C+" D="+D+" planetfcn="+
+ // (thePoint2.x*thePoint2.x*planetModel.inverseAb*planetModel.inverseAb + thePoint2.y*thePoint2.y*planetModel.inverseAb*planetModel.inverseAb + thePoint2.z*thePoint2.z*planetModel.inverseC*planetModel.inverseC);
+ //assert evaluateIsZero(thePoint1): "Evaluation of point1: "+evaluate(thePoint1);
+ //assert evaluateIsZero(thePoint2): "Evaluation of point2: "+evaluate(thePoint2);
+ addPoint(boundsInfo, bounds, thePoint1);
+ addPoint(boundsInfo, bounds, thePoint2);
+ } else {
+ // No solutions
}
+ } else if (Math.abs(b) > MINIMUM_RESOLUTION_SQUARED) {
+ // a = 0, so m = - c / b
+ final double m = -c / b;
+ final double l = r * m + q;
+ // x = ( -l*A * ab^2 ) / (2 * m)
+ // y = ((1-l*B) * ab^2) / ( 2 * m)
+ // z = (-l*C * c^2)/ (2 * m)
+ final double denom0 = 0.5 / m;
+ final GeoPoint thePoint = new GeoPoint(-l*A * abSquared * denom0, (1.0-l*B) * abSquared * denom0, -l*C * cSquared * denom0);
+ //Math is not quite accurate enough for this
+ //assert planetModel.pointOnSurface(thePoint): "Point: "+thePoint+"; Planetmodel="+planetModel+"; A="+A+" B="+B+" C="+C+" D="+D+" planetfcn="+
+ // (thePoint.x*thePoint.x*planetModel.inverseAb*planetModel.inverseAb + thePoint.y*thePoint.y*planetModel.inverseAb*planetModel.inverseAb + thePoint.z*thePoint.z*planetModel.inverseC*planetModel.inverseC);
+ //assert evaluateIsZero(thePoint): "Evaluation of point: "+evaluate(thePoint);
+ addPoint(boundsInfo, bounds, thePoint);
} else {
- // Since a==c==0, any plane including the Y axis suffices.
- // It doesn't matter that we may discard the point due to bounds, because if there are bounds, we'll have endpoints
- // that will be tallied separately.
- //System.err.println(" Perpendicular to y");
- final GeoPoint[] points = findIntersections(planetModel, normalXPlane, NO_BOUNDS, NO_BOUNDS);
- boundsInfo.addYValue(points[0]);
+ // Something went very wrong; a = b = 0
}
}
}
Modified: lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/PlanetModel.java
URL: http://svn.apache.org/viewvc/lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/PlanetModel.java?rev=1701837&r1=1701836&r2=1701837&view=diff
==============================================================================
--- lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/PlanetModel.java (original)
+++ lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/PlanetModel.java Tue Sep 8 16:29:22 2015
@@ -167,7 +167,7 @@ public class PlanetModel {
public boolean pointOnSurface(final double x, final double y, final double z) {
// Equation of planet surface is:
// x^2 / a^2 + y^2 / b^2 + z^2 / c^2 - 1 = 0
- return Math.abs((x * x + y * y) * inverseAb * inverseAb + z * z * inverseC * inverseC - 1.0) < Vector.MINIMUM_RESOLUTION;
+ return Math.abs(x * x * inverseAb * inverseAb + y * y * inverseAb * inverseAb + z * z * inverseC * inverseC - 1.0) < Vector.MINIMUM_RESOLUTION;
}
/** Check if point is outside surface.
Modified: lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/Vector.java
URL: http://svn.apache.org/viewvc/lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/Vector.java?rev=1701837&r1=1701836&r2=1701837&view=diff
==============================================================================
--- lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/Vector.java (original)
+++ lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/Vector.java Tue Sep 8 16:29:22 2015
@@ -28,7 +28,7 @@ public class Vector {
* Values that are all considered to be essentially zero have a magnitude
* less than this.
*/
- public static final double MINIMUM_RESOLUTION = 5.0e-13;
+ public static final double MINIMUM_RESOLUTION = 1.0e-12;
/**
* For squared quantities, the bound is squared too.
*/
Modified: lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYZBounds.java
URL: http://svn.apache.org/viewvc/lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYZBounds.java?rev=1701837&r1=1701836&r2=1701837&view=diff
==============================================================================
--- lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYZBounds.java (original)
+++ lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYZBounds.java Tue Sep 8 16:29:22 2015
@@ -30,7 +30,7 @@ public class XYZBounds implements Bounds
* of the shape, and we cannot guarantee that without making MINIMUM_RESOLUTION
* unacceptably large.
*/
- protected static final double FUDGE_FACTOR = Vector.MINIMUM_RESOLUTION * 2e6;
+ protected static final double FUDGE_FACTOR = Vector.MINIMUM_RESOLUTION * 2.0;
/** Minimum x */
protected Double minX = null;
Modified: lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYZSolid.java
URL: http://svn.apache.org/viewvc/lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYZSolid.java?rev=1701837&r1=1701836&r2=1701837&view=diff
==============================================================================
--- lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYZSolid.java (original)
+++ lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYZSolid.java Tue Sep 8 16:29:22 2015
@@ -186,7 +186,13 @@ public class XYZSolid extends BaseXYZSol
// First construct a perpendicular plane that will allow us to find a sample point.
// This plane is vertical and goes through the points (0,0,0) and (1,0,0)
// Then use it to compute a sample point.
- minXEdges = new GeoPoint[]{minXPlane.getSampleIntersectionPoint(planetModel, xVerticalPlane)};
+ final GeoPoint intPoint = minXPlane.getSampleIntersectionPoint(planetModel, xVerticalPlane);
+ if (intPoint != null) {
+ minXEdges = new GeoPoint[]{intPoint};
+ } else {
+ // No intersection found?
+ minXEdges = EMPTY_POINTS;
+ }
} else {
minXEdges = EMPTY_POINTS;
}
@@ -199,7 +205,12 @@ public class XYZSolid extends BaseXYZSol
// First construct a perpendicular plane that will allow us to find a sample point.
// This plane is vertical and goes through the points (0,0,0) and (1,0,0)
// Then use it to compute a sample point.
- maxXEdges = new GeoPoint[]{maxXPlane.getSampleIntersectionPoint(planetModel, xVerticalPlane)};
+ final GeoPoint intPoint = maxXPlane.getSampleIntersectionPoint(planetModel, xVerticalPlane);
+ if (intPoint != null) {
+ maxXEdges = new GeoPoint[]{intPoint};
+ } else {
+ maxXEdges = EMPTY_POINTS;
+ }
} else {
maxXEdges = EMPTY_POINTS;
}
@@ -212,7 +223,12 @@ public class XYZSolid extends BaseXYZSol
// First construct a perpendicular plane that will allow us to find a sample point.
// This plane is vertical and goes through the points (0,0,0) and (0,1,0)
// Then use it to compute a sample point.
- minYEdges = new GeoPoint[]{minYPlane.getSampleIntersectionPoint(planetModel, yVerticalPlane)};
+ final GeoPoint intPoint = minYPlane.getSampleIntersectionPoint(planetModel, yVerticalPlane);
+ if (intPoint != null) {
+ minYEdges = new GeoPoint[]{intPoint};
+ } else {
+ minYEdges = EMPTY_POINTS;
+ }
} else {
minYEdges = EMPTY_POINTS;
}
@@ -225,7 +241,12 @@ public class XYZSolid extends BaseXYZSol
// First construct a perpendicular plane that will allow us to find a sample point.
// This plane is vertical and goes through the points (0,0,0) and (0,1,0)
// Then use it to compute a sample point.
- maxYEdges = new GeoPoint[]{maxYPlane.getSampleIntersectionPoint(planetModel, yVerticalPlane)};
+ final GeoPoint intPoint = maxYPlane.getSampleIntersectionPoint(planetModel, yVerticalPlane);
+ if (intPoint != null) {
+ maxYEdges = new GeoPoint[]{intPoint};
+ } else {
+ maxYEdges = EMPTY_POINTS;
+ }
} else {
maxYEdges = EMPTY_POINTS;
}
@@ -238,7 +259,12 @@ public class XYZSolid extends BaseXYZSol
// First construct a perpendicular plane that will allow us to find a sample point.
// This plane is vertical and goes through the points (0,0,0) and (1,0,0)
// Then use it to compute a sample point.
- minZEdges = new GeoPoint[]{minZPlane.getSampleIntersectionPoint(planetModel, xVerticalPlane)};
+ final GeoPoint intPoint = minZPlane.getSampleIntersectionPoint(planetModel, xVerticalPlane);
+ if (intPoint != null) {
+ minZEdges = new GeoPoint[]{intPoint};
+ } else {
+ minZEdges = EMPTY_POINTS;
+ }
} else {
minZEdges = EMPTY_POINTS;
}
@@ -251,7 +277,12 @@ public class XYZSolid extends BaseXYZSol
// First construct a perpendicular plane that will allow us to find a sample point.
// This plane is vertical and goes through the points (0,0,0) and (1,0,0) (that is, its orientation doesn't matter)
// Then use it to compute a sample point.
- maxZEdges = new GeoPoint[]{maxZPlane.getSampleIntersectionPoint(planetModel, xVerticalPlane)};
+ final GeoPoint intPoint = maxZPlane.getSampleIntersectionPoint(planetModel, xVerticalPlane);
+ if (intPoint != null) {
+ maxZEdges = new GeoPoint[]{intPoint};
+ } else {
+ maxZEdges = EMPTY_POINTS;
+ }
} else {
maxZEdges = EMPTY_POINTS;
}
@@ -291,6 +322,16 @@ public class XYZSolid extends BaseXYZSol
return OVERLAPS;
}
+ /*
+ for (GeoPoint p : getEdgePoints()) {
+ System.err.println(" Edge point "+p+" path.isWithin()? "+path.isWithin(p));
+ }
+
+ for (GeoPoint p : path.getEdgePoints()) {
+ System.err.println(" path edge point "+p+" isWithin()? "+isWithin(p)+" minx="+minXPlane.evaluate(p)+" maxx="+maxXPlane.evaluate(p)+" miny="+minYPlane.evaluate(p)+" maxy="+maxYPlane.evaluate(p)+" minz="+minZPlane.evaluate(p)+" maxz="+maxZPlane.evaluate(p));
+ }
+ */
+
//System.err.println(this+" getrelationship with "+path);
final int insideRectangle = isShapeInsideArea(path);
if (insideRectangle == SOME_INSIDE) {
Modified: lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYdZSolid.java
URL: http://svn.apache.org/viewvc/lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYdZSolid.java?rev=1701837&r1=1701836&r2=1701837&view=diff
==============================================================================
--- lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYdZSolid.java (original)
+++ lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYdZSolid.java Tue Sep 8 16:29:22 2015
@@ -113,7 +113,12 @@ public class XYdZSolid extends BaseXYZSo
// First construct a perpendicular plane that will allow us to find a sample point.
// This plane is vertical and goes through the points (0,0,0) and (1,0,0)
// Then use it to compute a sample point.
- zEdges = new GeoPoint[]{zPlane.getSampleIntersectionPoint(planetModel, xVerticalPlane)};
+ final GeoPoint intPoint = zPlane.getSampleIntersectionPoint(planetModel, xVerticalPlane);
+ if (intPoint != null) {
+ zEdges = new GeoPoint[]{intPoint};
+ } else {
+ zEdges = EMPTY_POINTS;
+ }
} else {
zEdges= EMPTY_POINTS;
}
Modified: lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XdYZSolid.java
URL: http://svn.apache.org/viewvc/lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XdYZSolid.java?rev=1701837&r1=1701836&r2=1701837&view=diff
==============================================================================
--- lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XdYZSolid.java (original)
+++ lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XdYZSolid.java Tue Sep 8 16:29:22 2015
@@ -113,7 +113,12 @@ public class XdYZSolid extends BaseXYZSo
// First construct a perpendicular plane that will allow us to find a sample point.
// This plane is vertical and goes through the points (0,0,0) and (0,1,0)
// Then use it to compute a sample point.
- yEdges = new GeoPoint[]{yPlane.getSampleIntersectionPoint(planetModel, yVerticalPlane)};
+ final GeoPoint intPoint = yPlane.getSampleIntersectionPoint(planetModel, yVerticalPlane);
+ if (intPoint != null) {
+ yEdges = new GeoPoint[]{intPoint};
+ } else {
+ yEdges = EMPTY_POINTS;
+ }
} else {
yEdges = EMPTY_POINTS;
}
Modified: lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/dXYZSolid.java
URL: http://svn.apache.org/viewvc/lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/dXYZSolid.java?rev=1701837&r1=1701836&r2=1701837&view=diff
==============================================================================
--- lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/dXYZSolid.java (original)
+++ lucene/dev/branches/branch_5x/lucene/spatial3d/src/java/org/apache/lucene/geo3d/dXYZSolid.java Tue Sep 8 16:29:22 2015
@@ -114,7 +114,12 @@ public class dXYZSolid extends BaseXYZSo
// First construct a perpendicular plane that will allow us to find a sample point.
// This plane is vertical and goes through the points (0,0,0) and (1,0,0)
// Then use it to compute a sample point.
- xEdges = new GeoPoint[]{xPlane.getSampleIntersectionPoint(planetModel, xVerticalPlane)};
+ final GeoPoint intPoint = xPlane.getSampleIntersectionPoint(planetModel, xVerticalPlane);
+ if (intPoint != null) {
+ xEdges = new GeoPoint[]{intPoint};
+ } else {
+ xEdges = EMPTY_POINTS;
+ }
} else {
xEdges = EMPTY_POINTS;
}
Modified: lucene/dev/branches/branch_5x/lucene/spatial3d/src/test/org/apache/lucene/bkdtree3d/TestGeo3DPointField.java
URL: http://svn.apache.org/viewvc/lucene/dev/branches/branch_5x/lucene/spatial3d/src/test/org/apache/lucene/bkdtree3d/TestGeo3DPointField.java?rev=1701837&r1=1701836&r2=1701837&view=diff
==============================================================================
--- lucene/dev/branches/branch_5x/lucene/spatial3d/src/test/org/apache/lucene/bkdtree3d/TestGeo3DPointField.java (original)
+++ lucene/dev/branches/branch_5x/lucene/spatial3d/src/test/org/apache/lucene/bkdtree3d/TestGeo3DPointField.java Tue Sep 8 16:29:22 2015
@@ -227,6 +227,21 @@ public class TestGeo3DPointField extends
}
}
+ private static PlanetModel getPlanetModel() {
+ if (random().nextBoolean()) {
+ // Use one of the earth models:
+ if (random().nextBoolean()) {
+ return PlanetModel.WGS84;
+ } else {
+ return PlanetModel.SPHERE;
+ }
+ } else {
+ // Make a randomly squashed planet:
+ double oblateness = random().nextDouble() * 0.5 - 0.25;
+ return new PlanetModel(1.0 + oblateness, 1.0 - oblateness);
+ }
+ }
+
public void testBKDRandom() throws Exception {
final List<Point> points = new ArrayList<>();
int numPoints = atLeast(10000);
@@ -236,12 +251,7 @@ public class TestGeo3DPointField extends
int maxPointsSortInHeap = TestUtil.nextInt(random(), maxPointsInLeaf, 1024*1024);
- PlanetModel planetModel;
- if (random().nextBoolean()) {
- planetModel = PlanetModel.WGS84;
- } else {
- planetModel = PlanetModel.SPHERE;
- }
+ PlanetModel planetModel = getPlanetModel();
final double planetMax = planetModel.getMaximumMagnitude();
BKD3DTreeWriter w = new BKD3DTreeWriter(maxPointsInLeaf, maxPointsSortInHeap);
@@ -439,12 +449,7 @@ public class TestGeo3DPointField extends
/** Tests consistency of GeoArea.getRelationship vs GeoShape.isWithin */
public void testGeo3DRelations() throws Exception {
- PlanetModel planetModel;
- if (random().nextBoolean()) {
- planetModel = PlanetModel.WGS84;
- } else {
- planetModel = PlanetModel.SPHERE;
- }
+ PlanetModel planetModel = getPlanetModel();
int numDocs = atLeast(1000);
if (VERBOSE) {
@@ -699,9 +704,7 @@ public class TestGeo3DPointField extends
}
if (fail) {
- if (VERBOSE) {
- System.out.print(sw.toString());
- }
+ System.out.print(sw.toString());
fail("invalid hits for shape=" + shape);
}
}
@@ -898,12 +901,7 @@ public class TestGeo3DPointField extends
int maxPointsSortInHeap = TestUtil.nextInt(random(), maxPointsInLeaf, 1024*1024);
IndexWriterConfig iwc = newIndexWriterConfig();
- final PlanetModel planetModel;
- if (random().nextBoolean()) {
- planetModel = PlanetModel.WGS84;
- } else {
- planetModel = PlanetModel.SPHERE;
- }
+ final PlanetModel planetModel = getPlanetModel();
// Else we can get O(N^2) merging:
int mbd = iwc.getMaxBufferedDocs();
Modified: lucene/dev/branches/branch_5x/lucene/spatial3d/src/test/org/apache/lucene/geo3d/GeoCircleTest.java
URL: http://svn.apache.org/viewvc/lucene/dev/branches/branch_5x/lucene/spatial3d/src/test/org/apache/lucene/geo3d/GeoCircleTest.java?rev=1701837&r1=1701836&r2=1701837&view=diff
==============================================================================
--- lucene/dev/branches/branch_5x/lucene/spatial3d/src/test/org/apache/lucene/geo3d/GeoCircleTest.java (original)
+++ lucene/dev/branches/branch_5x/lucene/spatial3d/src/test/org/apache/lucene/geo3d/GeoCircleTest.java Tue Sep 8 16:29:22 2015
@@ -102,6 +102,16 @@ public class GeoCircleTest extends Lucen
GeoPoint p2;
int relationship;
+ // ...
+ c = new GeoCircle(PlanetModel.WGS84, -0.005931145568901605, -0.001942031539653079, 1.2991918568260272E-4);
+ area = GeoAreaFactory.makeGeoArea(PlanetModel.WGS84, 1.001098377143621, 1.001100011578687, -0.00207467080358696, -0.0018136665346280983, -0.006067808248760161, -0.005807683665759485);
+ p1 = new GeoPoint(PlanetModel.WGS84, -0.00591253844632244, -0.0020069187259065093);
+ p2 = new GeoPoint(1.001099185736782, -0.0020091272069679327, -0.005919118245803968);
+ assertTrue(c.isWithin(p1));
+ assertTrue(area.isWithin(p1));
+ relationship = area.getRelationship(c);
+ assertTrue(relationship != GeoArea.DISJOINT);
+
// Twelfth BKD discovered failure
c = new GeoCircle(PlanetModel.WGS84,-0.00824379317765984,-0.0011677469001838581,0.0011530035396910402);
p1 = new GeoPoint(PlanetModel.WGS84,-0.006505092992723671,0.007654282718327381);
Modified: lucene/dev/branches/branch_5x/lucene/spatial3d/src/test/org/apache/lucene/geo3d/GeoPathTest.java
URL: http://svn.apache.org/viewvc/lucene/dev/branches/branch_5x/lucene/spatial3d/src/test/org/apache/lucene/geo3d/GeoPathTest.java?rev=1701837&r1=1701836&r2=1701837&view=diff
==============================================================================
--- lucene/dev/branches/branch_5x/lucene/spatial3d/src/test/org/apache/lucene/geo3d/GeoPathTest.java (original)
+++ lucene/dev/branches/branch_5x/lucene/spatial3d/src/test/org/apache/lucene/geo3d/GeoPathTest.java Tue Sep 8 16:29:22 2015
@@ -130,7 +130,23 @@ public class GeoPathTest {
public void testGetRelationship() {
GeoArea rect;
GeoPath p;
-
+ GeoPath c;
+ GeoPoint point;
+ GeoPoint pointApprox;
+ int relationship;
+ GeoArea area;
+ PlanetModel planetModel;
+
+ planetModel = new PlanetModel(1.151145876105594, 0.8488541238944061);
+ c = new GeoPath(planetModel, 0.008726646259971648);
+ c.addPoint(-0.6925658899376476, 0.6316613927914589);
+ c.addPoint(0.27828548161836364, 0.6785795524104564);
+ c.done();
+ point = new GeoPoint(planetModel,-0.49298555067758226, 0.9892440995026406);
+ pointApprox = new GeoPoint(0.5110940362119821, 0.7774603209946239, -0.49984312299556544);
+ area = GeoAreaFactory.makeGeoArea(planetModel, 0.49937141144985997, 0.5161765426256085, 0.3337218719537796,0.8544419570901649, -0.6347692823688085, 0.3069696588119369);
+ assertTrue(!c.isWithin(point));
+
// Start by testing the basic kinds of relationship, increasing in order of difficulty.
p = new GeoPath(PlanetModel.SPHERE, 0.1);
@@ -170,7 +186,26 @@ public class GeoPathTest {
GeoPoint point;
int relationship;
GeoArea area;
-
+ PlanetModel planetModel;
+
+ planetModel = new PlanetModel(0.751521665790406,1.248478334209594);
+ c = new GeoPath(planetModel, 0.7504915783575618);
+ c.addPoint(0.10869761172400265, 0.08895880215465272);
+ c.addPoint(0.22467878641991612, 0.10972973084229565);
+ c.addPoint(-0.7398772468744732, -0.4465812941383364);
+ c.addPoint(-0.18462055300079366, -0.6713857796763727);
+ c.done();
+ point = new GeoPoint(planetModel,-0.626645355125733,-1.409304625439381);
+ xyzb = new XYZBounds();
+ c.getBounds(xyzb);
+ area = GeoAreaFactory.makeGeoArea(planetModel,
+ xyzb.getMinimumX(), xyzb.getMaximumX(), xyzb.getMinimumY(), xyzb.getMaximumY(), xyzb.getMinimumZ(), xyzb.getMaximumZ());
+ relationship = area.getRelationship(c);
+ assertTrue(relationship == GeoArea.WITHIN || relationship == GeoArea.OVERLAPS);
+ assertTrue(area.isWithin(point));
+ // No longer true due to fixed GeoPath waypoints.
+ //assertTrue(c.isWithin(point));
+
c = new GeoPath(PlanetModel.WGS84, 0.6894050545377601);
c.addPoint(-0.0788176065762948, 0.9431251741731624);
c.addPoint(0.510387871458147, 0.5327078872484678);