lucene-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From nkn...@apache.org
Subject svn commit: r1718481 - in /lucene/dev/trunk/lucene: ./ core/src/java/org/apache/lucene/util/ sandbox/src/java/org/apache/lucene/search/ sandbox/src/java/org/apache/lucene/util/ sandbox/src/test/org/apache/lucene/search/ sandbox/src/test/org/apache/luce...
Date Mon, 07 Dec 2015 22:02:48 GMT
Author: nknize
Date: Mon Dec  7 22:02:48 2015
New Revision: 1718481

URL: http://svn.apache.org/viewvc?rev=1718481&view=rev
Log:
LUCENE-6908: Fix TestGeoUtils.testGeoRelations to handle irregular rectangles. Refactors relational methods to new GeoRelationUtils class.

Added:
    lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/util/GeoRelationUtils.java   (with props)
Modified:
    lucene/dev/trunk/lucene/CHANGES.txt
    lucene/dev/trunk/lucene/core/src/java/org/apache/lucene/util/SloppyMath.java
    lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/search/DimensionalPointInPolygonQuery.java
    lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/search/GeoPointDistanceQueryImpl.java
    lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/search/GeoPointInBBoxQueryImpl.java
    lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/search/GeoPointInPolygonQuery.java
    lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/search/GeoPointTermsEnum.java
    lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/util/GeoDistanceUtils.java
    lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/util/GeoProjectionUtils.java
    lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/util/GeoUtils.java
    lucene/dev/trunk/lucene/sandbox/src/test/org/apache/lucene/search/TestDimensionalQueries.java
    lucene/dev/trunk/lucene/sandbox/src/test/org/apache/lucene/search/TestGeoPointQuery.java
    lucene/dev/trunk/lucene/sandbox/src/test/org/apache/lucene/util/TestGeoUtils.java

Modified: lucene/dev/trunk/lucene/CHANGES.txt
URL: http://svn.apache.org/viewvc/lucene/dev/trunk/lucene/CHANGES.txt?rev=1718481&r1=1718480&r2=1718481&view=diff
==============================================================================
--- lucene/dev/trunk/lucene/CHANGES.txt (original)
+++ lucene/dev/trunk/lucene/CHANGES.txt Mon Dec  7 22:02:48 2015
@@ -116,6 +116,10 @@ New Features
 
 API Changes
 
+* LUCENE-6908: GeoUtils static relational methods have been refactored to new 
+  GeoRelationUtils and now correctly handle large irregular rectangles, and
+  pole crossing distance queries. (Nick Knize)
+
 * LUCENE-6900: Grouping sortWithinGroup variables used to allow null to mean
   Sort.RELEVANCE.  Null is no longer permitted.  (David Smiley)
 

Modified: lucene/dev/trunk/lucene/core/src/java/org/apache/lucene/util/SloppyMath.java
URL: http://svn.apache.org/viewvc/lucene/dev/trunk/lucene/core/src/java/org/apache/lucene/util/SloppyMath.java?rev=1718481&r1=1718480&r2=1718481&view=diff
==============================================================================
--- lucene/dev/trunk/lucene/core/src/java/org/apache/lucene/util/SloppyMath.java (original)
+++ lucene/dev/trunk/lucene/core/src/java/org/apache/lucene/util/SloppyMath.java Mon Dec  7 22:02:48 2015
@@ -87,7 +87,41 @@ public class SloppyMath {
     double indexSin = sinTab[index];
     return indexCos + delta * (-indexSin + delta * (-indexCos * ONE_DIV_F2 + delta * (indexSin * ONE_DIV_F3 + delta * indexCos * ONE_DIV_F4)));
   }
-  
+
+  /**
+   * Returns the trigonometric sine of an angle converted as a cos operation.
+   * <p>
+   * Error is around 1E-15.
+   * <p>
+   * Special cases:
+   * <ul>
+   *  <li>If the argument is {@code NaN} or an infinity, then the result is {@code NaN}.
+   * </ul>
+   * @param a an angle, in radians.
+   * @return the sine of the argument.
+   * @see Math#cos(double)
+   */
+  public static double sin(double a) {
+    return cos(a - PIO2);
+  }
+
+  /**
+   * Returns the trigonometric tangent of an angle converted in terms of a sin/cos operation.
+   * <p>
+   * Error is around 1E-15.
+   * <p>
+   * Special cases:
+   * <ul>
+   *  <li>If the argument is {@code NaN} or an infinity, then the result is {@code NaN}.
+   * </ul>
+   * @param a an angle, in radians.
+   * @return the tangent of the argument.
+   * @see Math#sin(double) aand Math#cos(double)
+   */
+  public static double tan(double a) {
+    return sin(a) / cos(a);
+  }
+
   /**
    * Returns the arc sine of a value.
    * <p>
@@ -146,13 +180,15 @@ public class SloppyMath {
   }
 
   // haversin
-  private static final double TO_RADIANS = Math.PI / 180D;
+  static final double TO_RADIANS = Math.PI / 180D;
+  static final double TO_DEGREES = 180D / Math.PI;
   
   // cos/asin
   private static final double ONE_DIV_F2 = 1/2.0;
   private static final double ONE_DIV_F3 = 1/6.0;
   private static final double ONE_DIV_F4 = 1/24.0;
-  
+
+  static final double PIO2 = Math.PI / 2D;
   private static final double PIO2_HI = Double.longBitsToDouble(0x3FF921FB54400000L); // 1.57079632673412561417e+00 first 33 bits of pi/2
   private static final double PIO2_LO = Double.longBitsToDouble(0x3DD0B4611A626331L); // 6.07710050650619224932e-11 pi/2 - PIO2_HI
   private static final double TWOPI_HI = 4*PIO2_HI;

Modified: lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/search/DimensionalPointInPolygonQuery.java
URL: http://svn.apache.org/viewvc/lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/search/DimensionalPointInPolygonQuery.java?rev=1718481&r1=1718480&r2=1718481&view=diff
==============================================================================
--- lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/search/DimensionalPointInPolygonQuery.java (original)
+++ lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/search/DimensionalPointInPolygonQuery.java Mon Dec  7 22:02:48 2015
@@ -27,6 +27,7 @@ import org.apache.lucene.index.Dimension
 import org.apache.lucene.index.LeafReader;
 import org.apache.lucene.index.LeafReaderContext;
 import org.apache.lucene.util.DocIdSetBuilder;
+import org.apache.lucene.util.GeoRelationUtils;
 import org.apache.lucene.util.GeoUtils;
 import org.apache.lucene.util.bkd.BKDUtil;
 
@@ -125,7 +126,7 @@ public class DimensionalPointInPolygonQu
                              assert packedValue.length == 8;
                              double lat = DimensionalLatLonField.decodeLat(BKDUtil.bytesToInt(packedValue, 0));
                              double lon = DimensionalLatLonField.decodeLon(BKDUtil.bytesToInt(packedValue, 1));
-                             if (GeoUtils.pointInPolygon(polyLons, polyLats, lat, lon)) {
+                             if (GeoRelationUtils.pointInPolygon(polyLons, polyLats, lat, lon)) {
                                hitCount[0]++;
                                result.add(docID);
                              }
@@ -141,11 +142,11 @@ public class DimensionalPointInPolygonQu
                              if (cellMinLat <= minLat && cellMaxLat >= maxLat && cellMinLon <= minLon && cellMaxLon >= maxLon) {
                                // Cell fully encloses the query
                                return Relation.CELL_CROSSES_QUERY;
-                             } else  if (GeoUtils.rectWithinPoly(cellMinLon, cellMinLat, cellMaxLon, cellMaxLat,
+                             } else  if (GeoRelationUtils.rectWithinPoly(cellMinLon, cellMinLat, cellMaxLon, cellMaxLat,
                                                                  polyLons, polyLats,
                                                                  minLon, minLat, maxLon, maxLat)) {
                                return Relation.CELL_INSIDE_QUERY;
-                             } else if (GeoUtils.rectCrossesPoly(cellMinLon, cellMinLat, cellMaxLon, cellMaxLat,
+                             } else if (GeoRelationUtils.rectCrossesPoly(cellMinLon, cellMinLat, cellMaxLon, cellMaxLat,
                                                                  polyLons, polyLats,
                                                                  minLon, minLat, maxLon, maxLat)) {
                                return Relation.CELL_CROSSES_QUERY;

Modified: lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/search/GeoPointDistanceQueryImpl.java
URL: http://svn.apache.org/viewvc/lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/search/GeoPointDistanceQueryImpl.java?rev=1718481&r1=1718480&r2=1718481&view=diff
==============================================================================
--- lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/search/GeoPointDistanceQueryImpl.java (original)
+++ lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/search/GeoPointDistanceQueryImpl.java Mon Dec  7 22:02:48 2015
@@ -24,6 +24,7 @@ import org.apache.lucene.index.Terms;
 import org.apache.lucene.index.TermsEnum;
 import org.apache.lucene.util.AttributeSource;
 import org.apache.lucene.util.GeoRect;
+import org.apache.lucene.util.GeoRelationUtils;
 import org.apache.lucene.util.GeoUtils;
 import org.apache.lucene.util.SloppyMath;
 
@@ -77,12 +78,14 @@ final class GeoPointDistanceQueryImpl ex
 
     @Override
     protected boolean cellCrosses(final double minLon, final double minLat, final double maxLon, final double maxLat) {
-      return GeoUtils.rectCrossesCircle(minLon, minLat, maxLon, maxLat, centerLon, query.centerLat, query.radiusMeters);
+      return GeoRelationUtils.rectCrossesCircle(minLon, minLat, maxLon, maxLat,
+          centerLon, query.centerLat, query.radiusMeters, true);
     }
 
     @Override
     protected boolean cellWithin(final double minLon, final double minLat, final double maxLon, final double maxLat) {
-      return GeoUtils.rectWithinCircle(minLon, minLat, maxLon, maxLat, centerLon, query.centerLat, query.radiusMeters);
+      return GeoRelationUtils.rectWithinCircle(minLon, minLat, maxLon, maxLat,
+          centerLon, query.centerLat, query.radiusMeters, true);
     }
 
     @Override

Modified: lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/search/GeoPointInBBoxQueryImpl.java
URL: http://svn.apache.org/viewvc/lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/search/GeoPointInBBoxQueryImpl.java?rev=1718481&r1=1718480&r2=1718481&view=diff
==============================================================================
--- lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/search/GeoPointInBBoxQueryImpl.java (original)
+++ lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/search/GeoPointInBBoxQueryImpl.java Mon Dec  7 22:02:48 2015
@@ -23,6 +23,7 @@ import org.apache.lucene.document.GeoPoi
 import org.apache.lucene.index.Terms;
 import org.apache.lucene.index.TermsEnum;
 import org.apache.lucene.util.AttributeSource;
+import org.apache.lucene.util.GeoRelationUtils;
 import org.apache.lucene.util.GeoUtils;
 import org.apache.lucene.util.SloppyMath;
 import org.apache.lucene.util.ToStringUtils;
@@ -83,7 +84,7 @@ class GeoPointInBBoxQueryImpl extends Ge
      */
     @Override
     protected boolean cellCrosses(final double minLon, final double minLat, final double maxLon, final double maxLat) {
-      return GeoUtils.rectCrosses(minLon, minLat, maxLon, maxLat, this.minLon, this.minLat, this.maxLon, this.maxLat);
+      return GeoRelationUtils.rectCrosses(minLon, minLat, maxLon, maxLat, this.minLon, this.minLat, this.maxLon, this.maxLat);
     }
 
     /**
@@ -91,7 +92,7 @@ class GeoPointInBBoxQueryImpl extends Ge
      */
     @Override
     protected boolean cellWithin(final double minLon, final double minLat, final double maxLon, final double maxLat) {
-      return GeoUtils.rectWithin(minLon, minLat, maxLon, maxLat, this.minLon, this.minLat, this.maxLon, this.maxLat);
+      return GeoRelationUtils.rectWithin(minLon, minLat, maxLon, maxLat, this.minLon, this.minLat, this.maxLon, this.maxLat);
     }
 
     @Override
@@ -101,7 +102,7 @@ class GeoPointInBBoxQueryImpl extends Ge
 
     @Override
     protected boolean postFilter(final double lon, final double lat) {
-      return GeoUtils.bboxContains(lon, lat, minLon, minLat, maxLon, maxLat);
+      return GeoRelationUtils.pointInRect(lon, lat, minLon, minLat, maxLon, maxLat);
     }
   }
 

Modified: lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/search/GeoPointInPolygonQuery.java
URL: http://svn.apache.org/viewvc/lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/search/GeoPointInPolygonQuery.java?rev=1718481&r1=1718480&r2=1718481&view=diff
==============================================================================
--- lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/search/GeoPointInPolygonQuery.java (original)
+++ lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/search/GeoPointInPolygonQuery.java Mon Dec  7 22:02:48 2015
@@ -24,6 +24,7 @@ import org.apache.lucene.index.Terms;
 import org.apache.lucene.index.TermsEnum;
 import org.apache.lucene.util.AttributeSource;
 import org.apache.lucene.util.GeoRect;
+import org.apache.lucene.util.GeoRelationUtils;
 import org.apache.lucene.util.GeoUtils;
 
 /** Implements a simple point in polygon query on a GeoPoint field. This is based on
@@ -147,13 +148,13 @@ public final class GeoPointInPolygonQuer
 
     @Override
     protected boolean cellCrosses(final double minLon, final double minLat, final double maxLon, final double maxLat) {
-      return GeoUtils.rectCrossesPoly(minLon, minLat, maxLon, maxLat, x, y, GeoPointInPolygonQuery.this.minLon,
+      return GeoRelationUtils.rectCrossesPoly(minLon, minLat, maxLon, maxLat, x, y, GeoPointInPolygonQuery.this.minLon,
           GeoPointInPolygonQuery.this.minLat, GeoPointInPolygonQuery.this.maxLon, GeoPointInPolygonQuery.this.maxLat);
     }
 
     @Override
     protected boolean cellWithin(final double minLon, final double minLat, final double maxLon, final double maxLat) {
-      return GeoUtils.rectWithinPoly(minLon, minLat, maxLon, maxLat, x, y, GeoPointInPolygonQuery.this.minLon,
+      return GeoRelationUtils.rectWithinPoly(minLon, minLat, maxLon, maxLat, x, y, GeoPointInPolygonQuery.this.minLon,
           GeoPointInPolygonQuery.this.minLat, GeoPointInPolygonQuery.this.maxLon, GeoPointInPolygonQuery.this.maxLat);
     }
 
@@ -168,11 +169,11 @@ public final class GeoPointInPolygonQuer
      * {@link org.apache.lucene.search.GeoPointTermsEnum#accept} method is called to match
      * encoded terms that fall within the bounding box of the polygon. Those documents that pass the initial
      * bounding box filter are then compared to the provided polygon using the
-     * {@link org.apache.lucene.util.GeoUtils#pointInPolygon} method.
+     * {@link org.apache.lucene.util.GeoRelationUtils#pointInPolygon} method.
      */
     @Override
     protected boolean postFilter(final double lon, final double lat) {
-      return GeoUtils.pointInPolygon(x, y, lat, lon);
+      return GeoRelationUtils.pointInPolygon(x, y, lat, lon);
     }
   }
 

Modified: lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/search/GeoPointTermsEnum.java
URL: http://svn.apache.org/viewvc/lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/search/GeoPointTermsEnum.java?rev=1718481&r1=1718480&r2=1718481&view=diff
==============================================================================
--- lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/search/GeoPointTermsEnum.java (original)
+++ lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/search/GeoPointTermsEnum.java Mon Dec  7 22:02:48 2015
@@ -26,6 +26,7 @@ import org.apache.lucene.index.FilteredT
 import org.apache.lucene.index.TermsEnum;
 import org.apache.lucene.util.BytesRef;
 import org.apache.lucene.util.BytesRefBuilder;
+import org.apache.lucene.util.GeoRelationUtils;
 import org.apache.lucene.util.GeoUtils;
 import org.apache.lucene.util.NumericUtils;
 
@@ -139,14 +140,14 @@ abstract class GeoPointTermsEnum extends
    * Primary driver for cells intersecting shape boundaries
    */
   protected boolean cellIntersectsMBR(final double minLon, final double minLat, final double maxLon, final double maxLat) {
-    return GeoUtils.rectIntersects(minLon, minLat, maxLon, maxLat, this.minLon, this.minLat, this.maxLon, this.maxLat);
+    return GeoRelationUtils.rectIntersects(minLon, minLat, maxLon, maxLat, this.minLon, this.minLat, this.maxLon, this.maxLat);
   }
 
   /**
    * Return whether quad-cell contains the bounding box of this shape
    */
   protected boolean cellContains(final double minLon, final double minLat, final double maxLon, final double maxLat) {
-    return GeoUtils.rectWithin(this.minLon, this.minLat, this.maxLon, this.maxLat, minLon, minLat, maxLon, maxLat);
+    return GeoRelationUtils.rectWithin(this.minLon, this.minLat, this.maxLon, this.maxLat, minLon, minLat, maxLon, maxLat);
   }
 
   public boolean boundaryTerm() {

Modified: lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/util/GeoDistanceUtils.java
URL: http://svn.apache.org/viewvc/lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/util/GeoDistanceUtils.java?rev=1718481&r1=1718480&r2=1718481&view=diff
==============================================================================
--- lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/util/GeoDistanceUtils.java (original)
+++ lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/util/GeoDistanceUtils.java Mon Dec  7 22:02:48 2015
@@ -17,6 +17,8 @@ package org.apache.lucene.util;
  * limitations under the License.
  */
 
+import static org.apache.lucene.util.SloppyMath.TO_RADIANS;
+
 /**
  * Reusable geo-spatial distance utility methods.
  *
@@ -25,6 +27,28 @@ package org.apache.lucene.util;
 public class GeoDistanceUtils {
 
   /**
+   * Compute the great-circle distance using original haversine implementation published by Sinnot in:
+   * R.W. Sinnott, "Virtues of the Haversine", Sky and Telescope, vol. 68, no. 2, 1984, p. 159
+   *
+   * NOTE: this differs from {@link org.apache.lucene.util.SloppyMath#haversin} in that it uses the semi-major axis
+   * of the earth instead of an approximation based on the average latitude of the two points (which can introduce an
+   * additional error up to .337%, or ~67.6 km at the equator)
+   */
+  public static double haversin(double lat1, double lon1, double lat2, double lon2) {
+    double dLat = TO_RADIANS * (lat2 - lat1);
+    double dLon = TO_RADIANS * (lon2 - lon1);
+    lat1 = TO_RADIANS * (lat1);
+    lat2 = TO_RADIANS * (lat2);
+
+    final double sinDLatO2 = SloppyMath.sin(dLat / 2);
+    final double sinDLonO2 = SloppyMath.sin(dLon / 2);
+
+    double a = sinDLatO2*sinDLatO2 + sinDLonO2 * sinDLonO2 * SloppyMath.cos(lat1) * SloppyMath.cos(lat2);
+    double c = 2 * SloppyMath.asin(Math.sqrt(a));
+    return (GeoProjectionUtils.SEMIMAJOR_AXIS * c);
+  }
+
+  /**
    * Compute the distance between two geo-points using vincenty distance formula
    * Vincenty uses the oblate spheroid whereas haversine uses unit sphere, this will give roughly
    * 22m better accuracy (in worst case) than haversine
@@ -87,6 +111,20 @@ public class GeoDistanceUtils {
   }
 
   /**
+   * Computes distance between two points in a cartesian (x, y, {z - optional}) coordinate system
+   */
+  public static double linearDistance(double[] pt1, double[] pt2) {
+    assert pt1 != null && pt2 != null && pt1.length == pt2.length && pt1.length > 1;
+    final double d0 = pt1[0] - pt2[0];
+    final double d1 = pt1[1] - pt2[1];
+    if (pt1.length == 3) {
+      final double d2 = pt1[2] - pt2[2];
+      return Math.sqrt(d0*d0 + d1*d1 + d2*d2);
+    }
+    return Math.sqrt(d0*d0 + d1*d1);
+  }
+
+  /**
    * Compute the inverse haversine to determine distance in degrees longitude for provided distance in meters
    * @param lat latitude to compute delta degrees lon
    * @param distance distance in meters to convert to degrees lon
@@ -152,9 +190,9 @@ public class GeoDistanceUtils {
   /** Returns the maximum distance/radius (in meters) from the point 'center' before overlapping */
   public static double maxRadialDistanceMeters(final double centerLon, final double centerLat) {
     if (Math.abs(centerLat) == GeoUtils.MAX_LAT_INCL) {
-      return SloppyMath.haversin(centerLat, centerLon, 0, centerLon)*1000.0;
+      return GeoDistanceUtils.haversin(centerLat, centerLon, 0, centerLon);
     }
-    return SloppyMath.haversin(centerLat, centerLon, centerLat, (GeoUtils.MAX_LON_INCL + centerLon) % 360)*1000.0;
+    return GeoDistanceUtils.haversin(centerLat, centerLon, centerLat, (GeoUtils.MAX_LON_INCL + centerLon) % 360);
   }
 
   /**

Modified: lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/util/GeoProjectionUtils.java
URL: http://svn.apache.org/viewvc/lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/util/GeoProjectionUtils.java?rev=1718481&r1=1718480&r2=1718481&view=diff
==============================================================================
--- lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/util/GeoProjectionUtils.java (original)
+++ lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/util/GeoProjectionUtils.java Mon Dec  7 22:02:48 2015
@@ -17,6 +17,10 @@ package org.apache.lucene.util;
  * limitations under the License.
  */
 
+import static org.apache.lucene.util.SloppyMath.PIO2;
+import static org.apache.lucene.util.SloppyMath.TO_DEGREES;
+import static org.apache.lucene.util.SloppyMath.TO_RADIANS;
+
 /**
  * Reusable geo-spatial projection utility methods.
  *
@@ -28,13 +32,14 @@ public class GeoProjectionUtils {
   public static final double FLATTENING = 1.0/298.257223563;
   public static final double SEMIMINOR_AXIS = SEMIMAJOR_AXIS * (1.0 - FLATTENING); //6_356_752.31420; // [m]
   public static final double ECCENTRICITY = StrictMath.sqrt((2.0 - FLATTENING) * FLATTENING);
-  static final double PI_OVER_2 = StrictMath.PI / 2.0D;
   static final double SEMIMAJOR_AXIS2 = SEMIMAJOR_AXIS * SEMIMAJOR_AXIS;
   static final double SEMIMINOR_AXIS2 = SEMIMINOR_AXIS * SEMIMINOR_AXIS;
-  public static final double MIN_LON_RADIANS = StrictMath.toRadians(GeoUtils.MIN_LON_INCL);
-  public static final double MIN_LAT_RADIANS = StrictMath.toRadians(GeoUtils.MIN_LAT_INCL);
-  public static final double MAX_LON_RADIANS = StrictMath.toRadians(GeoUtils.MAX_LON_INCL);
-  public static final double MAX_LAT_RADIANS = StrictMath.toRadians(GeoUtils.MAX_LAT_INCL);
+  public static final double MIN_LON_RADIANS = TO_RADIANS * GeoUtils.MIN_LON_INCL;
+  public static final double MIN_LAT_RADIANS = TO_RADIANS * GeoUtils.MIN_LAT_INCL;
+  public static final double MAX_LON_RADIANS = TO_RADIANS * GeoUtils.MAX_LON_INCL;
+  public static final double MAX_LAT_RADIANS = TO_RADIANS * GeoUtils.MAX_LAT_INCL;
+  private static final double E2 = (SEMIMAJOR_AXIS2 - SEMIMINOR_AXIS2)/(SEMIMAJOR_AXIS2);
+  private static final double EP2 = (SEMIMAJOR_AXIS2 - SEMIMINOR_AXIS2)/(SEMIMINOR_AXIS2);
 
   /**
    * Converts from geocentric earth-centered earth-fixed to geodesic lat/lon/alt
@@ -47,8 +52,6 @@ public class GeoProjectionUtils {
   public static final double[] ecfToLLA(final double x, final double y, final double z, double[] lla) {
     boolean atPole = false;
     final double ad_c = 1.0026000D;
-    final double e2 = (SEMIMAJOR_AXIS2 - SEMIMINOR_AXIS2)/(SEMIMAJOR_AXIS2);
-    final double ep2 = (SEMIMAJOR_AXIS2 - SEMIMINOR_AXIS2)/(SEMIMINOR_AXIS2);
     final double cos67P5 = 0.38268343236508977D;
 
     if (lla == null) {
@@ -59,18 +62,18 @@ public class GeoProjectionUtils {
       lla[0] = StrictMath.atan2(y,x);
     } else {
       if (y > 0) {
-        lla[0] = PI_OVER_2;
+        lla[0] = PIO2;
       } else if (y < 0) {
-        lla[0] = -PI_OVER_2;
+        lla[0] = -PIO2;
       } else {
         atPole = true;
         lla[0] = 0.0D;
         if (z > 0.0) {
-          lla[1] = PI_OVER_2;
+          lla[1] = PIO2;
         } else if (z < 0.0) {
-          lla[1] = -PI_OVER_2;
+          lla[1] = -PIO2;
         } else {
-          lla[1] = PI_OVER_2;
+          lla[1] = PIO2;
           lla[2] = -SEMIMINOR_AXIS;
           return lla;
         }
@@ -84,25 +87,25 @@ public class GeoProjectionUtils {
     final double sinB0 = t0 / s0;
     final double cosB0 = w / s0;
     final double sin3B0 = sinB0 * sinB0 * sinB0;
-    final double t1 = z + SEMIMINOR_AXIS * ep2 * sin3B0;
-    final double sum = w - SEMIMAJOR_AXIS * e2 * cosB0 * cosB0 * cosB0;
+    final double t1 = z + SEMIMINOR_AXIS * EP2 * sin3B0;
+    final double sum = w - SEMIMAJOR_AXIS * E2 * cosB0 * cosB0 * cosB0;
     final double s1 = StrictMath.sqrt(t1 * t1 + sum * sum);
     final double sinP1 = t1 / s1;
     final double cosP1 = sum / s1;
-    final double rn = SEMIMAJOR_AXIS / StrictMath.sqrt(1.0D - e2 * sinP1 * sinP1);
+    final double rn = SEMIMAJOR_AXIS / StrictMath.sqrt(1.0D - E2 * sinP1 * sinP1);
 
     if (cosP1 >= cos67P5) {
       lla[2] = w / cosP1 - rn;
     } else if (cosP1 <= -cos67P5) {
       lla[2] = w / -cosP1 - rn;
     } else {
-      lla[2] = z / sinP1 + rn * (e2 - 1.0);
+      lla[2] = z / sinP1 + rn * (E2 - 1.0);
     }
     if (!atPole) {
       lla[1] = StrictMath.atan(sinP1/cosP1);
     }
-    lla[0] = StrictMath.toDegrees(lla[0]);
-    lla[1] = StrictMath.toDegrees(lla[1]);
+    lla[0] = TO_DEGREES * lla[0];
+    lla[1] = TO_DEGREES * lla[1];
 
     return lla;
   }
@@ -116,33 +119,32 @@ public class GeoProjectionUtils {
    * @return either a new ecef array or the reusable ecf parameter
    */
   public static final double[] llaToECF(double lon, double lat, double alt, double[] ecf) {
-    lon = StrictMath.toRadians(lon);
-    lat = StrictMath.toRadians(lat);
+    lon = TO_RADIANS * lon;
+    lat = TO_RADIANS * lat;
 
-    final double sl = StrictMath.sin(lat);
+    final double sl = SloppyMath.sin(lat);
     final double s2 = sl*sl;
-    final double cl = StrictMath.cos(lat);
-    final double ge2 = (SEMIMAJOR_AXIS2 - SEMIMINOR_AXIS2)/(SEMIMAJOR_AXIS2);
+    final double cl = SloppyMath.cos(lat);
 
     if (ecf == null) {
       ecf = new double[3];
     }
 
-    if (lat < -PI_OVER_2 && lat > -1.001D * PI_OVER_2) {
-      lat = -PI_OVER_2;
-    } else if (lat > PI_OVER_2 && lat < 1.001D * PI_OVER_2) {
-      lat = PI_OVER_2;
+    if (lat < -PIO2 && lat > -1.001D * PIO2) {
+      lat = -PIO2;
+    } else if (lat > PIO2 && lat < 1.001D * PIO2) {
+      lat = PIO2;
     }
-    assert (lat >= -PI_OVER_2) || (lat <= PI_OVER_2);
+    assert (lat >= -PIO2) || (lat <= PIO2);
 
     if (lon > StrictMath.PI) {
       lon -= (2*StrictMath.PI);
     }
 
-    final double rn = SEMIMAJOR_AXIS / StrictMath.sqrt(1.0D - ge2 * s2);
-    ecf[0] = (rn+alt) * cl * StrictMath.cos(lon);
-    ecf[1] = (rn+alt) * cl * StrictMath.sin(lon);
-    ecf[2] = ((rn*(1.0-ge2))+alt)*sl;
+    final double rn = SEMIMAJOR_AXIS / StrictMath.sqrt(1.0D - E2 * s2);
+    ecf[0] = (rn+alt) * cl * SloppyMath.cos(lon);
+    ecf[1] = (rn+alt) * cl * SloppyMath.sin(lon);
+    ecf[2] = ((rn*(1.0-E2))+alt)*sl;
 
     return ecf;
   }
@@ -272,13 +274,13 @@ public class GeoProjectionUtils {
       phiMatrix = new double[3][3];
     }
 
-    originLon = StrictMath.toRadians(originLon);
-    originLat = StrictMath.toRadians(originLat);
+    originLon = TO_RADIANS * originLon;
+    originLat = TO_RADIANS * originLat;
 
-    final double sLon = StrictMath.sin(originLon);
-    final double cLon = StrictMath.cos(originLon);
-    final double sLat = StrictMath.sin(originLat);
-    final double cLat = StrictMath.cos(originLat);
+    final double sLon = SloppyMath.sin(originLon);
+    final double cLon = SloppyMath.cos(originLon);
+    final double sLat = SloppyMath.sin(originLat);
+    final double cLat = SloppyMath.cos(originLat);
 
     phiMatrix[0][0] = -sLon;
     phiMatrix[0][1] = cLon;
@@ -306,13 +308,13 @@ public class GeoProjectionUtils {
       phiMatrix = new double[3][3];
     }
 
-    originLon = StrictMath.toRadians(originLon);
-    originLat = StrictMath.toRadians(originLat);
+    originLon = TO_RADIANS * originLon;
+    originLat = TO_RADIANS * originLat;
 
-    final double sLat = StrictMath.sin(originLat);
-    final double cLat = StrictMath.cos(originLat);
-    final double sLon = StrictMath.sin(originLon);
-    final double cLon = StrictMath.cos(originLon);
+    final double sLat = SloppyMath.sin(originLat);
+    final double cLat = SloppyMath.cos(originLat);
+    final double sLon = SloppyMath.sin(originLon);
+    final double cLon = SloppyMath.cos(originLon);
 
     phiMatrix[0][0] = -sLon;
     phiMatrix[1][0] = cLon;
@@ -337,22 +339,22 @@ public class GeoProjectionUtils {
    * @param pt resulting point
    * @return the point along a bearing at a given distance in meters
    */
-  public static final double[] pointFromLonLatBearing(double lon, double lat, double bearing, double dist, double[] pt) {
+  public static final double[] pointFromLonLatBearingVincenty(double lon, double lat, double bearing, double dist, double[] pt) {
 
     if (pt == null) {
       pt = new double[2];
     }
 
-    final double alpha1 = StrictMath.toRadians(bearing);
-    final double cosA1 = StrictMath.cos(alpha1);
-    final double sinA1 = StrictMath.sin(alpha1);
-    final double tanU1 = (1-FLATTENING) * StrictMath.tan(StrictMath.toRadians(lat));
+    final double alpha1 = TO_RADIANS * bearing;
+    final double cosA1 = SloppyMath.cos(alpha1);
+    final double sinA1 = SloppyMath.sin(alpha1);
+    final double tanU1 = (1-FLATTENING) * SloppyMath.tan(TO_RADIANS * lat);
     final double cosU1 = 1 / StrictMath.sqrt((1+tanU1*tanU1));
     final double sinU1 = tanU1*cosU1;
     final double sig1 = StrictMath.atan2(tanU1, cosA1);
     final double sinAlpha = cosU1 * sinA1;
     final double cosSqAlpha = 1 - sinAlpha*sinAlpha;
-    final double uSq = cosSqAlpha * (SEMIMAJOR_AXIS2 - SEMIMINOR_AXIS2) / SEMIMINOR_AXIS2;
+    final double uSq = cosSqAlpha * EP2;
     final double A = 1 + uSq/16384D*(4096D + uSq * (-768D + uSq * (320D - 175D*uSq)));
     final double B = uSq/1024D * (256D + uSq * (-128D + uSq * (74D - 47D * uSq)));
 
@@ -361,9 +363,9 @@ public class GeoProjectionUtils {
     double sinSigma, cosSigma, cos2SigmaM, deltaSigma;
 
     do {
-      cos2SigmaM = StrictMath.cos(2*sig1 + sigma);
-      sinSigma = StrictMath.sin(sigma);
-      cosSigma = StrictMath.cos(sigma);
+      cos2SigmaM = SloppyMath.cos(2*sig1 + sigma);
+      sinSigma = SloppyMath.sin(sigma);
+      cosSigma = SloppyMath.cos(sigma);
 
       deltaSigma = B * sinSigma * (cos2SigmaM + (B/4D) * (cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-
           (B/6) * cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)));
@@ -379,9 +381,42 @@ public class GeoProjectionUtils {
 
     final double lam = lambda - (1-c) * FLATTENING * sinAlpha *
         (sigma + c * sinSigma * (cos2SigmaM + c * cosSigma * (-1 + 2* cos2SigmaM*cos2SigmaM)));
-    pt[0] = GeoUtils.normalizeLon(lon + StrictMath.toDegrees(lam));
-    pt[1] = GeoUtils.normalizeLat(StrictMath.toDegrees(lat2));
+    pt[0] = GeoUtils.normalizeLon(lon + TO_DEGREES * lam);
+    pt[1] = GeoUtils.normalizeLat(TO_DEGREES * lat2);
+
+    return pt;
+  }
+
+  /**
+   * Finds a point along a bearing from a given lon,lat geolocation using great circle arc
+   *
+   * @param lon origin longitude in degrees
+   * @param lat origin latitude in degrees
+   * @param bearing azimuthal bearing in degrees
+   * @param dist distance in meters
+   * @param pt resulting point
+   * @return the point along a bearing at a given distance in meters
+   */
+  public static final double[] pointFromLonLatBearingGreatCircle(double lon, double lat, double bearing, double dist, double[] pt) {
+
+    if (pt == null) {
+      pt = new double[2];
+    }
+
+    lon *= TO_RADIANS;
+    lat *= TO_RADIANS;
+    bearing *= TO_RADIANS;
+
+    final double cLat = SloppyMath.cos(lat);
+    final double sLat = SloppyMath.sin(lat);
+    final double sinDoR = SloppyMath.sin(dist / GeoProjectionUtils.SEMIMAJOR_AXIS);
+    final double cosDoR = SloppyMath.cos(dist / GeoProjectionUtils.SEMIMAJOR_AXIS);
+
+    pt[1] = SloppyMath.asin(sLat*cosDoR + cLat * sinDoR * SloppyMath.cos(bearing));
+    pt[0] = TO_DEGREES * (lon + Math.atan2(SloppyMath.sin(bearing) * sinDoR * cLat, cosDoR - sLat * SloppyMath.sin(pt[1])));
+    pt[1] *= TO_DEGREES;
 
     return pt;
   }
+
 }

Added: lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/util/GeoRelationUtils.java
URL: http://svn.apache.org/viewvc/lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/util/GeoRelationUtils.java?rev=1718481&view=auto
==============================================================================
--- lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/util/GeoRelationUtils.java (added)
+++ lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/util/GeoRelationUtils.java Mon Dec  7 22:02:48 2015
@@ -0,0 +1,356 @@
+package org.apache.lucene.util;
+
+/*
+ * 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.
+ */
+
+/**
+ * Reusable geo-relation utility methods
+ */
+public class GeoRelationUtils {
+
+  /**
+   * Determine if a bbox (defined by minLon, minLat, maxLon, maxLat) contains the provided point (defined by lon, lat)
+   * NOTE: this is a basic method that does not handle dateline or pole crossing. Unwrapping must be done before
+   * calling this method.
+   */
+  public static boolean pointInRect(final double lon, final double lat, final double minLon,
+                                    final double minLat, final double maxLon, final double maxLat) {
+    return (GeoUtils.compare(lon, minLon) >= 0 && GeoUtils.compare(lon, maxLon) <= 0
+        && GeoUtils.compare(lat, minLat) >= 0 && GeoUtils.compare(lat, maxLat) <= 0);
+  }
+
+  /**
+   * simple even-odd point in polygon computation
+   *    1.  Determine if point is contained in the longitudinal range
+   *    2.  Determine whether point crosses the edge by computing the latitudinal delta
+   *        between the end-point of a parallel vector (originating at the point) and the
+   *        y-component of the edge sink
+   *
+   * NOTE: Requires polygon point (x,y) order either clockwise or counter-clockwise
+   */
+  public static boolean pointInPolygon(double[] x, double[] y, double lat, double lon) {
+    assert x.length == y.length;
+    boolean inPoly = false;
+    /**
+     * Note: This is using a euclidean coordinate system which could result in
+     * upwards of 110KM error at the equator.
+     * TODO convert coordinates to cylindrical projection (e.g. mercator)
+     */
+    for (int i = 1; i < x.length; i++) {
+      if (x[i] < lon && x[i-1] >= lon || x[i-1] < lon && x[i] >= lon) {
+        if (y[i] + (lon - x[i]) / (x[i-1] - x[i]) * (y[i-1] - y[i]) < lat) {
+          inPoly = !inPoly;
+        }
+      }
+    }
+    return inPoly;
+  }
+
+  public static boolean rectDisjoint(final double aMinX, final double aMinY, final double aMaxX, final double aMaxY,
+                                     final double bMinX, final double bMinY, final double bMaxX, final double bMaxY) {
+    return (aMaxX < bMinX || aMinX > bMaxX || aMaxY < bMinY || aMinY > bMaxY);
+  }
+
+  /**
+   * Computes whether the first (a) rectangle is wholly within another (b) rectangle (shared boundaries allowed)
+   */
+  public static boolean rectWithin(final double aMinX, final double aMinY, final double aMaxX, final double aMaxY,
+                                   final double bMinX, final double bMinY, final double bMaxX, final double bMaxY) {
+    return !(aMinX < bMinX || aMinY < bMinY || aMaxX > bMaxX || aMaxY > bMaxY);
+  }
+
+  public static boolean rectCrosses(final double aMinX, final double aMinY, final double aMaxX, final double aMaxY,
+                                    final double bMinX, final double bMinY, final double bMaxX, final double bMaxY) {
+    return !(rectDisjoint(aMinX, aMinY, aMaxX, aMaxY, bMinX, bMinY, bMaxX, bMaxY) ||
+        rectWithin(aMinX, aMinY, aMaxX, aMaxY, bMinX, bMinY, bMaxX, bMaxY));
+  }
+
+  /**
+   * Computes whether rectangle a contains rectangle b (touching allowed)
+   */
+  public static boolean rectContains(final double aMinX, final double aMinY, final double aMaxX, final double aMaxY,
+                                     final double bMinX, final double bMinY, final double bMaxX, final double bMaxY) {
+    return !(bMinX < aMinX || bMinY < aMinY || bMaxX > aMaxX || bMaxY > aMaxY);
+  }
+
+  /**
+   * Computes whether a rectangle intersects another rectangle (crosses, within, touching, etc)
+   */
+  public static boolean rectIntersects(final double aMinX, final double aMinY, final double aMaxX, final double aMaxY,
+                                       final double bMinX, final double bMinY, final double bMaxX, final double bMaxY) {
+    return !((aMaxX < bMinX || aMinX > bMaxX || aMaxY < bMinY || aMinY > bMaxY) );
+  }
+
+  /**
+   * Computes whether a rectangle crosses a shape. (touching not allowed)
+   */
+  public static boolean rectCrossesPoly(final double rMinX, final double rMinY, final double rMaxX,
+                                        final double rMaxY, final double[] shapeX, final double[] shapeY,
+                                        final double sMinX, final double sMinY, final double sMaxX,
+                                        final double sMaxY) {
+    // short-circuit: if the bounding boxes are disjoint then the shape does not cross
+    if (rectDisjoint(rMinX, rMinY, rMaxX, rMaxY, sMinX, sMinY, sMaxX, sMaxY)) {
+      return false;
+    }
+
+    final double[][] bbox = new double[][] { {rMinX, rMinY}, {rMaxX, rMinY}, {rMaxX, rMaxY}, {rMinX, rMaxY}, {rMinX, rMinY} };
+    final int polyLength = shapeX.length-1;
+    double d, s, t, a1, b1, c1, a2, b2, c2;
+    double x00, y00, x01, y01, x10, y10, x11, y11;
+
+    // computes the intersection point between each bbox edge and the polygon edge
+    for (short b=0; b<4; ++b) {
+      a1 = bbox[b+1][1]-bbox[b][1];
+      b1 = bbox[b][0]-bbox[b+1][0];
+      c1 = a1*bbox[b+1][0] + b1*bbox[b+1][1];
+      for (int p=0; p<polyLength; ++p) {
+        a2 = shapeY[p+1]-shapeY[p];
+        b2 = shapeX[p]-shapeX[p+1];
+        // compute determinant
+        d = a1*b2 - a2*b1;
+        if (d != 0) {
+          // lines are not parallel, check intersecting points
+          c2 = a2*shapeX[p+1] + b2*shapeY[p+1];
+          s = (1/d)*(b2*c1 - b1*c2);
+          t = (1/d)*(a1*c2 - a2*c1);
+          x00 = StrictMath.min(bbox[b][0], bbox[b+1][0]) - GeoUtils.TOLERANCE;
+          x01 = StrictMath.max(bbox[b][0], bbox[b+1][0]) + GeoUtils.TOLERANCE;
+          y00 = StrictMath.min(bbox[b][1], bbox[b+1][1]) - GeoUtils.TOLERANCE;
+          y01 = StrictMath.max(bbox[b][1], bbox[b+1][1]) + GeoUtils.TOLERANCE;
+          x10 = StrictMath.min(shapeX[p], shapeX[p+1]) - GeoUtils.TOLERANCE;
+          x11 = StrictMath.max(shapeX[p], shapeX[p+1]) + GeoUtils.TOLERANCE;
+          y10 = StrictMath.min(shapeY[p], shapeY[p+1]) - GeoUtils.TOLERANCE;
+          y11 = StrictMath.max(shapeY[p], shapeY[p+1]) + GeoUtils.TOLERANCE;
+          // check whether the intersection point is touching one of the line segments
+          boolean touching = ((x00 == s && y00 == t) || (x01 == s && y01 == t))
+              || ((x10 == s && y10 == t) || (x11 == s && y11 == t));
+          // if line segments are not touching and the intersection point is within the range of either segment
+          if (!(touching || x00 > s || x01 < s || y00 > t || y01 < t || x10 > s || x11 < s || y10 > t || y11 < t)) {
+            return true;
+          }
+        }
+      } // for each poly edge
+    } // for each bbox edge
+    return false;
+  }
+
+  /**
+   * Computes whether a rectangle is within a given polygon (shared boundaries allowed)
+   */
+  public static boolean rectWithinPoly(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
+                                       final double[] shapeX, final double[] shapeY, final double sMinX,
+                                       final double sMinY, final double sMaxX, final double sMaxY) {
+    // check if rectangle crosses poly (to handle concave/pacman polys), then check that all 4 corners
+    // are contained
+    return !(rectCrossesPoly(rMinX, rMinY, rMaxX, rMaxY, shapeX, shapeY, sMinX, sMinY, sMaxX, sMaxY) ||
+        !pointInPolygon(shapeX, shapeY, rMinY, rMinX) || !pointInPolygon(shapeX, shapeY, rMinY, rMaxX) ||
+        !pointInPolygon(shapeX, shapeY, rMaxY, rMaxX) || !pointInPolygon(shapeX, shapeY, rMaxY, rMinX));
+  }
+
+  private static boolean rectAnyCornersInCircle(final double rMinX, final double rMinY, final double rMaxX,
+                                                final double rMaxY, final double centerLon, final double centerLat,
+                                                final double radiusMeters, final boolean approx) {
+    if (approx == true) {
+      return rectAnyCornersInCircleSloppy(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters);
+    }
+    double w = Math.abs(rMaxX - rMinX);
+    if (w <= 90.0) {
+      return GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, rMinX) <= radiusMeters
+          || GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, rMinX) <= radiusMeters
+          || GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, rMaxX) <= radiusMeters
+          || GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, rMaxX) <= radiusMeters;
+    }
+    // partition
+    w /= 4;
+    final double p1 = rMinX + w;
+    final double p2 = p1 + w;
+    final double p3 = p2 + w;
+
+    return GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, rMinX) <= radiusMeters
+        || GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, rMinX) <= radiusMeters
+        || GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, p1) <= radiusMeters
+        || GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, p1) <= radiusMeters
+        || GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, p2) <= radiusMeters
+        || GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, p2) <= radiusMeters
+        || GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, p3) <= radiusMeters
+        || GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, p3) <= radiusMeters
+        || GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, rMaxX) <= radiusMeters
+        || GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, rMaxX) <= radiusMeters;
+  }
+
+  private static boolean rectAnyCornersInCircleSloppy(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
+                                                      final double centerLon, final double centerLat, final double radiusMeters) {
+    return SloppyMath.haversin(centerLat, centerLon, rMinY, rMinX)*1000.0 <= radiusMeters
+        || SloppyMath.haversin(centerLat, centerLon, rMaxY, rMinX)*1000.0 <= radiusMeters
+        || SloppyMath.haversin(centerLat, centerLon, rMaxY, rMaxX)*1000.0 <= radiusMeters
+        || SloppyMath.haversin(centerLat, centerLon, rMinY, rMaxX)*1000.0 <= radiusMeters;
+  }
+
+  private static boolean rectAnyCornersOutsideCircle(final double rMinX, final double rMinY, final double rMaxX,
+                                                     final double rMaxY, final double centerLon, final double centerLat,
+                                                     final double radiusMeters, final boolean approx) {
+    if (approx == true) {
+      return rectAnyCornersOutsideCircleSloppy(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters);
+    }
+    double w = Math.abs(rMaxX - rMinX);
+    if (w <= 90.0) {
+      return GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, rMinX) > radiusMeters
+          || GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, rMinX) > radiusMeters
+          || GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, rMaxX) > radiusMeters
+          || GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, rMaxX) > radiusMeters;
+    }
+    // partition
+    w /= 4;
+    final double p1 = rMinX + w;
+    final double p2 = p1 + w;
+    final double p3 = p2 + w;
+
+    return GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, rMinX) > radiusMeters
+        || GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, rMinX) > radiusMeters
+        || GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, p1) > radiusMeters
+        || GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, p1) > radiusMeters
+        || GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, p2) > radiusMeters
+        || GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, p2) > radiusMeters
+        || GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, p3) > radiusMeters
+        || GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, p3) > radiusMeters
+        || GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, rMaxX) > radiusMeters
+        || GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, rMaxX) > radiusMeters;
+  }
+
+  private static boolean rectAnyCornersOutsideCircleSloppy(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
+                                                           final double centerLon, final double centerLat, final double radiusMeters) {
+    return SloppyMath.haversin(centerLat, centerLon, rMinY, rMinX)*1000.0 > radiusMeters
+        || SloppyMath.haversin(centerLat, centerLon, rMaxY, rMinX)*1000.0 > radiusMeters
+        || SloppyMath.haversin(centerLat, centerLon, rMaxY, rMaxX)*1000.0 > radiusMeters
+        || SloppyMath.haversin(centerLat, centerLon, rMinY, rMaxX)*1000.0 > radiusMeters;
+  }
+
+  public static boolean rectWithinCircle(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
+                                         final double centerLon, final double centerLat, final double radiusMeters) {
+    return rectWithinCircle(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters, false);
+  }
+
+  public static boolean rectWithinCircle(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
+                                         final double centerLon, final double centerLat, final double radiusMeters,
+                                         final boolean approx) {
+    return rectAnyCornersOutsideCircle(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters, approx) == false;
+  }
+
+  /**
+   * Determine if a bbox (defined by minLon, minLat, maxLon, maxLat) contains the provided point (defined by lon, lat)
+   * NOTE: this is basic method that does not handle dateline or pole crossing. Unwrapping must be done before
+   * calling this method.
+   */
+  public static boolean rectCrossesCircle(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
+                                          final double centerLon, final double centerLat, final double radiusMeters) {
+    return rectCrossesCircle(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters, false);
+  }
+
+  public static boolean rectCrossesCircle(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
+                                          final double centerLon, final double centerLat, final double radiusMeters,
+                                          final boolean approx) {
+    if (approx == true) {
+      return rectAnyCornersInCircle(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters, approx)
+          || isClosestPointOnRectWithinRange(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters, approx);
+    }
+
+    return (rectAnyCornersInCircle(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters, approx) &&
+        rectAnyCornersOutsideCircle(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters, approx))
+        || isClosestPointOnRectWithinRange(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters, approx);
+  }
+
+  private static boolean isClosestPointOnRectWithinRange(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
+                                                         final double centerLon, final double centerLat, final double radiusMeters,
+                                                         final boolean approx) {
+    double[] closestPt = {0, 0};
+    GeoDistanceUtils.closestPointOnBBox(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, closestPt);
+    boolean haverShortCut = GeoDistanceUtils.haversin(centerLat, centerLon, closestPt[1], closestPt[0]) <= radiusMeters;
+    if (approx == true || haverShortCut == true) {
+      return haverShortCut;
+    }
+    double lon1 = rMinX;
+    double lon2 = rMaxX;
+    double lat1 = rMinY;
+    double lat2 = rMaxY;
+    if (closestPt[0] == rMinX || closestPt[0] == rMaxX) {
+      lon1 = closestPt[0];
+      lon2 = lon1;
+    } else if (closestPt[1] == rMinY || closestPt[1] == rMaxY) {
+      lat1 = closestPt[1];
+      lat2 = lat1;
+    }
+
+    return lineCrossesSphere(lon1, lat1, 0, lon2, lat2, 0, centerLon, centerLat, 0, radiusMeters);
+  }
+
+  /**
+   * Computes whether or a 3dimensional line segment intersects or crosses a sphere
+   *
+   * @param lon1 longitudinal location of the line segment start point (in degrees)
+   * @param lat1 latitudinal location of the line segment start point (in degrees)
+   * @param alt1 altitude of the line segment start point (in degrees)
+   * @param lon2 longitudinal location of the line segment end point (in degrees)
+   * @param lat2 latitudinal location of the line segment end point (in degrees)
+   * @param alt2 altitude of the line segment end point (in degrees)
+   * @param centerLon longitudinal location of center search point (in degrees)
+   * @param centerLat latitudinal location of center search point (in degrees)
+   * @param centerAlt altitude of the center point (in meters)
+   * @param radiusMeters search sphere radius (in meters)
+   * @return whether the provided line segment is a secant of the
+   */
+  private static boolean lineCrossesSphere(double lon1, double lat1, double alt1, double lon2,
+                                           double lat2, double alt2, double centerLon, double centerLat,
+                                           double centerAlt, double radiusMeters) {
+    // convert to cartesian 3d (in meters)
+    double[] ecf1 = GeoProjectionUtils.llaToECF(lon1, lat1, alt1, null);
+    double[] ecf2 = GeoProjectionUtils.llaToECF(lon2, lat2, alt2, null);
+    double[] cntr = GeoProjectionUtils.llaToECF(centerLon, centerLat, centerAlt, null);
+
+    // convert radius from arc radius to cartesian radius
+    double[] oneEighty = GeoProjectionUtils.pointFromLonLatBearingGreatCircle(centerLon, centerLat, 180.0d, radiusMeters, new double[3]);
+    GeoProjectionUtils.llaToECF(oneEighty[0], oneEighty[1], 0, oneEighty);
+
+    radiusMeters = GeoDistanceUtils.linearDistance(oneEighty, cntr);//   Math.sqrt(oneEighty[0]*cntr[0] + oneEighty[1]*cntr[1] + oneEighty[2]*cntr[2]);
+
+    final double dX = ecf2[0] - ecf1[0];
+    final double dY = ecf2[1] - ecf1[1];
+    final double dZ = ecf2[2] - ecf1[2];
+    final double fX = ecf1[0] - cntr[0];
+    final double fY = ecf1[1] - cntr[1];
+    final double fZ = ecf1[2] - cntr[2];
+
+    final double a = dX*dX + dY*dY + dZ*dZ;
+    final double b = 2 * (fX*dX + fY*dY + fZ*dZ);
+    final double c = (fX*fX + fY*fY + fZ*fZ) - (radiusMeters*radiusMeters);
+
+    double discrim = (b*b)-(4*a*c);
+    if (discrim < 0) {
+      return false;
+    }
+
+    discrim = StrictMath.sqrt(discrim);
+    final double a2 = 2*a;
+    final double t1 = (-b - discrim)/a2;
+    final double t2 = (-b + discrim)/a2;
+
+    if ( (t1 < 0 || t1 > 1) ) {
+      return !(t2 < 0 || t2 > 1);
+    }
+
+    return true;
+  }
+}

Modified: lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/util/GeoUtils.java
URL: http://svn.apache.org/viewvc/lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/util/GeoUtils.java?rev=1718481&r1=1718480&r2=1718481&view=diff
==============================================================================
--- lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/util/GeoUtils.java (original)
+++ lucene/dev/trunk/lucene/sandbox/src/java/org/apache/lucene/util/GeoUtils.java Mon Dec  7 22:02:48 2015
@@ -19,6 +19,9 @@ package org.apache.lucene.util;
 
 import java.util.ArrayList;
 
+import static org.apache.lucene.util.SloppyMath.TO_DEGREES;
+import static org.apache.lucene.util.SloppyMath.TO_RADIANS;
+
 /**
  * Basic reusable geo-spatial utility methods
  *
@@ -74,6 +77,9 @@ public final class GeoUtils {
     return (val / LAT_SCALE) + MIN_LAT_INCL;
   }
 
+  /**
+   * Compare two position values within a {@link org.apache.lucene.util.GeoUtils#TOLERANCE} factor
+   */
   public static double compare(final double v1, final double v2) {
     final double delta = v1-v2;
     return Math.abs(delta) <= TOLERANCE ? 0 : delta;
@@ -107,44 +113,6 @@ public final class GeoUtils {
     return (off <= 180 ? off : 360-off) - 90;
   }
 
-  /**
-   * Determine if a bbox (defined by minLon, minLat, maxLon, maxLat) contains the provided point (defined by lon, lat)
-   * NOTE: this is a basic method that does not handle dateline or pole crossing. Unwrapping must be done before
-   * calling this method.
-   */
-  public static boolean bboxContains(final double lon, final double lat, final double minLon,
-                                     final double minLat, final double maxLon, final double maxLat) {
-    return (compare(lon, minLon) >= 0 && compare(lon, maxLon) <= 0
-          && compare(lat, minLat) >= 0 && compare(lat, maxLat) <= 0);
-  }
-
-  /**
-   * simple even-odd point in polygon computation
-   *    1.  Determine if point is contained in the longitudinal range
-   *    2.  Determine whether point crosses the edge by computing the latitudinal delta
-   *        between the end-point of a parallel vector (originating at the point) and the
-   *        y-component of the edge sink
-   *
-   * NOTE: Requires polygon point (x,y) order either clockwise or counter-clockwise
-   */
-  public static boolean pointInPolygon(double[] x, double[] y, double lat, double lon) {
-    assert x.length == y.length;
-    boolean inPoly = false;
-    /**
-     * Note: This is using a euclidean coordinate system which could result in
-     * upwards of 110KM error at the equator.
-     * TODO convert coordinates to cylindrical projection (e.g. mercator)
-     */
-    for (int i = 1; i < x.length; i++) {
-      if (x[i] < lon && x[i-1] >= lon || x[i-1] < lon && x[i] >= lon) {
-        if (y[i] + (lon - x[i]) / (x[i-1] - x[i]) * (y[i-1] - y[i]) < lat) {
-          inPoly = !inPoly;
-        }
-      }
-    }
-    return inPoly;
-  }
-
   public static String geoTermToString(long term) {
     StringBuilder s = new StringBuilder(64);
     final int numberOfLeadingZeros = Long.numberOfLeadingZeros(term);
@@ -157,95 +125,6 @@ public final class GeoUtils {
     return s.toString();
   }
 
-
-  public static boolean rectDisjoint(final double aMinX, final double aMinY, final double aMaxX, final double aMaxY,
-                                     final double bMinX, final double bMinY, final double bMaxX, final double bMaxY) {
-    return (aMaxX < bMinX || aMinX > bMaxX || aMaxY < bMinY || aMinY > bMaxY);
-  }
-
-  /**
-   * Computes whether the first (a) rectangle is wholly within another (b) rectangle (shared boundaries allowed)
-   */
-  public static boolean rectWithin(final double aMinX, final double aMinY, final double aMaxX, final double aMaxY,
-                                   final double bMinX, final double bMinY, final double bMaxX, final double bMaxY) {
-    return !(aMinX < bMinX || aMinY < bMinY || aMaxX > bMaxX || aMaxY > bMaxY);
-  }
-
-  public static boolean rectCrosses(final double aMinX, final double aMinY, final double aMaxX, final double aMaxY,
-                                    final double bMinX, final double bMinY, final double bMaxX, final double bMaxY) {
-    return !(rectDisjoint(aMinX, aMinY, aMaxX, aMaxY, bMinX, bMinY, bMaxX, bMaxY) ||
-        rectWithin(aMinX, aMinY, aMaxX, aMaxY, bMinX, bMinY, bMaxX, bMaxY));
-  }
-
-  /**
-   * Computes whether rectangle a contains rectangle b (touching allowed)
-   */
-  public static boolean rectContains(final double aMinX, final double aMinY, final double aMaxX, final double aMaxY,
-                                     final double bMinX, final double bMinY, final double bMaxX, final double bMaxY) {
-    return !(bMinX < aMinX || bMinY < aMinY || bMaxX > aMaxX || bMaxY > aMaxY);
-  }
-
-  /**
-   * Computes whether a rectangle intersects another rectangle (crosses, within, touching, etc)
-   */
-  public static boolean rectIntersects(final double aMinX, final double aMinY, final double aMaxX, final double aMaxY,
-                                       final double bMinX, final double bMinY, final double bMaxX, final double bMaxY) {
-    return !((aMaxX < bMinX || aMinX > bMaxX || aMaxY < bMinY || aMinY > bMaxY) );
-  }
-
-  /**
-   * Computes whether a rectangle crosses a shape. (touching not allowed)
-   */
-  public static boolean rectCrossesPoly(final double rMinX, final double rMinY, final double rMaxX,
-                                        final double rMaxY, final double[] shapeX, final double[] shapeY,
-                                        final double sMinX, final double sMinY, final double sMaxX,
-                                        final double sMaxY) {
-    // short-circuit: if the bounding boxes are disjoint then the shape does not cross
-    if (rectDisjoint(rMinX, rMinY, rMaxX, rMaxY, sMinX, sMinY, sMaxX, sMaxY)) {
-      return false;
-    }
-
-    final double[][] bbox = new double[][] { {rMinX, rMinY}, {rMaxX, rMinY}, {rMaxX, rMaxY}, {rMinX, rMaxY}, {rMinX, rMinY} };
-    final int polyLength = shapeX.length-1;
-    double d, s, t, a1, b1, c1, a2, b2, c2;
-    double x00, y00, x01, y01, x10, y10, x11, y11;
-
-    // computes the intersection point between each bbox edge and the polygon edge
-    for (short b=0; b<4; ++b) {
-      a1 = bbox[b+1][1]-bbox[b][1];
-      b1 = bbox[b][0]-bbox[b+1][0];
-      c1 = a1*bbox[b+1][0] + b1*bbox[b+1][1];
-      for (int p=0; p<polyLength; ++p) {
-        a2 = shapeY[p+1]-shapeY[p];
-        b2 = shapeX[p]-shapeX[p+1];
-        // compute determinant
-        d = a1*b2 - a2*b1;
-        if (d != 0) {
-          // lines are not parallel, check intersecting points
-          c2 = a2*shapeX[p+1] + b2*shapeY[p+1];
-          s = (1/d)*(b2*c1 - b1*c2);
-          t = (1/d)*(a1*c2 - a2*c1);
-          x00 = StrictMath.min(bbox[b][0], bbox[b+1][0]) - TOLERANCE;
-          x01 = StrictMath.max(bbox[b][0], bbox[b+1][0]) + TOLERANCE;
-          y00 = StrictMath.min(bbox[b][1], bbox[b+1][1]) - TOLERANCE;
-          y01 = StrictMath.max(bbox[b][1], bbox[b+1][1]) + TOLERANCE;
-          x10 = StrictMath.min(shapeX[p], shapeX[p+1]) - TOLERANCE;
-          x11 = StrictMath.max(shapeX[p], shapeX[p+1]) + TOLERANCE;
-          y10 = StrictMath.min(shapeY[p], shapeY[p+1]) - TOLERANCE;
-          y11 = StrictMath.max(shapeY[p], shapeY[p+1]) + TOLERANCE;
-          // check whether the intersection point is touching one of the line segments
-          boolean touching = ((x00 == s && y00 == t) || (x01 == s && y01 == t))
-              || ((x10 == s && y10 == t) || (x11 == s && y11 == t));
-          // if line segments are not touching and the intersection point is within the range of either segment
-          if (!(touching || x00 > s || x01 < s || y00 > t || y01 < t || x10 > s || x11 < s || y10 > t || y11 < t)) {
-            return true;
-          }
-        }
-      } // for each poly edge
-    } // for each bbox edge
-    return false;
-  }
-
   /**
    * Converts a given circle (defined as a point/radius) to an approximated line-segment polygon
    *
@@ -267,7 +146,7 @@ public final class GeoUtils {
     final int sidesLen = sides-1;
     for (int i=0; i<sidesLen; ++i) {
       angle = (i*360/sides);
-      pt = GeoProjectionUtils.pointFromLonLatBearing(lon, lat, angle, radiusMeters, pt);
+      pt = GeoProjectionUtils.pointFromLonLatBearingGreatCircle(lon, lat, angle, radiusMeters, pt);
       lons[i] = pt[0];
       lats[i] = pt[1];
     }
@@ -281,63 +160,11 @@ public final class GeoUtils {
   }
 
   /**
-   * Computes whether a rectangle is within a given polygon (shared boundaries allowed)
-   */
-  public static boolean rectWithinPoly(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
-                                       final double[] shapeX, final double[] shapeY, final double sMinX,
-                                       final double sMinY, final double sMaxX, final double sMaxY) {
-    // check if rectangle crosses poly (to handle concave/pacman polys), then check that all 4 corners
-    // are contained
-    return !(rectCrossesPoly(rMinX, rMinY, rMaxX, rMaxY, shapeX, shapeY, sMinX, sMinY, sMaxX, sMaxY) ||
-        !pointInPolygon(shapeX, shapeY, rMinY, rMinX) || !pointInPolygon(shapeX, shapeY, rMinY, rMaxX) ||
-        !pointInPolygon(shapeX, shapeY, rMaxY, rMaxX) || !pointInPolygon(shapeX, shapeY, rMaxY, rMinX));
-  }
-
-  private static boolean rectAnyCornersOutsideCircle(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
-                                                     final double centerLon, final double centerLat, final double radiusMeters) {
-    return SloppyMath.haversin(centerLat, centerLon, rMinY, rMinX)*1000.0 > radiusMeters
-      || SloppyMath.haversin(centerLat, centerLon, rMaxY, rMinX)*1000.0 > radiusMeters
-      || SloppyMath.haversin(centerLat, centerLon, rMaxY, rMaxX)*1000.0 > radiusMeters
-      || SloppyMath.haversin(centerLat, centerLon, rMinY, rMaxX)*1000.0 > radiusMeters;
-  }
-
-  private static boolean rectAnyCornersInCircle(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
-                                                final double centerLon, final double centerLat, final double radiusMeters) {
-    return SloppyMath.haversin(centerLat, centerLon, rMinY, rMinX)*1000.0 <= radiusMeters
-        || SloppyMath.haversin(centerLat, centerLon, rMaxY, rMinX)*1000.0 <= radiusMeters
-        || SloppyMath.haversin(centerLat, centerLon, rMaxY, rMaxX)*1000.0 <= radiusMeters
-        || SloppyMath.haversin(centerLat, centerLon, rMinY, rMaxX)*1000.0 <= radiusMeters;
-  }
-
-  public static boolean rectWithinCircle(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
-                                         final double centerLon, final double centerLat, final double radiusMeters) {
-    return rectAnyCornersOutsideCircle(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters) == false;
-  }
-
-  /**
-   * Determine if a bbox (defined by minLon, minLat, maxLon, maxLat) contains the provided point (defined by lon, lat)
-   * NOTE: this is basic method that does not handle dateline or pole crossing. Unwrapping must be done before
-   * calling this method.
-   */
-  public static boolean rectCrossesCircle(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
-                                          final double centerLon, final double centerLat, final double radiusMeters) {
-    return rectAnyCornersInCircle(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters)
-        || isClosestPointOnRectWithinRange(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters);
-  }
-
-  private static boolean isClosestPointOnRectWithinRange(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
-                                                         final double centerLon, final double centerLat, final double radiusMeters) {
-    double[] closestPt = {0, 0};
-    GeoDistanceUtils.closestPointOnBBox(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, closestPt);
-    return SloppyMath.haversin(centerLat, centerLon, closestPt[1], closestPt[0])*1000.0 <= radiusMeters;
-  }
-
-  /**
    * Compute Bounding Box for a circle using WGS-84 parameters
    */
   public static GeoRect circleToBBox(final double centerLon, final double centerLat, final double radiusMeters) {
-    final double radLat = StrictMath.toRadians(centerLat);
-    final double radLon = StrictMath.toRadians(centerLon);
+    final double radLat = TO_RADIANS * centerLat;
+    final double radLon = TO_RADIANS * centerLon;
     double radDistance = radiusMeters / GeoProjectionUtils.SEMIMAJOR_AXIS;
     double minLat = radLat - radDistance;
     double maxLat = radLat + radDistance;
@@ -345,7 +172,7 @@ public final class GeoUtils {
     double maxLon;
 
     if (minLat > GeoProjectionUtils.MIN_LAT_RADIANS && maxLat < GeoProjectionUtils.MAX_LAT_RADIANS) {
-      double deltaLon = StrictMath.asin(StrictMath.sin(radDistance) / StrictMath.cos(radLat));
+      double deltaLon = SloppyMath.asin(SloppyMath.sin(radDistance) / SloppyMath.cos(radLat));
       minLon = radLon - deltaLon;
       if (minLon < GeoProjectionUtils.MIN_LON_RADIANS) {
         minLon += 2d * StrictMath.PI;
@@ -362,10 +189,12 @@ public final class GeoUtils {
       maxLon = GeoProjectionUtils.MAX_LON_RADIANS;
     }
 
-    return new GeoRect(StrictMath.toDegrees(minLon), StrictMath.toDegrees(maxLon),
-        StrictMath.toDegrees(minLat), StrictMath.toDegrees(maxLat));
+    return new GeoRect(TO_DEGREES * minLon, TO_DEGREES * maxLon, TO_DEGREES * minLat, TO_DEGREES * maxLat);
   }
 
+  /**
+   * Compute Bounding Box for a polygon using WGS-84 parameters
+   */
   public static GeoRect polyToBBox(double[] polyLons, double[] polyLats) {
     if (polyLons.length != polyLats.length) {
       throw new IllegalArgumentException("polyLons and polyLats must be equal length");
@@ -393,64 +222,16 @@ public final class GeoUtils {
         GeoUtils.unscaleLat(GeoUtils.scaleLat(minLat)), GeoUtils.unscaleLat(GeoUtils.scaleLat(maxLat)));
   }
 
-  /*
   /**
-   * Computes whether or a 3dimensional line segment intersects or crosses a sphere
-   *
-   * @param lon1 longitudinal location of the line segment start point (in degrees)
-   * @param lat1 latitudinal location of the line segment start point (in degrees)
-   * @param alt1 altitude of the line segment start point (in degrees)
-   * @param lon2 longitudinal location of the line segment end point (in degrees)
-   * @param lat2 latitudinal location of the line segment end point (in degrees)
-   * @param alt2 altitude of the line segment end point (in degrees)
-   * @param centerLon longitudinal location of center search point (in degrees)
-   * @param centerLat latitudinal location of center search point (in degrees)
-   * @param centerAlt altitude of the center point (in meters)
-   * @param radiusMeters search sphere radius (in meters)
-   * @return whether the provided line segment is a secant of the
-   * /
-  // NOTE: not used for 2d at the moment. used for 3d w/ altitude (we can keep or add back)
-  private static boolean lineCrossesSphere(double lon1, double lat1, double alt1, double lon2,
-                                           double lat2, double alt2, double centerLon, double centerLat,
-                                           double centerAlt, double radiusMeters) {
-    // convert to cartesian 3d (in meters)
-    double[] ecf1 = GeoProjectionUtils.llaToECF(lon1, lat1, alt1, null);
-    double[] ecf2 = GeoProjectionUtils.llaToECF(lon2, lat2, alt2, null);
-    double[] cntr = GeoProjectionUtils.llaToECF(centerLon, centerLat, centerAlt, null);
-
-    final double dX = ecf2[0] - ecf1[0];
-    final double dY = ecf2[1] - ecf1[1];
-    final double dZ = ecf2[2] - ecf1[2];
-    final double fX = ecf1[0] - cntr[0];
-    final double fY = ecf1[1] - cntr[1];
-    final double fZ = ecf1[2] - cntr[2];
-
-    final double a = dX*dX + dY*dY + dZ*dZ;
-    final double b = 2 * (fX*dX + fY*dY + fZ*dZ);
-    final double c = (fX*fX + fY*fY + fZ*fZ) - (radiusMeters*radiusMeters);
-
-    double discrim = (b*b)-(4*a*c);
-    if (discrim < 0) {
-      return false;
-    }
-
-    discrim = StrictMath.sqrt(discrim);
-    final double a2 = 2*a;
-    final double t1 = (-b - discrim)/a2;
-    final double t2 = (-b + discrim)/a2;
-
-    if ( (t1 < 0 || t1 > 1) ) {
-      return !(t2 < 0 || t2 > 1);
-    }
-
-    return true;
-  }
-  */
-
+   * validates latitude value is within standard +/-90 coordinate bounds
+   */
   public static boolean isValidLat(double lat) {
     return Double.isNaN(lat) == false && lat >= MIN_LAT_INCL && lat <= MAX_LAT_INCL;
   }
 
+  /**
+   * validates longitude value is within standard +/-180 coordinate bounds
+   */
   public static boolean isValidLon(double lon) {
     return Double.isNaN(lon) == false && lon >= MIN_LON_INCL && lon <= MAX_LON_INCL;
   }

Modified: lucene/dev/trunk/lucene/sandbox/src/test/org/apache/lucene/search/TestDimensionalQueries.java
URL: http://svn.apache.org/viewvc/lucene/dev/trunk/lucene/sandbox/src/test/org/apache/lucene/search/TestDimensionalQueries.java?rev=1718481&r1=1718480&r2=1718481&view=diff
==============================================================================
--- lucene/dev/trunk/lucene/sandbox/src/test/org/apache/lucene/search/TestDimensionalQueries.java (original)
+++ lucene/dev/trunk/lucene/sandbox/src/test/org/apache/lucene/search/TestDimensionalQueries.java Mon Dec  7 22:02:48 2015
@@ -21,6 +21,7 @@ import org.apache.lucene.document.Dimens
 import org.apache.lucene.document.Document;
 import org.apache.lucene.store.Directory;
 import org.apache.lucene.util.BaseGeoPointTestCase;
+import org.apache.lucene.util.GeoDistanceUtils;
 import org.apache.lucene.util.GeoRect;
 import org.apache.lucene.util.SloppyMath;
 
@@ -96,15 +97,15 @@ public class TestDimensionalQueries exte
 
   @Override
   protected Boolean circleContainsPoint(double centerLat, double centerLon, double radiusMeters, double pointLat, double pointLon) {
-    double distanceKM = SloppyMath.haversin(centerLat, centerLon, pointLat, pointLon);
-    boolean result = distanceKM*1000.0 <= radiusMeters;
+    double distanceMeters = GeoDistanceUtils.haversin(centerLat, centerLon, pointLat, pointLon);
+    boolean result = distanceMeters <= radiusMeters;
     //System.out.println("  shouldMatch?  centerLon=" + centerLon + " centerLat=" + centerLat + " pointLon=" + pointLon + " pointLat=" + pointLat + " result=" + result + " distanceMeters=" + (distanceKM * 1000));
     return result;
   }
 
   @Override
   protected Boolean distanceRangeContainsPoint(double centerLat, double centerLon, double minRadiusMeters, double radiusMeters, double pointLat, double pointLon) {
-    final double d = SloppyMath.haversin(centerLat, centerLon, pointLat, pointLon)*1000.0;
+    final double d = GeoDistanceUtils.haversin(centerLat, centerLon, pointLat, pointLon);
     return d >= minRadiusMeters && d <= radiusMeters;
   }
 

Modified: lucene/dev/trunk/lucene/sandbox/src/test/org/apache/lucene/search/TestGeoPointQuery.java
URL: http://svn.apache.org/viewvc/lucene/dev/trunk/lucene/sandbox/src/test/org/apache/lucene/search/TestGeoPointQuery.java?rev=1718481&r1=1718480&r2=1718481&view=diff
==============================================================================
--- lucene/dev/trunk/lucene/sandbox/src/test/org/apache/lucene/search/TestGeoPointQuery.java (original)
+++ lucene/dev/trunk/lucene/sandbox/src/test/org/apache/lucene/search/TestGeoPointQuery.java Mon Dec  7 22:02:48 2015
@@ -28,6 +28,7 @@ import org.apache.lucene.index.RandomInd
 import org.apache.lucene.store.Directory;
 import org.apache.lucene.util.BaseGeoPointTestCase;
 import org.apache.lucene.util.GeoRect;
+import org.apache.lucene.util.GeoRelationUtils;
 import org.apache.lucene.util.GeoUtils;
 import org.apache.lucene.util.SloppyMath;
 import org.apache.lucene.util.TestUtil;
@@ -171,11 +172,11 @@ public class TestGeoPointQuery extends B
     }
 
     if (rect.minLon < rect.maxLon) {
-      return GeoUtils.bboxContains(pointLon, pointLat, rect.minLon, rect.minLat, rect.maxLon, rect.maxLat);
+      return GeoRelationUtils.pointInRect(pointLon, pointLat, rect.minLon, rect.minLat, rect.maxLon, rect.maxLat);
     } else {
       // Rect crosses dateline:
-      return GeoUtils.bboxContains(pointLon, pointLat, -180.0, rect.minLat, rect.maxLon, rect.maxLat)
-          || GeoUtils.bboxContains(pointLon, pointLat, rect.minLon, rect.minLat, 180.0, rect.maxLat);
+      return GeoRelationUtils.pointInRect(pointLon, pointLat, -180.0, rect.minLat, rect.maxLon, rect.maxLat)
+          || GeoRelationUtils.pointInRect(pointLon, pointLat, rect.minLon, rect.minLat, 180.0, rect.maxLat);
     }
   }
 
@@ -221,7 +222,7 @@ public class TestGeoPointQuery extends B
   }
 
   public void testRectCrossesCircle() throws Exception {
-    assertTrue(GeoUtils.rectCrossesCircle(-180, -90, 180, 0.0, 0.667, 0.0, 88000.0));
+    assertTrue(GeoRelationUtils.rectCrossesCircle(-180, -90, 180, 0.0, 0.667, 0.0, 88000.0));
   }
 
   private TopDocs geoDistanceRangeQuery(double lon, double lat, double minRadius, double maxRadius, int limit)
@@ -261,9 +262,9 @@ public class TestGeoPointQuery extends B
     double yMax = 1;//5;
 
     // test cell crossing poly
-    assertTrue(GeoUtils.rectCrossesPoly(xMin, yMin, xMax, yMax, px, py, xMinA, yMinA, xMaxA, yMaxA));
-    assertFalse(GeoUtils.rectCrossesPoly(-5, 0,  0.000001, 5, px, py, xMin, yMin, xMax, yMax));
-    assertTrue(GeoUtils.rectWithinPoly(-5, 0, -2, 5, px, py, xMin, yMin, xMax, yMax));
+    assertTrue(GeoRelationUtils.rectCrossesPoly(xMin, yMin, xMax, yMax, px, py, xMinA, yMinA, xMaxA, yMaxA));
+    assertFalse(GeoRelationUtils.rectCrossesPoly(-5, 0,  0.000001, 5, px, py, xMin, yMin, xMax, yMax));
+    assertTrue(GeoRelationUtils.rectWithinPoly(-5, 0, -2, 5, px, py, xMin, yMin, xMax, yMax));
   }
 
   public void testBBoxCrossDateline() throws Exception {

Modified: lucene/dev/trunk/lucene/sandbox/src/test/org/apache/lucene/util/TestGeoUtils.java
URL: http://svn.apache.org/viewvc/lucene/dev/trunk/lucene/sandbox/src/test/org/apache/lucene/util/TestGeoUtils.java?rev=1718481&r1=1718480&r2=1718481&view=diff
==============================================================================
--- lucene/dev/trunk/lucene/sandbox/src/test/org/apache/lucene/util/TestGeoUtils.java (original)
+++ lucene/dev/trunk/lucene/sandbox/src/test/org/apache/lucene/util/TestGeoUtils.java Mon Dec  7 22:02:48 2015
@@ -312,7 +312,7 @@ public class TestGeoUtils extends Lucene
         // Leaf cell: brute force check all docs that fall within this cell:
         for(int docID=0;docID<docLons.length;docID++) {
           if (cell.contains(docLons[docID], docLats[docID])) {
-            double distanceMeters = SloppyMath.haversin(centerLat, centerLon, docLats[docID], docLons[docID]) * 1000.0;
+            double distanceMeters = GeoDistanceUtils.haversin(centerLat, centerLon, docLats[docID], docLons[docID]);
             if (distanceMeters <= radiusMeters) {
               if (VERBOSE) {
                 log.println("    check doc=" + docID + ": match!");
@@ -327,7 +327,7 @@ public class TestGeoUtils extends Lucene
         }
       } else {
 
-        if (GeoUtils.rectWithinCircle(cell.minLon, cell.minLat, cell.maxLon, cell.maxLat, centerLon, centerLat, radiusMeters)) {
+        if (GeoRelationUtils.rectWithinCircle(cell.minLon, cell.minLat, cell.maxLon, cell.maxLat, centerLon, centerLat, radiusMeters)) {
           // Query circle fully contains this cell, just addAll:
           if (VERBOSE) {
             log.println("    circle fully contains cell: now addAll");
@@ -341,13 +341,13 @@ public class TestGeoUtils extends Lucene
             }
           }
           continue;
-        } else if (GeoUtils.rectWithin(root.minLon, root.minLat, root.maxLon, root.maxLat,
+        } else if (GeoRelationUtils.rectWithin(root.minLon, root.minLat, root.maxLon, root.maxLat,
                                        cell.minLon, cell.minLat, cell.maxLon, cell.maxLat)) {
           // Fall through below to "recurse"
           if (VERBOSE) {
             log.println("    cell fully contains circle: keep splitting");
           }
-        } else if (GeoUtils.rectCrossesCircle(cell.minLon, cell.minLat, cell.maxLon, cell.maxLat,
+        } else if (GeoRelationUtils.rectCrossesCircle(cell.minLon, cell.minLat, cell.maxLon, cell.maxLat,
                                               centerLon, centerLat, radiusMeters)) {
           // Fall through below to "recurse"
           if (VERBOSE) {
@@ -415,15 +415,11 @@ public class TestGeoUtils extends Lucene
   }
 
   /** Tests consistency of GeoUtils.rectWithinCircle, .rectCrossesCircle, .rectWithin and SloppyMath.haversine distance check */
-  @AwaitsFix(bugUrl = "https://issues.apache.org/jira/browse/LUCENE-6846")
   public void testGeoRelations() throws Exception {
 
     int numDocs = atLeast(1000);
     
-    // boolean useSmallRanges = random().nextBoolean();
-
-    // TODO: the GeoUtils APIs have bugs if you use large distances:
-    boolean useSmallRanges = true;
+    boolean useSmallRanges = random().nextBoolean();
 
     if (VERBOSE) {
       System.out.println("TEST: " + numDocs + " docs useSmallRanges=" + useSmallRanges);
@@ -475,22 +471,29 @@ public class TestGeoUtils extends Lucene
 
       if (bbox.maxLon < bbox.minLon) {
         // Crosses dateline
-        log.println("  circle crosses dateline; first right query");
+        log.println("  circle crosses dateline; first left query");
+        double unwrappedLon = centerLon;
+        if (unwrappedLon > bbox.maxLon) {
+          // unwrap left
+          unwrappedLon += -360.0D;
+        }
         findMatches(hits, log,
                     new Cell(null,
                              -180, bbox.minLat,
                              bbox.maxLon, bbox.maxLat,
                              0),
-                    centerLon, centerLat, radiusMeters,
-                    docLons, docLats);
-        log.println("  circle crosses dateline; now left query");
+                    unwrappedLon, centerLat, radiusMeters, docLons, docLats);
+        log.println("  circle crosses dateline; now right query");
+        if (unwrappedLon < bbox.maxLon) {
+          // unwrap right
+          unwrappedLon += 360.0D;
+        }
         findMatches(hits, log,
                     new Cell(null,
                              bbox.minLon, bbox.minLat,
                              180, bbox.maxLat,
                              0),
-                    centerLon, centerLat, radiusMeters,
-                    docLons, docLats);
+                    unwrappedLon, centerLat, radiusMeters, docLons, docLats);
       } else {
         // Start with the root cell that fully contains the shape:
         findMatches(hits, log,
@@ -510,15 +513,15 @@ public class TestGeoUtils extends Lucene
 
       // Done matching, now verify:
       for(int docID=0;docID<numDocs;docID++) {
-        double distanceMeters = SloppyMath.haversin(centerLat, centerLon, docLats[docID], docLons[docID]) * 1000.0;
+        double distanceMeters = GeoDistanceUtils.haversin(centerLat, centerLon, docLats[docID], docLons[docID]);
         boolean expected = distanceMeters <= radiusMeters;
 
         boolean actual = hits.contains(docID);
         if (actual != expected) {
           if (actual) {
-            log.println("doc=" + docID + " matched but should not");
+            log.println("doc=" + docID + " matched but should not on iteration " + iter);
           } else {
-            log.println("doc=" + docID + " did not match but should");
+            log.println("doc=" + docID + " did not match but should on iteration " + iter);
           }
           log.println("  lon=" + docLons[docID] + " lat=" + docLats[docID] + " distanceMeters=" + distanceMeters + " vs radiusMeters=" + radiusMeters);
           failCount++;



Mime
View raw message