lucene-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From mikemcc...@apache.org
Subject svn commit: r1701835 - in /lucene/dev/trunk/lucene: ./ spatial3d/src/java/org/apache/lucene/geo3d/ spatial3d/src/test/org/apache/lucene/bkdtree3d/ spatial3d/src/test/org/apache/lucene/geo3d/
Date Tue, 08 Sep 2015 16:25:30 GMT
Author: mikemccand
Date: Tue Sep  8 16:25:30 2015
New Revision: 1701835

URL: http://svn.apache.org/r1701835
Log:
LUCENE-6776: Fix geo3d math to handle randomly squashed planets

Modified:
    lucene/dev/trunk/lucene/CHANGES.txt
    lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/GeoPath.java
    lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/Plane.java
    lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/PlanetModel.java
    lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/Vector.java
    lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYZBounds.java
    lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYZSolid.java
    lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYdZSolid.java
    lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XdYZSolid.java
    lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/dXYZSolid.java
    lucene/dev/trunk/lucene/spatial3d/src/test/org/apache/lucene/bkdtree3d/TestGeo3DPointField.java
    lucene/dev/trunk/lucene/spatial3d/src/test/org/apache/lucene/geo3d/GeoCircleTest.java
    lucene/dev/trunk/lucene/spatial3d/src/test/org/apache/lucene/geo3d/GeoPathTest.java

Modified: lucene/dev/trunk/lucene/CHANGES.txt
URL: http://svn.apache.org/viewvc/lucene/dev/trunk/lucene/CHANGES.txt?rev=1701835&r1=1701834&r2=1701835&view=diff
==============================================================================
--- lucene/dev/trunk/lucene/CHANGES.txt (original)
+++ lucene/dev/trunk/lucene/CHANGES.txt Tue Sep  8 16:25:30 2015
@@ -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.

Modified: lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/GeoPath.java
URL: http://svn.apache.org/viewvc/lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/GeoPath.java?rev=1701835&r1=1701834&r2=1701835&view=diff
==============================================================================
--- lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/GeoPath.java (original)
+++ lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/GeoPath.java Tue Sep
 8 16:25:30 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/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/Plane.java
URL: http://svn.apache.org/viewvc/lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/Plane.java?rev=1701835&r1=1701834&r2=1701835&view=diff
==============================================================================
--- lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/Plane.java (original)
+++ lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/Plane.java Tue Sep
 8 16:25:30 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/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/PlanetModel.java
URL: http://svn.apache.org/viewvc/lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/PlanetModel.java?rev=1701835&r1=1701834&r2=1701835&view=diff
==============================================================================
--- lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/PlanetModel.java (original)
+++ lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/PlanetModel.java Tue
Sep  8 16:25:30 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/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/Vector.java
URL: http://svn.apache.org/viewvc/lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/Vector.java?rev=1701835&r1=1701834&r2=1701835&view=diff
==============================================================================
--- lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/Vector.java (original)
+++ lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/Vector.java Tue Sep
 8 16:25:30 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/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYZBounds.java
URL: http://svn.apache.org/viewvc/lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYZBounds.java?rev=1701835&r1=1701834&r2=1701835&view=diff
==============================================================================
--- lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYZBounds.java (original)
+++ lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYZBounds.java Tue
Sep  8 16:25:30 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/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYZSolid.java
URL: http://svn.apache.org/viewvc/lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYZSolid.java?rev=1701835&r1=1701834&r2=1701835&view=diff
==============================================================================
--- lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYZSolid.java (original)
+++ lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYZSolid.java Tue Sep
 8 16:25:30 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/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYdZSolid.java
URL: http://svn.apache.org/viewvc/lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYdZSolid.java?rev=1701835&r1=1701834&r2=1701835&view=diff
==============================================================================
--- lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYdZSolid.java (original)
+++ lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XYdZSolid.java Tue
Sep  8 16:25:30 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/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XdYZSolid.java
URL: http://svn.apache.org/viewvc/lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XdYZSolid.java?rev=1701835&r1=1701834&r2=1701835&view=diff
==============================================================================
--- lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XdYZSolid.java (original)
+++ lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/XdYZSolid.java Tue
Sep  8 16:25:30 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/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/dXYZSolid.java
URL: http://svn.apache.org/viewvc/lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/dXYZSolid.java?rev=1701835&r1=1701834&r2=1701835&view=diff
==============================================================================
--- lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/dXYZSolid.java (original)
+++ lucene/dev/trunk/lucene/spatial3d/src/java/org/apache/lucene/geo3d/dXYZSolid.java Tue
Sep  8 16:25:30 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/trunk/lucene/spatial3d/src/test/org/apache/lucene/bkdtree3d/TestGeo3DPointField.java
URL: http://svn.apache.org/viewvc/lucene/dev/trunk/lucene/spatial3d/src/test/org/apache/lucene/bkdtree3d/TestGeo3DPointField.java?rev=1701835&r1=1701834&r2=1701835&view=diff
==============================================================================
--- lucene/dev/trunk/lucene/spatial3d/src/test/org/apache/lucene/bkdtree3d/TestGeo3DPointField.java
(original)
+++ lucene/dev/trunk/lucene/spatial3d/src/test/org/apache/lucene/bkdtree3d/TestGeo3DPointField.java
Tue Sep  8 16:25:30 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 {
     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();
 
-    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();

Modified: lucene/dev/trunk/lucene/spatial3d/src/test/org/apache/lucene/geo3d/GeoCircleTest.java
URL: http://svn.apache.org/viewvc/lucene/dev/trunk/lucene/spatial3d/src/test/org/apache/lucene/geo3d/GeoCircleTest.java?rev=1701835&r1=1701834&r2=1701835&view=diff
==============================================================================
--- lucene/dev/trunk/lucene/spatial3d/src/test/org/apache/lucene/geo3d/GeoCircleTest.java
(original)
+++ lucene/dev/trunk/lucene/spatial3d/src/test/org/apache/lucene/geo3d/GeoCircleTest.java
Tue Sep  8 16:25:30 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/trunk/lucene/spatial3d/src/test/org/apache/lucene/geo3d/GeoPathTest.java
URL: http://svn.apache.org/viewvc/lucene/dev/trunk/lucene/spatial3d/src/test/org/apache/lucene/geo3d/GeoPathTest.java?rev=1701835&r1=1701834&r2=1701835&view=diff
==============================================================================
--- lucene/dev/trunk/lucene/spatial3d/src/test/org/apache/lucene/geo3d/GeoPathTest.java (original)
+++ lucene/dev/trunk/lucene/spatial3d/src/test/org/apache/lucene/geo3d/GeoPathTest.java Tue
Sep  8 16:25:30 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);



Mime
View raw message