LUCENE-6196: committing Karl's latest patch

https://reviews.apache.org/r/33476/ (diff #2)
Addresses getBoundingBox issues

git-svn-id: https://svn.apache.org/repos/asf/lucene/dev/branches/lucene6196@1675747 13f79535-47bb-0310-9956-ffa450edef68

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/Bounds.java

@@ -259,6 +259,7 @@ public class Bounds
         if (!noLongitudeBound) {
             // Get a longitude value
             double longitude = Math.atan2(y,x);
+            //System.err.println(" add longitude bound at "+longitude * 180.0/Math.PI);
             addLongitudeBound(longitude);
         }
         if (!noTopLatitudeBound || !noBottomLatitudeBound) {

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/GeoCircle.java

@@ -74,6 +74,14 @@ public class GeoCircle extends GeoBaseExtendedShape implements GeoDistanceShape,
         return cutoffAngle;
     }
 
+    /** Returns the center of a circle into which the area will be inscribed.
+    *@return the center.
+    */
+    @Override
+    public GeoPoint getCenter() {
+        return center;
+    }
+
     /** Compute an estimate of "distance" to the GeoPoint.
     * A return value of Double.MAX_VALUE should be returned for
     * points outside of the shape.
@@ -221,6 +229,7 @@ public class GeoCircle extends GeoBaseExtendedShape implements GeoDistanceShape,
     public Bounds getBounds(Bounds bounds)
     {
         bounds = super.getBounds(bounds);
+        bounds.addPoint(center);
         circlePlane.recordBounds(bounds);
         return bounds;
     }

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/GeoCompositeMembershipShape.java

@@ -39,9 +39,12 @@ public class GeoCompositeMembershipShape implements GeoMembershipShape
     @Override
     public boolean isWithin(final Vector point)
     {
+        //System.err.println("Checking whether point "+point+" is within Composite");
         for (GeoMembershipShape shape : shapes) {
-            if (shape.isWithin(point))
+            if (shape.isWithin(point)) {
+                //System.err.println(" Point is within "+shape);
                 return true;
+            }
         }
         return false;
     }

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/GeoConvexPolygon.java

@@ -37,6 +37,8 @@ public class GeoConvexPolygon extends GeoBaseExtendedShape implements GeoMembers
     
     protected GeoPoint[] edgePoints = null;
     
+    protected double fullDistance = 0.0;
+    
     /** Create a convex polygon from a list of points.  The first point must be on the
     * external edge.
     */
@@ -98,14 +100,17 @@ public class GeoConvexPolygon extends GeoBaseExtendedShape implements GeoMembers
         if (points.size() < 3)
             throw new IllegalArgumentException("Polygon needs at least three points.");
         // Time to construct the planes.  If the polygon is truly convex, then any adjacent point
+        // to a segment can provide an interior measurement.
         edges = new SidedPlane[points.size()];
         notableEdgePoints = new GeoPoint[points.size()][];
         internalEdges = new boolean[points.size()];
-        // to a segment can provide an interior measurement.
         for (int i = 0; i < points.size(); i++) {
             final GeoPoint start = points.get(i);
             final boolean isInternalEdge = (isInternalEdges!=null?(i == isInternalEdges.size()?isInternalReturnEdge:isInternalEdges.get(i)):false);
             final GeoPoint end = points.get(legalIndex(i+1));
+            final double distance = start.arcDistance(end);
+            if (distance > fullDistance)
+                fullDistance = distance;
             final GeoPoint check = points.get(legalIndex(i+2));
             final SidedPlane sp = new SidedPlane(check,start,end);
             //System.out.println("Created edge "+sp+" using start="+start+" end="+end+" check="+check);
@@ -204,7 +209,7 @@ public class GeoConvexPolygon extends GeoBaseExtendedShape implements GeoMembers
     public Bounds getBounds(Bounds bounds)
     {
         bounds = super.getBounds(bounds);
-        
+
         // Add all the points
         for (final GeoPoint point : points) {
             bounds.addPoint(point);
@@ -224,6 +229,10 @@ public class GeoConvexPolygon extends GeoBaseExtendedShape implements GeoMembers
             edge.recordBounds(bounds,membershipBounds);
         }
 
+        if (fullDistance >= Math.PI * 0.5) {
+            // We can't reliably assume that bounds did its longitude calculation right, so we force it to be unbounded.
+            bounds.noLongitudeBound();
+        }
         return bounds;
     }
 
@@ -250,7 +259,12 @@ public class GeoConvexPolygon extends GeoBaseExtendedShape implements GeoMembers
 
     @Override
     public String toString() {
-        return "GeoConvexPolygon: {" + points + "}";
+        StringBuilder edgeString = new StringBuilder("{");
+        for (int i = 0; i < edges.length; i++) {
+            edgeString.append(edges[i]).append(" internal? ").append(internalEdges[i]).append("; ");
+        }
+        edgeString.append("}");
+        return "GeoConvexPolygon: {points=" + points + " edges="+edgeString+"}";
     }
 }
   

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/GeoDegenerateHorizontalLine.java

@@ -131,7 +131,15 @@ public class GeoDegenerateHorizontalLine extends GeoBBoxBase
         double bottomAngle = centerPoint.arcDistance(LHC);
         return Math.max(topAngle,bottomAngle);
     }
-      
+
+    /** Returns the center of a circle into which the area will be inscribed.
+    *@return the center.
+    */
+    @Override
+    public GeoPoint getCenter() {
+        return centerPoint;
+    }
+
     @Override
     public GeoPoint[] getEdgePoints()
     {

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/GeoDegenerateLatitudeZone.java

@@ -68,6 +68,15 @@ public class GeoDegenerateLatitudeZone extends GeoBBoxBase
         return Math.PI;
     }
 
+    /** Returns the center of a circle into which the area will be inscribed.
+    *@return the center.
+    */
+    @Override
+    public GeoPoint getCenter() {
+        // Totally arbitrary
+        return interiorPoint;
+    }
+
     @Override
     public GeoPoint[] getEdgePoints()
     {

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/GeoDegenerateLongitudeSlice.java

@@ -83,7 +83,15 @@ public class GeoDegenerateLongitudeSlice extends GeoBBoxBase
     {
         return Math.PI * 0.5;
     }
-      
+
+    /** Returns the center of a circle into which the area will be inscribed.
+    *@return the center.
+    */
+    @Override
+    public GeoPoint getCenter() {
+        return interiorPoint;
+    }
+
     @Override
     public GeoPoint[] getEdgePoints()
     {

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/GeoDegeneratePoint.java

@@ -144,6 +144,14 @@ public class GeoDegeneratePoint extends GeoPoint implements GeoBBox
         return 0.0;
     }
 
+    /** Returns the center of a circle into which the area will be inscribed.
+    *@return the center.
+    */
+    @Override
+    public GeoPoint getCenter() {
+        return this;
+    }
+
     /** Find the spatial relationship between a shape and the current geo area.
      * Note: return value is how the GeoShape relates to the GeoArea, not the
      * other way around. For example, if this GeoArea is entirely within the

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/GeoDegenerateVerticalLine.java

@@ -127,7 +127,15 @@ public class GeoDegenerateVerticalLine extends GeoBBoxBase
         final double bottomAngle = centerPoint.arcDistance(LHC);
         return Math.max(topAngle,bottomAngle);
     }
-      
+
+    /** Returns the center of a circle into which the area will be inscribed.
+    *@return the center.
+    */
+    @Override
+    public GeoPoint getCenter() {
+        return centerPoint;
+    }
+
     @Override
     public GeoPoint[] getEdgePoints()
     {

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/GeoLatitudeZone.java

@@ -100,6 +100,15 @@ public class GeoLatitudeZone extends GeoBBoxBase
         return maxCosLat * Math.PI;
     }
 
+    /** Returns the center of a circle into which the area will be inscribed.
+    *@return the center.
+    */
+    @Override
+    public GeoPoint getCenter() {
+        // This is totally arbitrary and only a cartesian could agree with it.
+        return interiorPoint;
+    }
+
     @Override
     public GeoPoint[] getEdgePoints()
     {

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/GeoLongitudeSlice.java

@@ -109,7 +109,15 @@ public class GeoLongitudeSlice extends GeoBBoxBase
             extent += Math.PI * 2.0;
         return Math.max(Math.PI * 0.5, extent * 0.5);
     }
-      
+
+    /** Returns the center of a circle into which the area will be inscribed.
+    *@return the center.
+    */
+    @Override
+    public GeoPoint getCenter() {
+        return centerPoint;
+    }
+
     @Override
     public GeoPoint[] getEdgePoints()
     {

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/GeoNorthLatitudeZone.java

@@ -86,6 +86,14 @@ public class GeoNorthLatitudeZone extends GeoBBoxBase
         return maxCosLat * Math.PI;
     }
 
+    /** Returns the center of a circle into which the area will be inscribed.
+    *@return the center.
+    */
+    @Override
+    public GeoPoint getCenter() {
+        return interiorPoint;
+    }
+
     @Override
     public GeoPoint[] getEdgePoints()
     {

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/GeoNorthRectangle.java

@@ -152,7 +152,15 @@ public class GeoNorthRectangle extends GeoBBoxBase
     {
         return edgePoints;
     }
-      
+
+    /** Returns the center of a circle into which the area will be inscribed.
+    *@return the center.
+    */
+    @Override
+    public GeoPoint getCenter() {
+        return centerPoint;
+    }
+    
     @Override
     public boolean intersects(final Plane p, final GeoPoint[] notablePoints, final Membership... bounds)
     {

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/GeoPath.java

@@ -335,12 +335,15 @@ public class GeoPath extends GeoBaseExtendedShape implements GeoDistanceShape
     public Bounds getBounds(Bounds bounds)
     {
         bounds = super.getBounds(bounds);
-        for (SegmentEndpoint pathPoint : points) {
+        // For building bounds, order matters.  We want to traverse
-            pathPoint.getBounds(bounds);
+        // never more than 180 degrees longitude at a pop or we risk having the
-        }
+        // bounds object get itself inverted.  So do the edges first.
         for (PathSegment pathSegment : segments) {
             pathSegment.getBounds(bounds);
         }
+        for (SegmentEndpoint pathPoint : points) {
+            pathPoint.getBounds(bounds);
+        }
         return bounds;
     }
 
@@ -465,6 +468,7 @@ public class GeoPath extends GeoBaseExtendedShape implements GeoDistanceShape
 
         public void getBounds(Bounds bounds)
         {
+            bounds.addPoint(point);
             circlePlane.recordBounds(bounds);
         }
 
@@ -664,6 +668,7 @@ public class GeoPath extends GeoBaseExtendedShape implements GeoDistanceShape
         public void getBounds(Bounds bounds)
         {
             // We need to do all bounding planes as well as corner points
+            bounds.addPoint(start).addPoint(end);
             upperConnectingPlane.recordBounds(startCutoffPlane, bounds, lowerConnectingPlane, endCutoffPlane);
             startCutoffPlane.recordBounds(lowerConnectingPlane, bounds, endCutoffPlane, upperConnectingPlane);
             lowerConnectingPlane.recordBounds(endCutoffPlane, bounds, upperConnectingPlane, startCutoffPlane);
@@ -672,6 +677,12 @@ public class GeoPath extends GeoBaseExtendedShape implements GeoDistanceShape
             lowerConnectingPlane.recordBounds(bounds, upperConnectingPlane, startCutoffPlane, endCutoffPlane);
             startCutoffPlane.recordBounds(bounds, endCutoffPlane, upperConnectingPlane, lowerConnectingPlane);
             endCutoffPlane.recordBounds(bounds, startCutoffPlane, upperConnectingPlane, lowerConnectingPlane);
+            if (fullDistance >= Math.PI * 0.5) {
+                // Too large a segment basically means that we can confuse the Bounds object.  Specifically, if our span exceeds 180 degrees
+                // in longitude (which even a segment whose actual length is less than that might if it goes close to a pole).
+                // Unfortunately, we can get arbitrarily close to the pole, so this may still not work in all cases.
+                bounds.noLongitudeBound();
+            }
         }
         
     }

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/GeoRectangle.java

@@ -165,11 +165,18 @@ public class GeoRectangle extends GeoBBoxBase
     }
       
     @Override
-    public GeoPoint[] getEdgePoints()
+    public GeoPoint[] getEdgePoints() {
-    {
         return edgePoints;
     }
-      
+
+    /** Returns the center of a circle into which the area will be inscribed.
+    *@return the center.
+    */
+    @Override
+    public GeoPoint getCenter() {
+        return centerPoint;
+    }
+
     @Override
     public boolean intersects(final Plane p, final GeoPoint[] notablePoints, final Membership... bounds)
     {

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/GeoSizeable.java

@@ -27,4 +27,9 @@ public interface GeoSizeable
      */
     public double getRadius();
 
+    /** Returns the center of a circle into which the area will be inscribed.
+    *@return the center.
+    */
+    public GeoPoint getCenter();
+    
 }

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/GeoSouthLatitudeZone.java

@@ -84,6 +84,14 @@ public class GeoSouthLatitudeZone extends GeoBBoxBase
         return maxCosLat * Math.PI;
     }
 
+    /** Returns the center of a circle into which the area will be inscribed.
+    *@return the center.
+    */
+    @Override
+    public GeoPoint getCenter() {
+        return interiorPoint;
+    }
+
     @Override
     public GeoPoint[] getEdgePoints()
     {

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/GeoSouthRectangle.java

@@ -149,7 +149,15 @@ public class GeoSouthRectangle extends GeoBBoxBase
     {
         return edgePoints;
     }
-      
+
+    /** Returns the center of a circle into which the area will be inscribed.
+    *@return the center.
+    */
+    @Override
+    public GeoPoint getCenter() {
+        return centerPoint;
+    }
+
     @Override
     public boolean intersects(final Plane p, final GeoPoint[] notablePoints, final Membership... bounds)
     {

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/GeoWideDegenerateHorizontalLine.java

@@ -142,7 +142,15 @@ public class GeoWideDegenerateHorizontalLine extends GeoBBoxBase
         final double bottomAngle = centerPoint.arcDistance(LHC);
         return Math.max(topAngle,bottomAngle);
     }
-      
+
+    /** Returns the center of a circle into which the area will be inscribed.
+    *@return the center.
+    */
+    @Override
+    public GeoPoint getCenter() {
+        return centerPoint;
+    }
+
     @Override
     public GeoPoint[] getEdgePoints()
     {

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/GeoWideLongitudeSlice.java

@@ -109,7 +109,15 @@ public class GeoWideLongitudeSlice extends GeoBBoxBase
             extent += Math.PI * 2.0;
         return Math.max(Math.PI * 0.5, extent * 0.5);
     }
-      
+
+    /** Returns the center of a circle into which the area will be inscribed.
+    *@return the center.
+    */
+    @Override
+    public GeoPoint getCenter() {
+        return centerPoint;
+    }
+
     @Override
     public GeoPoint[] getEdgePoints()
     {

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/GeoWideNorthRectangle.java

@@ -150,6 +150,14 @@ public class GeoWideNorthRectangle extends GeoBBoxBase
         return Math.max(centerAngle,bottomAngle);
     }
       
+    /** Returns the center of a circle into which the area will be inscribed.
+    *@return the center.
+    */
+    @Override
+    public GeoPoint getCenter() {
+        return centerPoint;
+    }
+
     @Override
     public GeoPoint[] getEdgePoints()
     {

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/GeoWideRectangle.java

@@ -174,7 +174,15 @@ public class GeoWideRectangle extends GeoBBoxBase
     {
         return edgePoints;
     }
-      
+
+    /** Returns the center of a circle into which the area will be inscribed.
+    *@return the center.
+    */
+    @Override
+    public GeoPoint getCenter() {
+        return centerPoint;
+    }
+
     @Override
     public boolean intersects(final Plane p, final GeoPoint[] notablePoints, final Membership... bounds)
     {

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/GeoWideSouthRectangle.java

@@ -147,7 +147,15 @@ public class GeoWideSouthRectangle extends GeoBBoxBase
         final double topAngle = centerPoint.arcDistance(URHC);
         return Math.max(centerAngle,topAngle);
     }
-      
+
+    /** Returns the center of a circle into which the area will be inscribed.
+    *@return the center.
+    */
+    @Override
+    public GeoPoint getCenter() {
+        return centerPoint;
+    }
+
     @Override
     public GeoPoint[] getEdgePoints()
     {

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/GeoWorld.java

@@ -40,6 +40,15 @@ public class GeoWorld extends GeoBBoxBase
         return Math.PI;
     }
 
+    /** Returns the center of a circle into which the area will be inscribed.
+    *@return the center.
+    */
+    @Override
+    public GeoPoint getCenter() {
+        // Totally arbitrary
+        return originPoint;
+    }
+
     @Override
     public boolean isWithin(final Vector point)
     {

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/Plane.java

@@ -68,10 +68,41 @@ public class Plane extends Vector
      *@param v is the vector.
      *@return the result of the evaluation.
      */
+    @Override
     public double evaluate(final Vector v) {
         return super.evaluate(v) + D;
     }
-  
+
+    /** Evaluate the plane equation for a given point, as represented
+     * by a vector.
+     *@param x,y,z is the vector.
+     *@return the result of the evaluation.
+     */
+    @Override
+    public double evaluate(final double x, final double y, final double z) {
+        return super.evaluate(x,y,z) + D;
+    }
+
+    /** Evaluate the plane equation for a given point, as represented
+     * by a vector.
+     *@param v is the vector.
+     *@return true if the result is on the plane.
+     */
+    @Override
+    public boolean evaluateIsZero(final Vector v) {
+        return Math.abs(evaluate(v)) < MINIMUM_RESOLUTION;
+    }
+
+    /** Evaluate the plane equation for a given point, as represented
+     * by a vector.
+     *@param x,y,z is the vector.
+     *@return true if the result is on the plane.
+     */
+    @Override
+    public boolean evaluateIsZero(final double x, final double y, final double z) {
+        return Math.abs(evaluate(x,y,z)) < MINIMUM_RESOLUTION;
+    }
+    
     /** Build a normalized plane, so that the vector is normalized.
      *@return the normalized plane object, or null if the plane is indeterminate.
      */
@@ -90,8 +121,9 @@ public class Plane extends Vector
     */
     protected GeoPoint[] findIntersections(final Plane q, final Membership[] bounds, final Membership[] moreBounds) {
         final Vector lineVector = new Vector(this,q);
-        if (lineVector.x == 0.0 && lineVector.y == 0.0 && lineVector.z == 0.0) {
+        if (Math.abs(lineVector.x) < MINIMUM_RESOLUTION && Math.abs(lineVector.y) < MINIMUM_RESOLUTION && Math.abs(lineVector.z) < MINIMUM_RESOLUTION) {
             // Degenerate case: parallel planes
+            //System.err.println(" planes are parallel - no intersection");
             return NO_POINTS;
         }
 
@@ -123,24 +155,30 @@ public class Plane extends Vector
         final double denomXY = this.x*q.y - this.y*q.x;
         if (Math.abs(denomYZ) >= Math.abs(denomXZ) && Math.abs(denomYZ) >= Math.abs(denomXY)) {
             // X is the biggest, so our point will have x0 = 0.0
-            if (Math.abs(denomYZ) < 1.0e-35)
+            if (Math.abs(denomYZ) < MINIMUM_RESOLUTION_SQUARED) {
+                //System.err.println(" Denominator is zero: no intersection");
                 return NO_POINTS;
+            }
             final double denom = 1.0 / denomYZ;
             x0 = 0.0;
             y0 = (-this.D * q.z - this.z * -q.D) * denom;
             z0 = (this.y * -q.D + this.D * q.y) * denom;
         } else if (Math.abs(denomXZ) >= Math.abs(denomXY) && Math.abs(denomXZ) >= Math.abs(denomYZ)) {
             // Y is the biggest, so y0 = 0.0
-            if (Math.abs(denomXZ) < 1.0e-35)
+            if (Math.abs(denomXZ) < MINIMUM_RESOLUTION_SQUARED) {
+                //System.err.println(" Denominator is zero: no intersection");
                 return NO_POINTS;
+            }
             final double denom = 1.0 / denomXZ;
             x0 = (-this.D * q.z - this.z * -q.D) * denom;
             y0 = 0.0;
             z0 = (this.x * -q.D + this.D * q.x) * denom;
         } else {
             // Z is the biggest, so Z0 = 0.0
-            if (Math.abs(denomXY) < 1.0e-35)
+            if (Math.abs(denomXY) < MINIMUM_RESOLUTION_SQUARED) {
+                //System.err.println(" Denominator is zero: no intersection");
                 return NO_POINTS;
+            }
             final double denom = 1.0 / denomXY;
             x0 = (-this.D * q.y - this.y * -q.D) * denom;
             y0 = (this.x * -q.D + this.D * q.x) * denom;
@@ -159,23 +197,25 @@ public class Plane extends Vector
         final double C = x0*x0 + y0*y0 + z0*z0 - 1.0;
 
         final double BsquaredMinus = B * B - 4.0 * A * C;
-        if (BsquaredMinus < 0.0)
+        if (Math.abs(BsquaredMinus) < MINIMUM_RESOLUTION_SQUARED) {
-            return NO_POINTS;
+            //System.err.println(" One point of intersection");
-        final double inverse2A = 1.0 / (2.0 * A);
+            final double inverse2A = 1.0 / (2.0 * A);
-        if (BsquaredMinus == 0.0) {
             // One solution only
             final double t = -B * inverse2A;
             GeoPoint point = new GeoPoint(lineVector.x * t + x0, lineVector.y * t + y0, lineVector.z * t + z0);
             if (point.isWithin(bounds,moreBounds))
                 return new GeoPoint[]{point};
             return NO_POINTS;
-        } else {
+        } else if (BsquaredMinus > 0.0) {
+            //System.err.println(" Two points of intersection");
+            final double inverse2A = 1.0 / (2.0 * A);
             // Two solutions
             final double sqrtTerm = Math.sqrt(BsquaredMinus);
             final double t1 = (-B + sqrtTerm) * inverse2A;
             final double t2 = (-B - sqrtTerm) * inverse2A;
             GeoPoint point1 = new GeoPoint(lineVector.x * t1 + x0, lineVector.y * t1 + y0, lineVector.z * t1 + z0);
             GeoPoint point2 = new GeoPoint(lineVector.x * t2 + x0, lineVector.y * t2 + y0, lineVector.z * t2 + z0);
+            //System.err.println("  "+point1+" and "+point2);
             if (point1.isWithin(bounds,moreBounds)) {
                 if (point2.isWithin(bounds,moreBounds))
                     return new GeoPoint[]{point1,point2};
@@ -184,6 +224,9 @@ public class Plane extends Vector
             if (point2.isWithin(bounds,moreBounds))
                 return new GeoPoint[]{point2};
             return NO_POINTS;
+        } else {
+            //System.err.println(" no solutions - no intersection");
+            return NO_POINTS;
         }
     }
     
@@ -210,12 +253,13 @@ public class Plane extends Vector
     */
     public void recordBounds(final Bounds boundsInfo, final Membership... bounds) {
         // For clarity, load local variables with good names
-        double A = this.x;
+        final double A = this.x;
-        double B = this.y;
+        final double B = this.y;
-        double C = this.z;
+        final double C = this.z;
 
         // Now compute latitude min/max points
         if (!boundsInfo.checkNoTopLatitudeBound() || !boundsInfo.checkNoBottomLatitudeBound()) {
+            //System.err.println("Looking at latitude for plane "+this);
             if ((Math.abs(A) >= MINIMUM_RESOLUTION || Math.abs(B) >= MINIMUM_RESOLUTION)) {
                 //System.out.println("A = "+A+" B = "+B+" C = "+C+" D = "+D);
                 // sin (phi) = z
@@ -225,9 +269,9 @@ public class Plane extends Vector
                 //
                 // cos (theta-phi) = cos(theta)cos(phi) + sin(theta)sin(phi) = D
 
-                if (Math.abs(C) < 1.0e-10) {
+                if (Math.abs(C) < MINIMUM_RESOLUTION) {
                     // Special case: circle is vertical.
-                    //System.out.println("Degenerate case; it's vertical circle");
+                    //System.err.println(" Degenerate case; it's vertical circle");
                     // cos(phi) = D, and we want sin(phi) = z
                     // There are two solutions for phi given cos(phi) = D: a positive solution and a negative solution.
                     // So, when we compute z = sqrt(1-D^2), it's really z = +/- sqrt(1-D^2) .
@@ -236,7 +280,7 @@ public class Plane extends Vector
                     double x;
                     double y;
 
-                    double denom = 1.0 / (A*A + B*B);
+                    final double denom = 1.0 / (A*A + B*B);
 
                     z = Math.sqrt(1.0 - D*D);
                     y = -B * D * denom;
@@ -245,7 +289,34 @@ public class Plane extends Vector
 
                     z = -z;
                     addPoint(boundsInfo, bounds, x, y, z);
+                } else if (Math.abs(D) < MINIMUM_RESOLUTION) {
+                    //System.err.println(" Plane through origin case");
+                    // The general case is degenerate when the plane goes through the origin.
+                    // Luckily there's a pretty good way to figure out the max and min for that case though.
+                    // We find the two z values by computing the angle of the plane's inclination with the normal.
+                    // E.g., if this.z == 1, then our z value is 0, and if this.z == 0, our z value is 1.
+                    // Also if this.z == -1, then z value is 0 again.
+                    // Another way of putting this is that our z = sqrt(this.x^2 + this.y^2).
+                    //
+                    // The only tricky part is computing x and y.
+                    double z;
+                    double x;
+                    double y;
+                    
+                    final double denom = 1.0 / (A*A + B*B);
+                    
+                    z = Math.sqrt((A * A + B * B)/(A*A+B*B+C*C));
+                    y = -B * (C*z) * denom;
+                    x = -A * (C*z) * denom;
+                    addPoint(boundsInfo, bounds, x, y, z);
+                    
+                    z = -z;
+                    y = -B * (C*z) * denom;
+                    x = -A * (C*z) * denom;
+                    addPoint(boundsInfo, bounds, x, y, z);
+
                 } else {
+                    //System.err.println(" General latitude case");
                     // We might be able to identify a specific new latitude maximum or minimum.
                     //
                     // cos (theta-phi) = cos(theta)cos(phi) + sin(theta)sin(phi) = D
@@ -261,80 +332,153 @@ public class Plane extends Vector
                     //
                     // D = cos(theta)cos(phi) + sin(theta)sin(phi)
                     // Substitute:
-                    // D = sqrt(1-C^2) * sqrt(1-z^2) + C * z
+                    // D = sqrt(1-C^2) * sqrt(1-z^2) -/+ C * z
                     // Solve for z...
-                    // D-Cz = sqrt(1-C^2)*sqrt(1-z^2) = sqrt(1 - z^2 - C^2 + z^2*C^2)
+                    // D +/- Cz = sqrt(1-C^2)*sqrt(1-z^2) = sqrt(1 - z^2 - C^2 + z^2*C^2)
                     // Square both sides.
-                    // (D-Cz)^2 = 1 - z^2 - C^2 + z^2*C^2
+                    // (D +/- Cz)^2 = 1 - z^2 - C^2 + z^2*C^2
-                    // D^2 - 2DCz + C^2*z^2 = 1 - z^2 - C^2 + z^2*C^2
+                    // D^2 +/- 2DCz + C^2*z^2 = 1 - z^2 - C^2 + z^2*C^2
-                    // D^2 - 2DCz  = 1 - C^2 - z^2
+                    // D^2 +/- 2DCz  = 1 - C^2 - z^2
-                    // 0 = z^2 - 2DCz + (C^2 +D^2-1) = 0
+                    // 0 = z^2 +/- 2DCz + (C^2 +D^2-1) = 0
                     //
-                    // z = (2DC +/- sqrt(4*D^2*C^2 - 4*(C^2+D^2-1))) / (2)
+                    // z = (+/- 2DC +/- sqrt(4*D^2*C^2 - 4*(C^2+D^2-1))) / (2)
-                    // z  = DC +/- sqrt(D^2*C^2 + 1 - C^2 - D^2 )
+                    // z  = +/- DC +/- sqrt(D^2*C^2 + 1 - C^2 - D^2 )
-                    //    = DC +/- sqrt(D^2*C^2 + 1 - C^2 - D^2)
+                    //    = +/- DC +/- sqrt(D^2*C^2 + 1 - C^2 - D^2)
-
+                    //
+                    // NOTE WELL: The above is clearly degenerate when D = 0.  So we'll have to
+                    // code a different solution for that case!
+                    
+                    // To get x and y, we need to plug z into the equations, as follows:
+                    //
+                    // Ax + By = -Cz - D
+                    // x^2 + y^2 = 1 - z^2
+                    //
+                    // x = (-Cz -D -By) /A
+                    // y = (-Cz -D -Ax) /B
+                    //
+                    // [(-Cz -D -By) /A]^2 + y^2 = 1 - z^2
+                    // [-Cz -D -By]^2 + A^2*y^2 = A^2 - A^2*z^2
+                    // C^2*z^2 + D^2 + B^2*y^2 + 2CDz + 2CBzy + 2DBy + A^2*y^2 - A^2 + A^2*z^2 = 0
+                    // y^2 [A^2 + B^2]  + y [2DB + 2CBz] + [C^2*z^2 + D^2 + 2CDz - A^2 + A^2*z^2] = 0
+                    //
+                    // 
+                    // Use quadratic formula, where:
+                    // a = [A^2 + B^2]
+                    // b = [2BD + 2CBz]
+                    // c = [C^2*z^2 + D^2 + 2CDz - A^2 + A^2*z^2]
+                    //
+                    // y = (-[2BD + 2CBz] +/- sqrt([2BD + 2CBz]^2 - 4 * [A^2 + B^2] * [C^2*z^2 + D^2 + 2CDz - A^2 + A^2*z^2]) ) / (2 * [A^2 + B^2])
+                    // Take out a 2:
+                    // y = (-[DB +CBz] +/- sqrt([DB + CBz]^2 - [A^2 + B^2] * [C^2*z^2 + D^2 + 2CDz - A^2 + A^2*z^2]) ) / [A^2 + B^2]
+                    //
+                    // The sqrt term simplifies:
+                    //
+                    // B^2*D^2 + C^2*B^2*z^2 + 2C*D*B^2*z - [A^2 + B^2] * [C^2*z^2 + D^2 + 2CDz - A^2 + A^2*z^2] = ?
+                    // B^2*D^2 + C^2*B^2*z^2 + 2C*D*B^2*z - [A^2 * C^2 * z^2 + A^2 * D^2 + 2 * A^2 * CDz - A^4 + A^4*z^2
+                    //                  + B^2 * C^2 * z^2 + B^2 * D^2 + 2 * B^2 * CDz - A^2 * B^2 + B^2 * A^2 * z^2] =?
+                    // C^2*B^2*z^2 + 2C*D*B^2*z - [A^2 * C^2 * z^2 + A^2 * D^2 + 2 * A^2 * CDz - A^4 + A^4*z^2
+                    //                  + B^2 * C^2 * z^2 + 2 * B^2 * CDz - A^2 * B^2 + B^2 * A^2 * z^2] =?
+                    // 2C*D*B^2*z - [A^2 * C^2 * z^2 + A^2 * D^2 + 2 * A^2 * CDz - A^4 + A^4*z^2
+                    //                  + 2 * B^2 * CDz - A^2 * B^2 + B^2 * A^2 * z^2] =?
+                    // - [A^2 * C^2 * z^2 + A^2 * D^2 + 2 * A^2 * CDz - A^4 + A^4*z^2
+                    //                  - A^2 * B^2 + B^2 * A^2 * z^2] =?
+                    // - A^2 * [C^2 * z^2 + D^2 + 2 * CDz - A^2 + A^2*z^2
+                    //                  - B^2 + B^2 * z^2] =?
+                    // - A^2 * [z^2[A^2 + B^2 + C^2] - [A^2 + B^2 - D^2] + 2CDz] =?
+                    // - A^2 * [z^2 - [A^2 + B^2 - D^2] + 2CDz] =?
+                    //
+                    // y = (-[DB +CBz] +/- A*sqrt([A^2 + B^2 - D^2] - z^2 - 2CDz) ) / [A^2 + B^2]
+                    //
+                    // correspondingly:
+                    // x = (-[DA +CAz] +/- B*sqrt([A^2 + B^2 - D^2] - z^2 - 2CDz) ) / [A^2 + B^2]
+                    //
+                    // However, for the maximum or minimum we seek, the clause inside the sqrt should be zero.  If
+                    // it is NOT zero, then we aren't looking at the right z value.
+                    
                     double z;
                     double x;
                     double y;
                     
                     double sqrtValue = D*D*C*C + 1.0 - C*C - D*D;
-                    if (sqrtValue >= 0.0) {
+                    double denom = 1.0 / (A*A + B*B);
-                        // y = -B[D+Cz] / [A^2 + B^2]
+                    if (Math.abs(sqrtValue) < MINIMUM_RESOLUTION_SQUARED) {
-                        // x = -A[D+Cz] / [A^2 + B^2]
+                        //System.err.println(" One latitude solution");
-                        double denom = 1.0 / (A*A + B*B);
+                        double insideValue;
-                        if (sqrtValue == 0.0) {
+                        double sqrtTerm;
-                            //System.out.println("Zero sqrt term");
+                            
-                            z = D*C;
+                        z = D*C;
-                            // Since we squared both sides of the equation, we may have introduced spurious solutions, so we have to check.
+                        // Since we squared both sides of the equation, we may have introduced spurious solutions, so we have to check.
-                            if (Math.abs(D-C*z - Math.sqrt(1.0 - z*z - C*C + z*z*C*C)) < 1.0e-10) {
+                        // But the same check applies to BOTH solutions -- the +z one as well as the -z one.
-                                y = -B * (D + C*z) * denom;
+                        insideValue = A * A + B * B - D * D - z*z - 2.0 * C * D * z;
-                                x = -A * (D + C*z) * denom;
+                        if (Math.abs(insideValue) < MINIMUM_RESOLUTION) {
+                            y = -B * (D + C*z) * denom;
+                            x = -A * (D + C*z) * denom;
+                            if (evaluateIsZero(x,y,z)) {
                                 addPoint(boundsInfo, bounds, x, y, z);
                             }
-                            z = -D*C;
+                        }
-                            // Since we squared both sides of the equation, we may have introduced spurious solutions, so we have to check.
+                        // Check the solution on the other side of the x-y plane
-                            if (Math.abs(D+C*z + Math.sqrt(1.0 - z*z - C*C + z*z*C*C)) < 1.0e-10) {
+                        z = -z;
-                                y = -B * (D + C*z) * denom;
+                        insideValue = A * A + B * B - D * D - z*z - 2.0 * C * D * z;
-                                x = -A * (D + C*z) * denom;
+                        if (Math.abs(insideValue) < MINIMUM_RESOLUTION) {
+                            y = -B * (D + C*z) * denom;
+                            x = -A * (D + C*z) * denom;
+                            if (evaluateIsZero(x,y,z)) {
                                 addPoint(boundsInfo, bounds, x, y, z);
                             }
-                        } else {
+                        }
-                            double sqrtResult = Math.sqrt(sqrtValue);
+                    } else if (sqrtValue > 0.0) {
-                            z = D*C + sqrtResult;
+                        //System.err.println(" Two latitude solutions");
-                            //System.out.println("z= "+z+" D-C*z = " + (D-C*z) + " Math.sqrt(1.0 - z*z - C*C + z*z*C*C) = "+(Math.sqrt(1.0 - z*z - C*C + z*z*C*C)));
+                        double sqrtResult = Math.sqrt(sqrtValue);
-                            // Since we squared both sides of the equation, we may have introduced spurios solutions, so we have to check.
+                            
-                            if (Math.abs(D-C*z - Math.sqrt(1.0 - z*z - C*C + z*z*C*C)) < 1.0e-10) {
+                        double insideValue;
-                                //System.out.println("found a point; z = "+z);
+                        double sqrtTerm;
-                                y = -B * (D + C*z) * denom;
+                            
-                                x = -A * (D + C*z) * denom;
+                        z = D*C + sqrtResult;
+                        //System.out.println("z= "+z+" D-C*z = " + (D-C*z) + " Math.sqrt(1.0 - z*z - C*C + z*z*C*C) = "+(Math.sqrt(1.0 - z*z - C*C + z*z*C*C)));
+                        // Since we squared both sides of the equation, we may have introduced spurios solutions, so we have to check.
+                        // But the same check applies to BOTH solutions -- the +z one as well as the -z one.
+                        insideValue = A * A + B * B - D * D - z*z - 2.0 * C * D * z;
+                        //System.err.println(" z="+z+" C="+C+" D="+D+" inside value "+insideValue);
+                        if (Math.abs(insideValue) < MINIMUM_RESOLUTION) {
+                            y = -B * (D + C*z) * denom;
+                            x = -A * (D + C*z) * denom;
+                            if (evaluateIsZero(x,y,z)) {
                                 addPoint(boundsInfo, bounds, x, y, z);
                             }
-                            z = D*C - sqrtResult;
+                        }
-                            //System.out.println("z= "+z+" D-C*z = " + (D-C*z) + " Math.sqrt(1.0 - z*z - C*C + z*z*C*C) = "+(Math.sqrt(1.0 - z*z - C*C + z*z*C*C)));
+                        // Check the solution on the other side of the x-y plane
-                            // Since we squared both sides of the equation, we may have introduced spurios solutions, so we have to check.
+                        z = -z;
-                            if (Math.abs(D-C*z - Math.sqrt(1.0 - z*z - C*C + z*z*C*C)) < 1.0e-10) {
+                        insideValue = A * A + B * B - D * D - z*z - 2.0 * C * D * z;
-                                //System.out.println("found a point; z="+z);
+                        //System.err.println(" z="+z+" C="+C+" D="+D+" inside value "+insideValue);
-                                y = -B * (D + C*z) * denom;
+                        if (Math.abs(insideValue) < MINIMUM_RESOLUTION) {
-                                x = -A * (D + C*z) * denom;
+                            y = -B * (D + C*z) * denom;
+                            x = -A * (D + C*z) * denom;
+                            if (evaluateIsZero(x,y,z)) {
                                 addPoint(boundsInfo, bounds, x, y, z);
                             }
-                            z = -(D*C + sqrtResult);
+                        }
-                            //System.out.println("z= "+z+" D+C*z = " + (D+C*z) + " -Math.sqrt(1.0 - z*z - C*C + z*z*C*C) = "+(-Math.sqrt(1.0 - z*z - C*C + z*z*C*C)));
+                        z = D*C - sqrtResult;
-                            // Since we squared both sides of the equation, we may have introduced spurios solutions, so we have to check.
+                        //System.out.println("z= "+z+" D-C*z = " + (D-C*z) + " Math.sqrt(1.0 - z*z - C*C + z*z*C*C) = "+(Math.sqrt(1.0 - z*z - C*C + z*z*C*C)));
-                            if (Math.abs(D+C*z + Math.sqrt(1.0 - z*z - C*C + z*z*C*C)) < 1.0e-10) {
+                        // Since we squared both sides of the equation, we may have introduced spurios solutions, so we have to check.
-                                //System.out.println("found a point; z = "+z);
+                        // But the same check applies to BOTH solutions -- the +z one as well as the -z one.
-                                y = -B * (D + C*z) * denom;
+                        insideValue = A * A + B * B - D * D - z*z - 2.0 * C * D * z;
-                                x = -A * (D + C*z) * denom;
+                        //System.err.println(" z="+z+" C="+C+" D="+D+" inside value "+insideValue);
+                        if (Math.abs(insideValue) < MINIMUM_RESOLUTION) {
+                            y = -B * (D + C*z) * denom;
+                            x = -A * (D + C*z) * denom;
+                            if (evaluateIsZero(x,y,z)) {
                                 addPoint(boundsInfo, bounds, x, y, z);
                             }
-                            z = -(D*C - sqrtResult);
+                        }
-                            //System.out.println("z= "+z+" D+C*z = " + (D+C*z) + " -Math.sqrt(1.0 - z*z - C*C + z*z*C*C) = "+(-Math.sqrt(1.0 - z*z - C*C + z*z*C*C)));
+                        // Check the solution on the other side of the x-y plane
-                            // Since we squared both sides of the equation, we may have introduced spurios solutions, so we have to check.
+                        z = -z;
-                            if (Math.abs(D+C*z + Math.sqrt(1 - z*z - C*C + z*z*C*C)) < 1.0e-10) {
+                        insideValue = A * A + B * B - D * D - z*z - 2.0 * C * D * z;
-                                //System.out.println("found a point; z="+z);
+                        //System.err.println(" z="+z+" C="+C+" D="+D+" inside value "+insideValue);
-                                y = -B * (D + C*z) * denom;
+                        if (Math.abs(insideValue) < MINIMUM_RESOLUTION) {
-                                x = -A * (D + C*z) * denom;
+                            y = -B * (D + C*z) * denom;
+                            x = -A * (D + C*z) * denom;
+                            if (evaluateIsZero(x,y,z)) {
                                 addPoint(boundsInfo, bounds, x, y, z);
                             }
                         }
@@ -347,11 +491,12 @@ public class Plane extends Vector
                 // to check Membership objects.
                 boundsInfo.addHorizontalCircle(-D * C);
             }
+            //System.err.println("Done latitude bounds");
         }
         
         // First, figure out our longitude bounds, unless we no longer need to consider that
         if (!boundsInfo.checkNoLongitudeBound()) {
-            
+            //System.err.println("Computing longitude bounds for "+this);
             //System.out.println("A = "+A+" B = "+B+" C = "+C+" D = "+D);
             // Compute longitude bounds
             
@@ -396,29 +541,27 @@ public class Plane extends Vector
                         
                         double sqrtClause = b * b - 4.0 * a * c;
                             
-                        if (sqrtClause >= 0.0) {
+                        if (Math.abs(sqrtClause) < MINIMUM_RESOLUTION_SQUARED) {
-                            if (sqrtClause == 0.0) {
+                            double y0 = -b / (2.0 * a);
-                                double y0 = -b / (2.0 * a);
+                            double x0 = (-D - B * y0) / A;
-                                double x0 = (-D - B * y0) / A;
+                            double z0 = 0.0;
-                                double z0 = 0.0;
+                            addPoint(boundsInfo, bounds, x0, y0, z0);
-                                addPoint(boundsInfo, bounds, x0, y0, z0);
+                        } else if (sqrtClause > 0.0) {
-                            } else {
+                            double sqrtResult = Math.sqrt(sqrtClause);
-                                double sqrtResult = Math.sqrt(sqrtClause);
+                            double denom = 1.0 / (2.0 * a);
-                                double denom = 1.0 / (2.0 * a);
+                            double Hdenom = 1.0 / A;
-                                double Hdenom = 1.0 / A;
                                 
-                                double y0a = (-b + sqrtResult ) * denom;
+                            double y0a = (-b + sqrtResult ) * denom;
-                                double y0b = (-b - sqrtResult ) * denom;
+                            double y0b = (-b - sqrtResult ) * denom;
 
-                                double x0a = (-D - B * y0a) * Hdenom;
+                            double x0a = (-D - B * y0a) * Hdenom;
-                                double x0b = (-D - B * y0b) * Hdenom;
+                            double x0b = (-D - B * y0b) * Hdenom;
                                 
-                                double z0a = 0.0;
+                            double z0a = 0.0;
-                                double z0b = 0.0;
+                            double z0b = 0.0;
 
-                                addPoint(boundsInfo, bounds, x0a, y0a, z0a);
+                            addPoint(boundsInfo, bounds, x0a, y0a, z0a);
-                                addPoint(boundsInfo, bounds, x0b, y0b, z0b);
+                            addPoint(boundsInfo, bounds, x0b, y0b, z0b);
-                            }
                         }
 
                     } else {
@@ -431,33 +574,36 @@ public class Plane extends Vector
                         
                         double sqrtClause = b * b - 4.0 * a * c;
 
-                        if (sqrtClause >= 0.0) {
+                        if (Math.abs(sqrtClause) < MINIMUM_RESOLUTION_SQUARED) {
-                            if (sqrtClause == 0.0) {
+                            double x0 = -b / (2.0 * a);
-                                double x0 = -b / (2.0 * a);
+                            double y0 = (- D - A * x0) / B;
-                                double y0 = (- D - A * x0) / B;
+                            double z0 = 0.0;
-                                double z0 = 0.0;
+                            addPoint(boundsInfo, bounds, x0, y0, z0);
-                                addPoint(boundsInfo, bounds, x0, y0, z0);
+                        } else if (sqrtClause > 0.0) {
-                            } else {
+                            double sqrtResult = Math.sqrt(sqrtClause);
-                                double sqrtResult = Math.sqrt(sqrtClause);
+                            double denom = 1.0 / (2.0 * a);
-                                double denom = 1.0 / (2.0 * a);
+                            double Idenom = 1.0 / B;
-                                double Idenom = 1.0 / B;
                                 
-                                double x0a = (-b + sqrtResult ) * denom;
+                            double x0a = (-b + sqrtResult ) * denom;
-                                double x0b = (-b - sqrtResult ) * denom;
+                            double x0b = (-b - sqrtResult ) * denom;
-                                double y0a = (- D - A * x0a) * Idenom;
+                            double y0a = (- D - A * x0a) * Idenom;
-                                double y0b = (- D - A * x0b) * Idenom;
+                            double y0b = (- D - A * x0b) * Idenom;
-                                double z0a = 0.0;
+                            double z0a = 0.0;
-                                double z0b = 0.0;
+                            double z0b = 0.0;
 
-                                addPoint(boundsInfo, bounds, x0a, y0a, z0a);
+                            addPoint(boundsInfo, bounds, x0a, y0a, z0a);
-                                addPoint(boundsInfo, bounds, x0b, y0b, z0b);
+                            addPoint(boundsInfo, bounds, x0b, y0b, z0b);
-                            }
                         }
                     }
                 }
                 
             } else {
-
+                //System.err.println("General longitude bounds...");
+                
+                // NOTE WELL: The x,y,z values generated here are NOT on the unit sphere.
+                // They are for lat/lon calculation purposes only.  x-y is meant to be used for longitude determination,
+                // and z for latitude, and that's all the values are good for.
+                
                 // (1) Intersect the plane and the unit sphere, and project the results into the x-y plane:
                 // From plane:
                 // z = (-Ax - By - D) / C
@@ -478,10 +624,15 @@ public class Plane extends Vector
                 double I = 2.0 * B * D;
                 double J = D * D - C * C;
 
-                //System.out.println("E = " + E + " F = " + F + " G = " + G + " H = "+ H + " I = " + I + " J = " + J);
+                //System.err.println("E = " + E + " F = " + F + " G = " + G + " H = "+ H + " I = " + I + " J = " + J);
+                
+                double trialX = 2.0;
+                double trialY = 2.0;
+                
+                //System.err.println("Trial point evaluates to: "+(E*trialX*trialX + F*trialY*trialY + G*trialX*trialY + H*trialX + I*trialY + J));
                 
                 // Check if the origin is within, by substituting x = 0, y = 0 and seeing if less than zero
-                if (J > 0.0) {
+                if (Math.abs(J) >= MINIMUM_RESOLUTION && J > 0.0) {
                     // The derivative of the curve above is:
                     // 2Exdx + 2Fydy + G(xdy+ydx) + Hdx + Idy = 0
                     // (2Ex + Gy + H)dx + (2Fy + Gx + I)dy = 0
@@ -503,7 +654,7 @@ public class Plane extends Vector
                     // But we will need to base this on which coefficient is non-zero
                     
                     if (Math.abs(H) > Math.abs(I)) {
-                        //System.out.println("Using the y quadratic");
+                        //System.err.println(" Using the y quadratic");
                         // x = (-2J - Iy)/H
                         
                         // Plug into the original equation:
@@ -534,35 +685,33 @@ public class Plane extends Vector
                         double sqrtClause = b * b - 4.0 * a * c;
                         //System.out.println("sqrtClause="+sqrtClause);
                         
-                        if (sqrtClause >= 0.0) {
+                        if (Math.abs(sqrtClause) < MINIMUM_RESOLUTION_SQUARED) {
-                            if (sqrtClause == 0.0) {
+                            //System.err.println(" One solution");
-                                //System.out.println("One solution");
+                            double y0 = -b / (2.0 * a);
-                                double y0 = -b / (2.0 * a);
+                            double x0 = (-2.0 * J - I * y0) / H;
-                                double x0 = (-2.0 * J - I * y0) / H;
+                            double z0 = (-A*x0 - B*y0 - D)/C;
-                                double z0 = (-A*x0 - B*y0 - D)/C;
 
-                                addPoint(boundsInfo, bounds, x0, y0, z0);
+                            addPoint(boundsInfo, bounds, x0, y0, z0);
-                            } else {
+                        } else if (sqrtClause > 0.0) {
-                                //System.out.println("Two solutions");
+                            //System.err.println(" Two solutions");
-                                double sqrtResult = Math.sqrt(sqrtClause);
+                            double sqrtResult = Math.sqrt(sqrtClause);
-                                double denom = 1.0 / (2.0 * a);
+                            double denom = 1.0 / (2.0 * a);
-                                double Hdenom = 1.0 / H;
+                            double Hdenom = 1.0 / H;
-                                double Cdenom = 1.0 / C;
+                            double Cdenom = 1.0 / C;
                                     
-                                double y0a = (-b + sqrtResult ) * denom;
+                            double y0a = (-b + sqrtResult ) * denom;
-                                double y0b = (-b - sqrtResult ) * denom;
+                            double y0b = (-b - sqrtResult ) * denom;
-                                double x0a = (-2.0 * J - I * y0a) * Hdenom;
+                            double x0a = (-2.0 * J - I * y0a) * Hdenom;
-                                double x0b = (-2.0 * J - I * y0b) * Hdenom;
+                            double x0b = (-2.0 * J - I * y0b) * Hdenom;
-                                double z0a = (-A*x0a - B*y0a - D) * Cdenom;
+                            double z0a = (-A*x0a - B*y0a - D) * Cdenom;
-                                double z0b = (-A*x0b - B*y0b - D) * Cdenom;
+                            double z0b = (-A*x0b - B*y0b - D) * Cdenom;
 
-                                addPoint(boundsInfo, bounds, x0a, y0a, z0a);
+                            addPoint(boundsInfo, bounds, x0a, y0a, z0a);
-                                addPoint(boundsInfo, bounds, x0b, y0b, z0b);
+                            addPoint(boundsInfo, bounds, x0b, y0b, z0b);
-                            }
                         }
 
                     } else {
-                        //System.out.println("Using the x quadratic");
+                        //System.err.println(" Using the x quadratic");
                         // y = (-2J - Hx)/I
                         
                         // Plug into the original equation:
@@ -578,38 +727,39 @@ public class Plane extends Vector
                         // Group:
                         // x^2 [E*I^2 - GHI + F*H^2] + x [4FJH - 2GIJ] + [4FJ^2 - J*I^2] = 0
 
+                        // E x^2 + F y^2 + G xy + H x + I y + J = 0
+
                         a = E * I * I - G * H * I + F*H*H;
                         b = 4.0 * F * H * J - 2.0 * G * I * J;
                         c = 4.0 * F * J * J - J * I * I;
-                        
+
                         //System.out.println("a="+a+" b="+b+" c="+c);
                         double sqrtClause = b * b - 4.0 * a * c;
                         //System.out.println("sqrtClause="+sqrtClause);
-                        
+                        if (Math.abs(sqrtClause) < MINIMUM_RESOLUTION_SQUARED) {
-                        if (sqrtClause >= 0.0) {
+                            //System.err.println(" One solution; sqrt clause was "+sqrtClause);
-                            if (sqrtClause == 0.0) {
+                            double x0 = -b / (2.0 * a);
-                                //System.out.println("One solution");
+                            double y0 = (-2.0 * J - H * x0) / I;
-                                double x0 = -b / (2.0 * a);
+                            double z0 = (-A*x0 - B*y0 - D)/C;
-                                double y0 = (-2.0 * J - H * x0) / I;
+                            // Verify that x&y fulfill the equation
-                                double z0 = (-A*x0 - B*y0 - D)/C;
+                            // 2Ex^2 + 2Fy^2 + 2Gxy + Hx + Iy = 0
-                                addPoint(boundsInfo, bounds, x0, y0, z0);
+                            addPoint(boundsInfo, bounds, x0, y0, z0);
-                            } else {
+                        } else if (sqrtClause > 0.0) {
-                                //System.out.println("Two solutions");
+                            //System.err.println(" Two solutions");
-                                double sqrtResult = Math.sqrt(sqrtClause);
+                            double sqrtResult = Math.sqrt(sqrtClause);
-                                double denom = 1.0 / (2.0 * a);
+                            double denom = 1.0 / (2.0 * a);
-                                double Idenom = 1.0 / I;
+                            double Idenom = 1.0 / I;
-                                double Cdenom = 1.0 / C;
+                            double Cdenom = 1.0 / C;
                                     
-                                double x0a = (-b + sqrtResult ) * denom;
+                            double x0a = (-b + sqrtResult ) * denom;
-                                double x0b = (-b - sqrtResult ) * denom;
+                            double x0b = (-b - sqrtResult ) * denom;
-                                double y0a = (-2.0 * J - H * x0a) * Idenom;
+                            double y0a = (-2.0 * J - H * x0a) * Idenom;
-                                double y0b = (-2.0 * J - H * x0b) * Idenom;
+                            double y0b = (-2.0 * J - H * x0b) * Idenom;
-                                double z0a = (-A*x0a - B*y0a - D) * Cdenom;
+                            double z0a = (-A*x0a - B*y0a - D) * Cdenom;
-                                double z0b = (-A*x0b - B*y0b - D) * Cdenom;
+                            double z0b = (-A*x0b - B*y0b - D) * Cdenom;
 
-                                addPoint(boundsInfo, bounds, x0a, y0a, z0a);
+                            addPoint(boundsInfo, bounds, x0a, y0a, z0a);
-                                addPoint(boundsInfo, bounds, x0b, y0b, z0b);
+                            addPoint(boundsInfo, bounds, x0b, y0b, z0b);
-                            }
                         }
                     }
                 }
@@ -619,12 +769,14 @@ public class Plane extends Vector
     }
     
     protected static void addPoint(final Bounds boundsInfo, final Membership[] bounds, final double x, final double y, final double z) {
+        //System.err.println(" Want to add point x="+x+" y="+y+" z="+z);
         // Make sure the discovered point is within the bounds
         for (Membership bound : bounds) {
             if (!bound.isWithin(x,y,z))
                 return;
         }
         // Add the point
+        //System.err.println("  point added");
         //System.out.println("Adding point x="+x+" y="+y+" z="+z);
         boundsInfo.addPoint(x,y,z);
     }
@@ -639,26 +791,50 @@ public class Plane extends Vector
      *@return true if there's an intersection.
      */
     public boolean intersects(final Plane q, final GeoPoint[] notablePoints, final GeoPoint[] moreNotablePoints, final Membership[] bounds, final Membership... moreBounds) {
+        //System.err.println("Does plane "+this+" intersect with plane "+q);
         // If the two planes are identical, then the math will find no points of intersection.
         // So a special case of this is to check for plane equality.  But that is not enough, because
         // what we really need at that point is to determine whether overlap occurs between the two parts of the intersection
         // of plane and circle.  That is, are there *any* points on the plane that are within the bounds described?
-	if (equals(q)) {
+	if (isNumericallyIdentical(q)) {
+            //System.err.println(" Identical plane");
             // The only way to efficiently figure this out will be to have a list of trial points available to evaluate.
             // We look for any point that fulfills all the bounds.
             for (GeoPoint p : notablePoints) {
-                if (meetsAllBounds(p,bounds,moreBounds))
+                if (meetsAllBounds(p,bounds,moreBounds)) {
+                    //System.err.println("  found a notable point in bounds, so intersects");
                     return true;
+                }
             }
             for (GeoPoint p : moreNotablePoints) {
-                if (meetsAllBounds(p,bounds,moreBounds))
+                if (meetsAllBounds(p,bounds,moreBounds)) {
+                    //System.err.println("  found a notable point in bounds, so intersects");
                     return true;
+                }
             }
+            //System.err.println("  no notable points inside found; no intersection");
             return false;
         }
         return findIntersections(q,bounds,moreBounds).length > 0;
     }
 
+    /** Returns true if this plane and the other plane are identical within the margin of error.
+    */
+    protected boolean isNumericallyIdentical(final Plane p) {
+        // We can get the correlation by just doing a parallel plane check.  If that passes, then compute a point on the plane
+        // (using D) and see if it also on the other plane.
+        if (Math.abs(this.y * p.z - this.z * p.y) >= MINIMUM_RESOLUTION_SQUARED)
+            return false;
+        if (Math.abs(this.z * p.x - this.x * p.z) >= MINIMUM_RESOLUTION_SQUARED)
+            return false;
+        if (Math.abs(this.x * p.y - this.y * p.x) >= MINIMUM_RESOLUTION_SQUARED)
+            return false;
+
+        // Now, see whether the parallel planes are in fact on top of one another.  
+        final double denom = 1.0 / (p.x * p.x + p.y * p.y + p.z * p.z);
+        return evaluateIsZero(- p.x * p.D * denom, - p.y * p.D * denom, - p.z * p.D * denom);
+    }
+
     protected static boolean meetsAllBounds(final GeoPoint p, final Membership[] bounds, final Membership[] moreBounds) {
         for (final Membership bound : bounds) {
             if (!bound.isWithin(p))

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/SidedPlane.java

@@ -95,7 +95,7 @@ public class SidedPlane extends Plane implements Membership
     @Override
     public boolean isWithin(double x, double y, double z)
     {
-        double evalResult = this.x * x + this.y * y + this.z * z;
+        double evalResult = evaluate(x,y,z);
         if (Math.abs(evalResult) < MINIMUM_RESOLUTION)
             return true;
         double sigNum = Math.signum(evalResult);

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/Vector.java

@@ -23,7 +23,10 @@ public class Vector
 {
     /** Values that are all considered to be essentially zero have a magnitude
     * less than this. */
-    public static final double MINIMUM_RESOLUTION = 1e-15;
+    public static final double MINIMUM_RESOLUTION = 1e-12;
+    /** For squared quantities, the bound is squared too.
+    */
+    public static final double MINIMUM_RESOLUTION_SQUARED = MINIMUM_RESOLUTION * MINIMUM_RESOLUTION;
     
     public final double x;
     public final double y;
@@ -179,7 +182,7 @@ public class Vector
      *@return the square of the normal distance.
      */
     public double normalDistanceSquared(final Vector v) {
-        double t = this.evaluate(v);
+        double t = x*v.x + y*v.y + z*v.z;
         double deltaX = this.x * t - v.x;
         double deltaY = this.y * t - v.y;
         double deltaZ = this.z * t - v.z;
@@ -195,7 +198,7 @@ public class Vector
      *@return the square of the normal distance.
      */
     public double normalDistanceSquared(final double x, final double y, final double z) {
-        double t = this.evaluate(x,y,z);
+        double t = this.x*x + this.y*y + this.z*z;
         double deltaX = this.x * t - x;
         double deltaY = this.y * t - y;
         double deltaZ = this.z * t - z;

lucene/spatial/src/test/org/apache/lucene/spatial/spatial4j/Geo3dRptTest.java

@@ -24,6 +24,7 @@ import java.util.List;
 import com.carrotsearch.randomizedtesting.annotations.Repeat;
 import com.spatial4j.core.context.SpatialContext;
 import com.spatial4j.core.shape.Point;
+import com.spatial4j.core.shape.Rectangle;
 import com.spatial4j.core.shape.Shape;
 import org.apache.lucene.spatial.composite.CompositeSpatialStrategy;
 import org.apache.lucene.spatial.prefix.RandomSpatialOpStrategyTestCase;
@@ -76,6 +77,36 @@ public class Geo3dRptTest extends RandomSpatialOpStrategyTestCase {
         rptStrategy, serializedDVStrategy);
   }
 
+  @Test
+  public void testFailure() throws IOException {
+    setupStrategy();
+    // [junit4]    > Throwable #1: java.lang.AssertionError: [Intersects] qIdx:25 Shouldn't match I#1:Rect(minX=-49.0,maxX=-45.0,minY=73.0,maxY=86.0) 
+    // Q:Geo3dShape{GeoCompositeMembershipShape: {[GeoCompositeMembershipShape: {[GeoConvexPolygon: {[
+    // [X=-0.8606462131055999, Y=0.4385211485883089, Z=-0.25881904510252074], 
+    //  [X=-0.4668467715008339, Y=0.28050984011500923, Z=-0.838670567945424],
+    // [X=-0.9702957262759965, Y=1.1882695554102554E-16, Z=0.24192189559966773]]}]}, 
+    // GeoConvexPolygon: {[[X=0.8473975608908426, Y=-0.43177062311338915, Z=0.3090169943749474],
+    //[X=-0.4668467715008339, Y=0.28050984011500923, Z=-0.838670567945424], 
+    // [X=-0.8606462131055999, Y=0.4385211485883089, Z=-0.25881904510252074]]}]}}
+    //
+    // Points in order (I think):
+    // 
+    // [X=0.8473975608908426, Y=-0.43177062311338915, Z=0.3090169943749474],
+    //[X=-0.4668467715008339, Y=0.28050984011500923, Z=-0.838670567945424], 
+    // [X=-0.9702957262759965, Y=1.1882695554102554E-16, Z=0.24192189559966773],
+    // [X=-0.8606462131055999, Y=0.4385211485883089, Z=-0.25881904510252074], 
+    // Index: 0
+    final List<GeoPoint> points = new ArrayList<GeoPoint>();
+    points.add(new GeoPoint(18 * DEGREES_TO_RADIANS, -27 * DEGREES_TO_RADIANS));
+    points.add(new GeoPoint(-57 * DEGREES_TO_RADIANS, 146 * DEGREES_TO_RADIANS));
+    points.add(new GeoPoint(14 * DEGREES_TO_RADIANS, -180 * DEGREES_TO_RADIANS));
+    points.add(new GeoPoint(-15 * DEGREES_TO_RADIANS, 153 * DEGREES_TO_RADIANS));
+    
+    final Shape triangle = new Geo3dShape(GeoPolygonFactory.makeGeoPolygon(points,0),ctx);
+    final Rectangle rect = ctx.makeRectangle(-49, -45, 73, 86);
+    testOperation(rect,SpatialOperation.Intersects,triangle, false);
+  }
+
   @Test
   @Repeat(iterations = 2000)
   //@Seed("B808B88D6F8E285C")

lucene/spatial/src/test/org/apache/lucene/spatial/spatial4j/Geo3dShapeRectRelationTest.java

@@ -25,6 +25,9 @@ import com.carrotsearch.randomizedtesting.RandomizedContext;
 import com.spatial4j.core.context.SpatialContext;
 import com.spatial4j.core.distance.DistanceUtils;
 import com.spatial4j.core.shape.Point;
+import org.apache.lucene.spatial.spatial4j.geo3d.Bounds;
+import org.apache.lucene.spatial.spatial4j.geo3d.GeoArea;
+import org.apache.lucene.spatial.spatial4j.geo3d.GeoBBox;
 import org.apache.lucene.spatial.spatial4j.geo3d.GeoBBoxFactory;
 import org.apache.lucene.spatial.spatial4j.geo3d.GeoCircle;
 import org.apache.lucene.spatial.spatial4j.geo3d.GeoPath;
@@ -48,6 +51,73 @@ public class Geo3dShapeRectRelationTest extends RandomizedShapeTest {
     ctx = SpatialContext.GEO;
   }
 
+  protected final static double RADIANS_PER_DEGREE = Math.PI/180.0;
+  
+  @Test
+  public void testFailure() {
+      /*
+   [junit4]   1> S-R Rel: {}, Shape {}, Rectangle {} [INTERSECTS, Geo3dShape{GeoCompositeMembershipShape: {[
+   GeoConvexPolygon: {
+   points=[
+     [X=0.03206699943821901, Y=-0.7556330442094724, Z=0.6542097599743943], 
+     [X=-0.2848733212046893, Y=-0.9533780638748927, Z=0.09958643576296423], 
+     [X=0.37929990916639644, Y=0.9241954620264722, Z=0.044657887053005746]] 
+   edges={
+     [A=0.5484584327149066, B=-0.18956034526809354, C=-0.2458316687546487, D=0.0, side=1.0] internal? false; 
+     [A=-0.13461318190686059, B=0.05049496664187115, C=0.09833758231919826, D=0.0, side=1.0] internal? false; 
+     [A=0.6383626665235883, B=-0.246709658095017, C=-0.31624772039338794, D=0.0, side=1.0] internal? false; }}]}},
+     X=0.03206699943821901, Y=-0.7556330442094724, Z=0.6542097599743943
+   Rect(minX=-52.0,maxX=50.0,minY=58.0,maxY=68.0)](no slf4j subst; sorry)
+
+    */
+      final GeoBBox rect = GeoBBoxFactory.makeGeoBBox(68 * RADIANS_PER_DEGREE, 58 * RADIANS_PER_DEGREE, -52 * RADIANS_PER_DEGREE, 50 * RADIANS_PER_DEGREE);
+      final List<GeoPoint> points = new ArrayList<GeoPoint>();
+      points.add(new GeoPoint(40.8597568993 * RADIANS_PER_DEGREE, -87.5699819016 * RADIANS_PER_DEGREE));
+      points.add(new GeoPoint(5.71535611517 * RADIANS_PER_DEGREE, -106.636363741 * RADIANS_PER_DEGREE));
+      points.add(new GeoPoint(2.55955969779 * RADIANS_PER_DEGREE, 67.6862179901 * RADIANS_PER_DEGREE));
+      final GeoShape path = GeoPolygonFactory.makeGeoPolygon(points,0);
+      
+      System.err.println("Rectangle = "+rect+"; path = "+path);
+
+      // Edges intersect == OVERLAP.  This seems reasonable... between points 2 and 3 the path could well cross.
+      assertFalse(GeoArea.DISJOINT == rect.getRelationship(path));
+      
+      final GeoBBox pathBounds = getBoundingBox(path);
+      // Path bounds go around the back side of the world rather than the front.  The actual path goes around the front.  This is I think what the problem is.
+      System.err.println("Path bounds = "+pathBounds);
+      assertFalse(GeoArea.DISJOINT == rect.getRelationship(pathBounds));
+      
+      final GeoBBox rectBounds = getBoundingBox(rect);
+      assertFalse(GeoArea.DISJOINT == rectBounds.getRelationship(pathBounds));
+  }
+
+  protected static GeoBBox getBoundingBox(final GeoShape path) {
+      Bounds bounds = path.getBounds(null);
+
+      double leftLon;
+      double rightLon;
+      if (bounds.checkNoLongitudeBound()) {
+        leftLon = -Math.PI;
+        rightLon = Math.PI;
+      } else {
+        leftLon = bounds.getLeftLongitude().doubleValue();
+        rightLon = bounds.getRightLongitude().doubleValue();
+      }
+      double minLat;
+      if (bounds.checkNoBottomLatitudeBound()) {
+        minLat = -Math.PI * 0.5;
+      } else {
+        minLat = bounds.getMinLatitude().doubleValue();
+      }
+      double maxLat;
+      if (bounds.checkNoTopLatitudeBound()) {
+        maxLat = Math.PI * 0.5;
+      } else {
+        maxLat = bounds.getMaxLatitude().doubleValue();
+      }
+      return GeoBBoxFactory.makeGeoBBox(maxLat, minLat, leftLon, rightLon);
+  }
+
   @Test
   //@Seed("FAD1BAB12B6DCCFE")
   public void testGeoCircleRect() {

lucene/spatial/src/test/org/apache/lucene/spatial/spatial4j/geo3d/PlaneTest.java

@@ -0,0 +1,48 @@
+package org.apache.lucene.spatial.spatial4j.geo3d;
+
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+import org.junit.Test;
+
+import static org.junit.Assert.assertFalse;
+import static org.junit.Assert.assertTrue;
+
+/** Test basic plane functionality.
+*/
+public class PlaneTest {
+
+
+    @Test
+    public void testIdenticalPlanes() {
+        final GeoPoint p = new GeoPoint(0.123,-0.456);
+        final Plane plane1 = new Plane(p,0.0);
+        final Plane plane2 = new Plane(p,0.0);
+        assertTrue(plane1.isNumericallyIdentical(plane2));
+        final Plane plane3 = new Plane(p,0.1);
+        assertFalse(plane1.isNumericallyIdentical(plane3));
+        final Vector v1 = new Vector(0.1,-0.732,0.9);
+        final double constant = 0.432;
+        final Vector v2 = new Vector(v1.x*constant,v1.y*constant,v1.z*constant);
+        final Plane p1 = new Plane(v1,0.2);
+        final Plane p2 = new Plane(v2,0.2*constant);
+        assertTrue(p1.isNumericallyIdentical(p2));
+    }
+    
+
+}
+