LUCENE-7221: Fix broken bounds for GeoCircles.

This commit is contained in:
Karl Wright 2016-04-15 03:01:11 -04:00
parent 7c098913e2
commit 5f0245febe
3 changed files with 80 additions and 25 deletions

View File

@ -926,17 +926,23 @@ public class Plane extends Vector {
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);
// Valid?
if (Math.abs(m) >= MINIMUM_RESOLUTION) {
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 {
// This is a plane of the form A=n B=0 C=0. We can set a bound only by noting the D value.
boundsInfo.addXValue(-D/A);
}
} else if (sqrtTerm > 0.0) {
// Two solutions
final double sqrtResult = Math.sqrt(sqrtTerm);
@ -1089,17 +1095,23 @@ public class Plane extends Vector {
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);
// Valid?
if (Math.abs(m) >= MINIMUM_RESOLUTION) {
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 {
// This is a plane of the form A=0 B=n C=0. We can set a bound only by noting the D value.
boundsInfo.addYValue(-D/B);
}
} else if (sqrtTerm > 0.0) {
// Two solutions
final double sqrtResult = Math.sqrt(sqrtTerm);

View File

@ -195,7 +195,14 @@ public class XYZBounds implements Bounds {
@Override
public Bounds addXValue(final GeoPoint point) {
final double x = point.x;
return addXValue(point.x);
}
/** Add a specific X value.
* @param x is the value to add.
* @return the bounds object.
*/
public Bounds addXValue(final double x) {
final double small = x - FUDGE_FACTOR;
if (minX == null || minX > small) {
minX = new Double(small);
@ -209,7 +216,14 @@ public class XYZBounds implements Bounds {
@Override
public Bounds addYValue(final GeoPoint point) {
final double y = point.y;
return addYValue(point.y);
}
/** Add a specific Y value.
* @param y is the value to add.
* @return the bounds object.
*/
public Bounds addYValue(final double y) {
final double small = y - FUDGE_FACTOR;
if (minY == null || minY > small) {
minY = new Double(small);
@ -223,7 +237,14 @@ public class XYZBounds implements Bounds {
@Override
public Bounds addZValue(final GeoPoint point) {
final double z = point.z;
return addZValue(point.z);
}
/** Add a specific Z value.
* @param z is the value to add.
* @return the bounds object.
*/
public Bounds addZValue(final double z) {
final double small = z - FUDGE_FACTOR;
if (minZ == null || minZ > small) {
minZ = new Double(small);
@ -264,4 +285,9 @@ public class XYZBounds implements Bounds {
return this;
}
@Override
public String toString() {
return "XYZBounds: [xmin="+minX+" xmax="+maxX+" ymin="+minY+" ymax="+maxY+" zmin="+minZ+" zmax="+maxZ+"]";
}
}

View File

@ -392,4 +392,21 @@ public class GeoCircleTest extends LuceneTestCase {
}
@Test
public void testBoundsFailureCase1() {
// lat=2.7399499693409367E-13, lon=-3.141592653589793([X=-1.0011188539924791, Y=-1.226017000107956E-16, Z=2.743015573303327E-13])], radius=2.1814042682464985
final GeoCircle gc = GeoCircleFactory.makeGeoCircle(PlanetModel.WGS84, 2.7399499693409367E-13, -3.141592653589793, 2.1814042682464985);
// With a circle like this, zmin should equal zmax, and xmin should be PlanetModel.minimumX.
final GeoPoint gp = new GeoPoint(0.0054866241253590815, -0.004009749293376541, 0.997739304376186);
final GeoPoint gpOnSurface = PlanetModel.WGS84.createSurfacePoint(gp);
final XYZBounds bounds = new XYZBounds();
gc.getBounds(bounds);
//System.out.println("Bounds: "+bounds);
final XYZSolid solid = XYZSolidFactory.makeXYZSolid(PlanetModel.WGS84, bounds.getMinimumX(), bounds.getMaximumX(), bounds.getMinimumY(), bounds.getMaximumY(), bounds.getMinimumZ(), bounds.getMaximumZ());
assertTrue(gc.isWithin(gpOnSurface));
assertTrue(gc.isWithin(gp));
assertTrue(solid.isWithin(gpOnSurface)); // This fails
assertTrue(solid.isWithin(gp));
}
}