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

LUCENE-6776: Fix geo3d math to handle randomly squashed planets 

git-svn-id: https://svn.apache.org/repos/asf/lucene/dev/trunk@1701835 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Michael McCandless 2015-09-08 16:25:30 +00:00
parent f92df1c93b
commit 9dcef8c8c3
13 changed files with 475 additions and 111 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
lucene/CHANGES.txt
View File

@ -118,6 +118,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.

18
lucene/spatial3d/src/java/org/apache/lucene/geo3d/GeoPath.java
View File

@ -475,7 +475,14 @@ public class GeoPath extends GeoBaseDistanceShape {
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 GeoBaseDistanceShape {
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.

396
lucene/spatial3d/src/java/org/apache/lucene/geo3d/Plane.java
View File

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

2
lucene/spatial3d/src/java/org/apache/lucene/geo3d/PlanetModel.java
View File

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

2
lucene/spatial3d/src/java/org/apache/lucene/geo3d/Vector.java
View File

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

2
lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYZBounds.java
View File

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

53
lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYZSolid.java
View File

@ -186,7 +186,13 @@ public class XYZSolid extends BaseXYZSolid {
// 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 BaseXYZSolid {
// 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 BaseXYZSolid {
// 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 BaseXYZSolid {
// 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 BaseXYZSolid {
// 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 BaseXYZSolid {
// 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 BaseXYZSolid {
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) {

7
lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYdZSolid.java
View File

@ -113,7 +113,12 @@ public class XYdZSolid extends BaseXYZSolid {
// 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;
}

7
lucene/spatial3d/src/java/org/apache/lucene/geo3d/XdYZSolid.java
View File

@ -113,7 +113,12 @@ public class XdYZSolid extends BaseXYZSolid {
// 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;
}

7
lucene/spatial3d/src/java/org/apache/lucene/geo3d/dXYZSolid.java
View File

@ -114,7 +114,12 @@ public class dXYZSolid extends BaseXYZSolid {
// 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;
}

40
lucene/spatial3d/src/test/org/apache/lucene/bkdtree3d/TestGeo3DPointField.java
View File

@ -227,6 +227,21 @@ public class TestGeo3DPointField extends LuceneTestCase {
}
}
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 {
List<Point> points = new ArrayList<>();
int numPoints = atLeast(10000);
@ -236,12 +251,7 @@ public class TestGeo3DPointField extends LuceneTestCase {
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 LuceneTestCase {
/** 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 LuceneTestCase {
}
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 LuceneTestCase {
int maxPointsSortInHeap = TestUtil.nextInt(random(), maxPointsInLeaf, 1024*1024);
IndexWriterConfig iwc = newIndexWriterConfig();
PlanetModel planetModel;
if (random().nextBoolean()) {
planetModel = PlanetModel.WGS84;
} else {
planetModel = PlanetModel.SPHERE;
}
PlanetModel planetModel = getPlanetModel();
// Else we can get O(N^2) merging:
int mbd = iwc.getMaxBufferedDocs();

10
lucene/spatial3d/src/test/org/apache/lucene/geo3d/GeoCircleTest.java
View File

@ -102,6 +102,16 @@ public class GeoCircleTest extends LuceneTestCase {
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);

39
lucene/spatial3d/src/test/org/apache/lucene/geo3d/GeoPathTest.java
View File

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