LUCENE-7159: Speed up LatLonPoint point-in-polygon performance

lucene/CHANGES.txt

@@ -56,6 +56,9 @@ Optimizations
   multiple polygons and holes, with memory usage independent of
   polygon complexity. (Karl Wright, Mike McCandless, Robert Muir)
 
+* LUCENE-7159: Speed up LatLonPoint polygon performance for complex
+  polygons. (Robert Muir)
+
 Bug Fixes
 
 * LUCENE-7127: Fix corner case bugs in GeoPointDistanceQuery. (Robert Muir)

lucene/core/src/java/org/apache/lucene/geo/Polygon.java

@@ -18,6 +18,8 @@ package org.apache.lucene.geo;
 
 import java.util.Arrays;
 
+import org.apache.lucene.index.PointValues.Relation;
+
 /**
  * Represents a closed polygon on the earth's surface.
  * <p>
@@ -151,89 +153,126 @@ public final class Polygon {
     }
   }
   
-  /**
+  /** Returns relation to the provided rectangle */
-   * Computes whether a rectangle is within a polygon (shared boundaries not allowed)
+  public Relation relate(double minLat, double maxLat, double minLon, double maxLon) {
-   */
-  public boolean contains(double minLat, double maxLat, double minLon, double maxLon) {
-    // check if rectangle crosses poly (to handle concave/pacman polys), then check that all 4 corners
-    // are contained
-    boolean contains = crosses(minLat, maxLat, minLon, maxLon) == false &&
-                       contains(minLat, minLon) &&
-                       contains(minLat, maxLon) &&
-                       contains(maxLat, maxLon) &&
-                       contains(maxLat, minLon);
-
-    if (contains) {
-      // if we intersect with any hole, game over
-      for (Polygon hole : holes) {
-        if (hole.crosses(minLat, maxLat, minLon, maxLon) || hole.contains(minLat, maxLat, minLon, maxLon)) {
-          return false;
-        }
-      }
-      return true;
-    } else {
-      return false;
-    }
-  }
-
-  /**
-   * Convenience method for accurately computing whether a rectangle crosses a poly.
-   */
-  public boolean crosses(double minLat, double maxLat, final double minLon, final double maxLon) {
     // if the bounding boxes are disjoint then the shape does not cross
     if (maxLon < this.minLon || minLon > this.maxLon || maxLat < this.minLat || minLat > this.maxLat) {
-      return false;
+      return Relation.CELL_OUTSIDE_QUERY;
     }
     // if the rectangle fully encloses us, we cross.
     if (minLat <= this.minLat && maxLat >= this.maxLat && minLon <= this.minLon && maxLon >= this.maxLon) {
-      return true;
+      return Relation.CELL_CROSSES_QUERY;
     }
-    // if we cross any hole, we cross
+    // check any holes
     for (Polygon hole : holes) {
-      if (hole.crosses(minLat, maxLat, minLon, maxLon)) {
+      Relation holeRelation = hole.relate(minLat, maxLat, minLon, maxLon);
-        return true;
+      if (holeRelation == Relation.CELL_CROSSES_QUERY) {
+        return Relation.CELL_CROSSES_QUERY;
+      } else if (holeRelation == Relation.CELL_INSIDE_QUERY) {
+        return Relation.CELL_OUTSIDE_QUERY;
       }
     }
+    // check each corner: if < 4 are present, its cheaper than crossesSlowly
+    int numCorners = numberOfCorners(minLat, maxLat, minLon, maxLon);
+    if (numCorners == 4) {
+      if (crossesSlowly(minLat, maxLat, minLon, maxLon)) {
+        return Relation.CELL_CROSSES_QUERY;
+      }
+      return Relation.CELL_INSIDE_QUERY;
+    } else if (numCorners > 0) {
+      return Relation.CELL_CROSSES_QUERY;
+    }
     
+    // we cross
+    if (crossesSlowly(minLat, maxLat, minLon, maxLon)) {
+      return Relation.CELL_CROSSES_QUERY;
+    }
+    
+    return Relation.CELL_OUTSIDE_QUERY;
+  }
+  
+  // returns 0, 4, or something in between
+  private int numberOfCorners(double minLat, double maxLat, double minLon, double maxLon) {
+    int containsCount = 0;
+    if (contains(minLat, minLon)) {
+      containsCount++;
+    }
+    if (contains(minLat, maxLon)) {
+      containsCount++;
+    }
+    if (containsCount == 1) {
+      return containsCount;
+    }
+    if (contains(maxLat, maxLon)) {
+      containsCount++;
+    }
+    if (containsCount == 2) {
+      return containsCount;
+    }
+    if (contains(maxLat, minLon)) {
+      containsCount++;
+    }
+    return containsCount;
+  }
+
+  private boolean crossesSlowly(double minLat, double maxLat, final double minLon, final double maxLon) {
     /*
      * Accurately compute (within restrictions of cartesian decimal degrees) whether a rectangle crosses a polygon
      */
-    final double[][] bbox = new double[][] { {minLon, minLat}, {maxLon, minLat}, {maxLon, maxLat}, {minLon, maxLat}, {minLon, minLat} };
+    final double[] boxLats = new double[] { minLat, minLat, maxLat, maxLat, minLat };
-    final int polyLength = polyLons.length-1;
+    final double[] boxLons = new double[] { minLon, maxLon, maxLon, minLon, minLon };
-    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) {
+    for (int b=0; b<4; ++b) {
-      a1 = bbox[b+1][1]-bbox[b][1];
+      double a1 = boxLats[b+1]-boxLats[b];
-      b1 = bbox[b][0]-bbox[b+1][0];
+      double b1 = boxLons[b]-boxLons[b+1];
-      c1 = a1*bbox[b+1][0] + b1*bbox[b+1][1];
+      double c1 = a1*boxLons[b+1] + b1*boxLats[b+1];
-      for (int p=0; p<polyLength; ++p) {
+      for (int p=0; p<polyLons.length-1; ++p) {
-        a2 = polyLats[p+1]-polyLats[p];
+        double a2 = polyLats[p+1]-polyLats[p];
-        b2 = polyLons[p]-polyLons[p+1];
+        double b2 = polyLons[p]-polyLons[p+1];
         // compute determinant
-        d = a1*b2 - a2*b1;
+        double d = a1*b2 - a2*b1;
         if (d != 0) {
           // lines are not parallel, check intersecting points
-          c2 = a2*polyLons[p+1] + b2*polyLats[p+1];
+          double c2 = a2*polyLons[p+1] + b2*polyLats[p+1];
-          s = (1/d)*(b2*c1 - b1*c2);
+          double s = (1/d)*(b2*c1 - b1*c2);
-          t = (1/d)*(a1*c2 - a2*c1);
           // todo TOLERANCE SHOULD MATCH EVERYWHERE this is currently blocked by LUCENE-7165
-          x00 = Math.min(bbox[b][0], bbox[b+1][0]) - ENCODING_TOLERANCE;
+          double x00 = Math.min(boxLons[b], boxLons[b+1]) - ENCODING_TOLERANCE;
-          x01 = Math.max(bbox[b][0], bbox[b+1][0]) + ENCODING_TOLERANCE;
+          if (x00 > s) {
-          y00 = Math.min(bbox[b][1], bbox[b+1][1]) - ENCODING_TOLERANCE;
+            continue; // out of range
-          y01 = Math.max(bbox[b][1], bbox[b+1][1]) + ENCODING_TOLERANCE;
-          x10 = Math.min(polyLons[p], polyLons[p+1]) - ENCODING_TOLERANCE;
-          x11 = Math.max(polyLons[p], polyLons[p+1]) + ENCODING_TOLERANCE;
-          y10 = Math.min(polyLats[p], polyLats[p+1]) - ENCODING_TOLERANCE;
-          y11 = Math.max(polyLats[p], polyLats[p+1]) + ENCODING_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;
           }
+          double x01 = Math.max(boxLons[b], boxLons[b+1]) + ENCODING_TOLERANCE;
+          if (x01 < s) {
+            continue; // out of range
+          }
+          double x10 = Math.min(polyLons[p], polyLons[p+1]) - ENCODING_TOLERANCE;
+          if (x10 > s) {
+            continue; // out of range
+          }
+          double x11 = Math.max(polyLons[p], polyLons[p+1]) + ENCODING_TOLERANCE;
+          if (x11 < s) {
+            continue; // out of range
+          }
+
+          double t = (1/d)*(a1*c2 - a2*c1);
+          double y00 = Math.min(boxLats[b], boxLats[b+1]) - ENCODING_TOLERANCE;
+          if (y00 > t || (x00 == s && y00 == t)) {
+            continue; // out of range or touching
+          }
+          double y01 = Math.max(boxLats[b], boxLats[b+1]) + ENCODING_TOLERANCE;
+          if (y01 < t || (x01 == s && y01 == t)) {
+            continue; // out of range or touching
+          }
+          double y10 = Math.min(polyLats[p], polyLats[p+1]) - ENCODING_TOLERANCE;
+          if (y10 > t || (x10 == s && y10 == t)) {
+            continue; // out of range or touching
+          }
+          double y11 = Math.max(polyLats[p], polyLats[p+1]) + ENCODING_TOLERANCE;
+          if (y11 < t || (x11 == s && y11 == t)) {
+            continue; // out of range or touching
+          }
+          // if line segments are not touching and the intersection point is within the range of either segment
+          return true;
         }
       } // for each poly edge
     } // for each bbox edge
@@ -265,24 +304,17 @@ public final class Polygon {
     return false;
   }
 
-  /** Helper for multipolygon logic: returns true if any of the supplied polygons contain the rectangle */
+  /** Returns the multipolygon relation for the rectangle */
-  public static boolean contains(Polygon[] polygons, double minLat, double maxLat, double minLon, double maxLon) {
+  public static Relation relate(Polygon[] polygons, double minLat, double maxLat, double minLon, double maxLon) {
     for (Polygon polygon : polygons) {
-      if (polygon.contains(minLat, maxLat, minLon, maxLon)) {
+      Relation relation = polygon.relate(minLat, maxLat, minLon, maxLon);
-        return true;
+      if (relation != Relation.CELL_OUTSIDE_QUERY) {
+        // note: we optimize for non-overlapping multipolygons. so if we cross one,
+        // we won't keep iterating to try to find a contains.
+        return relation;
       }
     }
-    return false;
+    return Relation.CELL_OUTSIDE_QUERY;
-  }
-
-  /** Helper for multipolygon logic: returns true if any of the supplied polygons crosses the rectangle */
-  public static boolean crosses(Polygon[] polygons, double minLat, double maxLat, double minLon, double maxLon) {
-    for (Polygon polygon : polygons) {
-      if (polygon.crosses(minLat, maxLat, minLon, maxLon)) {
-        return true;
-      }
-    }
-    return false;
   }
 
   @Override

lucene/core/src/test/org/apache/lucene/geo/TestPolygon.java

@@ -17,6 +17,7 @@
 package org.apache.lucene.geo;
 
 import org.apache.lucene.geo.Polygon;
+import org.apache.lucene.index.PointValues.Relation;
 import org.apache.lucene.util.LuceneTestCase;
 
 import static org.apache.lucene.geo.GeoTestUtil.nextLatitude;
@@ -80,21 +81,15 @@ public class TestPolygon extends LuceneTestCase {
     assertTrue(Polygon.contains(polygons, -25, 25)); // on the mainland
     assertFalse(Polygon.contains(polygons, -51, 51)); // in the ocean
     
-    // contains(box): this can conservatively return false
+    // relate(box): this can conservatively return CELL_CROSSES_QUERY
-    assertTrue(Polygon.contains(polygons, -2, 2, -2, 2)); // on the island
+    assertEquals(Relation.CELL_INSIDE_QUERY, Polygon.relate(polygons, -2, 2, -2, 2)); // on the island
-    assertFalse(Polygon.contains(polygons, 6, 7, 6, 7)); // in the hole
+    assertEquals(Relation.CELL_OUTSIDE_QUERY, Polygon.relate(polygons, 6, 7, 6, 7)); // in the hole
-    assertTrue(Polygon.contains(polygons, 24, 25, 24, 25)); // on the mainland
+    assertEquals(Relation.CELL_INSIDE_QUERY, Polygon.relate(polygons, 24, 25, 24, 25)); // on the mainland
-    assertFalse(Polygon.contains(polygons, 51, 52, 51, 52)); // in the ocean
+    assertEquals(Relation.CELL_OUTSIDE_QUERY, Polygon.relate(polygons, 51, 52, 51, 52)); // in the ocean
-    assertFalse(Polygon.contains(polygons, -60, 60, -60, 60)); // enclosing us completely
+    assertEquals(Relation.CELL_CROSSES_QUERY, Polygon.relate(polygons, -60, 60, -60, 60)); // enclosing us completely
-    assertFalse(Polygon.contains(polygons, 49, 51, 49, 51)); // overlapping the mainland
+    assertEquals(Relation.CELL_CROSSES_QUERY, Polygon.relate(polygons, 49, 51, 49, 51)); // overlapping the mainland
-    assertFalse(Polygon.contains(polygons, 9, 11, 9, 11)); // overlapping the hole
+    assertEquals(Relation.CELL_CROSSES_QUERY, Polygon.relate(polygons, 9, 11, 9, 11)); // overlapping the hole
-    assertFalse(Polygon.contains(polygons, 5, 6, 5, 6)); // overlapping the island
+    assertEquals(Relation.CELL_CROSSES_QUERY, Polygon.relate(polygons, 5, 6, 5, 6)); // overlapping the island
-
-    // crosses(box): this can conservatively return true
-    assertTrue(Polygon.crosses(polygons, -60, 60, -60, 60)); // enclosing us completely
-    assertTrue(Polygon.crosses(polygons, 49, 51, 49, 51)); // overlapping the mainland and ocean
-    assertTrue(Polygon.crosses(polygons, 9, 11, 9, 11)); // overlapping the hole and mainland
-    assertTrue(Polygon.crosses(polygons, 5, 6, 5, 6)); // overlapping the island
   }
   
   public void testPacMan() throws Exception {
@@ -110,8 +105,7 @@ public class TestPolygon extends LuceneTestCase {
 
     // test cell crossing poly
     Polygon polygon = new Polygon(py, px);
-    assertTrue(polygon.crosses(yMin, yMax, xMin, xMax));
+    assertEquals(Relation.CELL_CROSSES_QUERY, polygon.relate(yMin, yMax, xMin, xMax));
-    assertFalse(polygon.contains(yMin, yMax, xMin, xMax));
   }
   
   public void testBoundingBox() throws Exception {
@@ -154,7 +148,7 @@ public class TestPolygon extends LuceneTestCase {
       for (int j = 0; j < 100; j++) {
         Rectangle rectangle = GeoTestUtil.nextSimpleBox();
         // allowed to conservatively return false
-        if (polygon.contains(rectangle.minLat, rectangle.maxLat, rectangle.minLon, rectangle.maxLon)) {
+        if (polygon.relate(rectangle.minLat, rectangle.maxLat, rectangle.minLon, rectangle.maxLon) == Relation.CELL_INSIDE_QUERY) {
           for (int k = 0; k < 1000; k++) {
             // this tests in our range but sometimes outside! so we have to double-check its really in other box
             double latitude = nextLatitudeAround(rectangle.minLat, rectangle.maxLat);
@@ -183,7 +177,7 @@ public class TestPolygon extends LuceneTestCase {
         for (int j = 0; j < 10; j++) {
           Rectangle rectangle = GeoTestUtil.nextSimpleBoxNear(polyLats[vertex], polyLons[vertex]);
           // allowed to conservatively return false
-          if (polygon.contains(rectangle.minLat, rectangle.maxLat, rectangle.minLon, rectangle.maxLon)) {
+          if (polygon.relate(rectangle.minLat, rectangle.maxLat, rectangle.minLon, rectangle.maxLon) == Relation.CELL_INSIDE_QUERY) {
             for (int k = 0; k < 100; k++) {
               // this tests in our range but sometimes outside! so we have to double-check its really in other box
               double latitude = nextLatitudeAround(rectangle.minLat, rectangle.maxLat);
@@ -207,8 +201,7 @@ public class TestPolygon extends LuceneTestCase {
       for (int j = 0; j < 100; j++) {
         Rectangle rectangle = GeoTestUtil.nextSimpleBox();
         // allowed to conservatively return true.
-        if (polygon.contains(rectangle.minLat, rectangle.maxLat, rectangle.minLon, rectangle.maxLon) == false &&
+        if (polygon.relate(rectangle.minLat, rectangle.maxLat, rectangle.minLon, rectangle.maxLon) == Relation.CELL_OUTSIDE_QUERY) {
-            polygon.crosses(rectangle.minLat, rectangle.maxLat, rectangle.minLon, rectangle.maxLon) == false) {
           for (int k = 0; k < 1000; k++) {
             // this tests in our range but sometimes outside! so we have to double-check its really in other box
             double latitude = nextLatitudeAround(rectangle.minLat, rectangle.maxLat);
@@ -237,8 +230,7 @@ public class TestPolygon extends LuceneTestCase {
         for (int j = 0; j < 10; j++) {
           Rectangle rectangle = GeoTestUtil.nextSimpleBoxNear(polyLats[vertex], polyLons[vertex]);
           // allowed to conservatively return true.
-          if (polygon.contains(rectangle.minLat, rectangle.maxLat, rectangle.minLon, rectangle.maxLon) == false &&
+          if (polygon.relate(rectangle.minLat, rectangle.maxLat, rectangle.minLon, rectangle.maxLon) == Relation.CELL_OUTSIDE_QUERY) {
-              polygon.crosses(rectangle.minLat, rectangle.maxLat, rectangle.minLon, rectangle.maxLon) == false) {
             for (int k = 0; k < 100; k++) {
               // this tests in our range but sometimes outside! so we have to double-check its really in other box
               double latitude = nextLatitudeAround(rectangle.minLat, rectangle.maxLat);

lucene/sandbox/src/java/org/apache/lucene/document/LatLonGrid.java

@@ -0,0 +1,155 @@
+/*
+ * 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.
+ */
+package org.apache.lucene.document;
+
+import org.apache.lucene.geo.Polygon;
+import org.apache.lucene.index.PointValues.Relation;
+import org.apache.lucene.util.FixedBitSet;
+
+/**
+ * This is a temporary hack, until some polygon methods have better performance!
+ * <p>
+ * When this file is removed then we have made good progress! In general we don't call
+ * the point-in-polygon algorithm that much, because of how BKD divides up the data. But
+ * today the method is very slow (general to all polygons, linear with the number of vertices).
+ * At the same time polygon-rectangle relation operations are also slow in the same way, this
+ * just really ensures they are the bottleneck by removing most of the point-in-polygon calls.
+ * <p>
+ * See the "grid" algorithm description here: http://erich.realtimerendering.com/ptinpoly/
+ * A few differences:
+ * <ul>
+ *   <li> We work in an integer encoding, so edge cases are simpler.
+ *   <li> We classify each grid cell as "contained", "not contained", or "don't know".
+ *   <li> We form a grid over a potentially complex multipolygon with holes.
+ *   <li> Construction is less efficient because we do not do anything "smart" such
+ *        as following polygon edges. 
+ *   <li> Instead we construct a baby tree to reduce the number of relation operations,
+ *        which are currently expensive.
+ * </ul>
+ */
+// TODO: just make a more proper tree (maybe in-ram BKD)? then we can answer most 
+// relational operations as rectangle <-> rectangle relations in integer space in log(n) time..
+final class LatLonGrid {
+  // must be a power of two!
+  static final int GRID_SIZE = 1<<5;
+  final int minLat;
+  final int maxLat;
+  final int minLon;
+  final int maxLon;
+  // TODO: something more efficient than parallel bitsets? maybe one bitset?
+  final FixedBitSet haveAnswer = new FixedBitSet(GRID_SIZE * GRID_SIZE);
+  final FixedBitSet answer = new FixedBitSet(GRID_SIZE * GRID_SIZE);
+  
+  final long latPerCell;
+  final long lonPerCell;
+  
+  final Polygon[] polygons;
+  
+  LatLonGrid(int minLat, int maxLat, int minLon, int maxLon, Polygon... polygons) {
+    this.minLat = minLat;
+    this.maxLat = maxLat;
+    this.minLon = minLon;
+    this.maxLon = maxLon;
+    this.polygons = polygons;
+    if (minLon > maxLon) {
+      // maybe make 2 grids if you want this? 
+      throw new IllegalArgumentException("Grid cannot cross the dateline");
+    }
+    if (minLat > maxLat) {
+      throw new IllegalArgumentException("bogus grid");
+    }
+    long latitudeRange = maxLat - (long) minLat;
+    long longitudeRange = maxLon - (long) minLon;
+    // we spill over the edge of the bounding box in each direction a bit,
+    // but it prevents edge case bugs.
+    latPerCell = latitudeRange / (GRID_SIZE - 1);
+    lonPerCell = longitudeRange / (GRID_SIZE - 1);
+    fill(polygons, 0, GRID_SIZE, 0, GRID_SIZE);
+  }
+  
+  /** fills a 2D range of grid cells [minLatIndex .. maxLatIndex) X [minLonIndex .. maxLonIndex) */
+  void fill(Polygon[] polygons, int minLatIndex, int maxLatIndex, int minLonIndex, int maxLonIndex) {
+    // grid cells at the edge of the bounding box are typically smaller than normal, because we spill over.
+    long cellMinLat = minLat + (minLatIndex * latPerCell);
+    long cellMaxLat = Math.min(maxLat, minLat + (maxLatIndex * latPerCell) - 1);
+    long cellMinLon = minLon + (minLonIndex * lonPerCell);
+    long cellMaxLon = Math.min(maxLon, minLon + (maxLonIndex * lonPerCell) - 1);
+
+    assert cellMinLat <= maxLat && cellMinLon <= maxLon;
+    assert cellMaxLat >= cellMinLat;
+    assert cellMaxLon >= cellMinLon;
+
+    Relation relation = Polygon.relate(polygons, LatLonPoint.decodeLatitude((int) cellMinLat), 
+                                                 LatLonPoint.decodeLatitude((int) cellMaxLat), 
+                                                 LatLonPoint.decodeLongitude((int) cellMinLon), 
+                                                 LatLonPoint.decodeLongitude((int) cellMaxLon));
+    if (relation != Relation.CELL_CROSSES_QUERY) {
+      // we know the answer for this region, fill the cell range
+      for (int i = minLatIndex; i < maxLatIndex; i++) {
+        for (int j = minLonIndex; j < maxLonIndex; j++) {
+          int index = i * GRID_SIZE + j;
+          assert haveAnswer.get(index) == false;
+          haveAnswer.set(index);
+          if (relation == Relation.CELL_INSIDE_QUERY) {
+            answer.set(index);
+          }
+        }
+      }
+    } else if (minLatIndex == maxLatIndex - 1) {
+      // nothing more to do: this is a single grid cell (leaf node) and
+      // is an edge case for the polygon.
+    } else {
+      // grid range crosses our polygon, keep recursing.
+      int midLatIndex = (minLatIndex + maxLatIndex) >>> 1;
+      int midLonIndex = (minLonIndex + maxLonIndex) >>> 1;
+      fill(polygons, minLatIndex, midLatIndex, minLonIndex, midLonIndex);
+      fill(polygons, minLatIndex, midLatIndex, midLonIndex, maxLonIndex);
+      fill(polygons, midLatIndex, maxLatIndex, minLonIndex, midLonIndex);
+      fill(polygons, midLatIndex, maxLatIndex, midLonIndex, maxLonIndex);
+    }
+  }
+  
+  /** Returns true if inside one of our polygons, false otherwise */
+  boolean contains(int latitude, int longitude) {
+    // first see if the grid knows the answer
+    int index = index(latitude, longitude);
+    if (index == -1) {
+      return false; // outside of bounding box range
+    } else if (haveAnswer.get(index)) {
+      return answer.get(index);
+    }
+
+    // the grid is unsure (boundary): do a real test.
+    double docLatitude = LatLonPoint.decodeLatitude(latitude);
+    double docLongitude = LatLonPoint.decodeLongitude(longitude);
+    return Polygon.contains(polygons, docLatitude, docLongitude);
+  }
+  
+  /** Returns grid index of lat/lon, or -1 if the value is outside of the bounding box. */
+  private int index(int latitude, int longitude) {
+    if (latitude < minLat || latitude > maxLat || longitude < minLon || longitude > maxLon) {
+      return -1; // outside of bounding box range
+    }
+    
+    long latRel = latitude - (long) minLat;
+    long lonRel = longitude - (long) minLon;
+    
+    int latIndex = (int) (latRel / latPerCell);
+    int lonIndex = (int) (lonRel / lonPerCell);
+    return latIndex * GRID_SIZE + lonIndex;
+  }
+}

lucene/sandbox/src/java/org/apache/lucene/document/LatLonPointInPolygonQuery.java

@@ -100,6 +100,11 @@ final class LatLonPointInPolygonQuery extends Query {
     }
     final float matchCost = cumulativeCost;
 
+    final LatLonGrid grid = new LatLonGrid(LatLonPoint.encodeLatitude(box.minLat), 
+                                           LatLonPoint.encodeLatitude(box.maxLat),
+                                           LatLonPoint.encodeLongitude(box.minLon),
+                                           LatLonPoint.encodeLongitude(box.maxLon), polygons);
+
     return new ConstantScoreWeight(this) {
 
       @Override
@@ -164,13 +169,7 @@ final class LatLonPointInPolygonQuery extends Query {
                              double cellMaxLat = LatLonPoint.decodeLatitude(maxPackedValue, 0);
                              double cellMaxLon = LatLonPoint.decodeLongitude(maxPackedValue, Integer.BYTES);
 
-                             if (Polygon.contains(polygons, cellMinLat, cellMaxLat, cellMinLon, cellMaxLon)) {
+                             return Polygon.relate(polygons, cellMinLat, cellMaxLat, cellMinLon, cellMaxLon);
-                               return Relation.CELL_INSIDE_QUERY;
-                             } else if (Polygon.crosses(polygons, cellMinLat, cellMaxLat, cellMinLon, cellMaxLon)) {
-                               return Relation.CELL_CROSSES_QUERY;
-                             } else {
-                               return Relation.CELL_OUTSIDE_QUERY;
-                             }
                            }
                          });
 
@@ -194,9 +193,9 @@ final class LatLonPointInPolygonQuery extends Query {
               int count = docValues.count();
               for (int i = 0; i < count; i++) {
                 long encoded = docValues.valueAt(i);
-                double docLatitude = LatLonPoint.decodeLatitude((int)(encoded >> 32));
+                int latitudeBits = (int)(encoded >> 32);
-                double docLongitude = LatLonPoint.decodeLongitude((int)(encoded & 0xFFFFFFFF));
+                int longitudeBits = (int)(encoded & 0xFFFFFFFF);
-                if (Polygon.contains(polygons, docLatitude, docLongitude)) {
+                if (grid.contains(latitudeBits, longitudeBits)) {
                   return true;
                 }
               }

lucene/sandbox/src/test/org/apache/lucene/document/TestLatLonGrid.java

@@ -0,0 +1,50 @@
+/*
+ * 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.
+ */
+package org.apache.lucene.document;
+
+import org.apache.lucene.geo.GeoTestUtil;
+import org.apache.lucene.geo.Polygon;
+import org.apache.lucene.geo.Rectangle;
+import org.apache.lucene.util.LuceneTestCase;
+import org.apache.lucene.util.TestUtil;
+
+/** tests against LatLonGrid (avoiding indexing/queries) */
+public class TestLatLonGrid extends LuceneTestCase {
+
+  /** If the grid returns true, then any point in that cell should return true as well */
+  public void testRandom() throws Exception {
+    for (int i = 0; i < 100; i++) {
+      Polygon polygon = GeoTestUtil.nextPolygon();
+      Rectangle box = Rectangle.fromPolygon(new Polygon[] { polygon });
+      int minLat = LatLonPoint.encodeLatitude(box.minLat);
+      int maxLat = LatLonPoint.encodeLatitude(box.maxLat);
+      int minLon = LatLonPoint.encodeLongitude(box.minLon);
+      int maxLon = LatLonPoint.encodeLongitude(box.maxLon);
+      LatLonGrid grid = new LatLonGrid(minLat, maxLat, minLon, maxLon, polygon);
+      // we are in integer space... but exhaustive testing is slow!
+      for (int j = 0; j < 10000; j++) {
+        int lat = TestUtil.nextInt(random(), minLat, maxLat);
+        int lon = TestUtil.nextInt(random(), minLon, maxLon);
+
+        boolean expected = polygon.contains(LatLonPoint.decodeLatitude(lat), 
+                                            LatLonPoint.decodeLongitude(lon));
+        boolean actual = grid.contains(lat, lon);
+        assertEquals(expected, actual);
+      }
+    }
+  }
+}

lucene/spatial/src/java/org/apache/lucene/spatial/geopoint/search/GeoPointInPolygonQueryImpl.java

@@ -21,6 +21,7 @@ import java.util.Objects;
 import org.apache.lucene.search.MultiTermQuery;
 import org.apache.lucene.spatial.geopoint.document.GeoPointField.TermEncoding;
 import org.apache.lucene.geo.Polygon;
+import org.apache.lucene.index.PointValues.Relation;
 
 /** Package private implementation for the public facing GeoPointInPolygonQuery delegate class.
  *
@@ -58,18 +59,17 @@ final class GeoPointInPolygonQueryImpl extends GeoPointInBBoxQueryImpl {
 
     @Override
     protected boolean cellCrosses(final double minLat, final double maxLat, final double minLon, final double maxLon) {
-      return Polygon.crosses(polygons, minLat, maxLat, minLon, maxLon);
+      return Polygon.relate(polygons, minLat, maxLat, minLon, maxLon) == Relation.CELL_CROSSES_QUERY;
     }
 
     @Override
     protected boolean cellWithin(final double minLat, final double maxLat, final double minLon, final double maxLon) {
-      return Polygon.contains(polygons, minLat, maxLat, minLon, maxLon);
+      return Polygon.relate(polygons, minLat, maxLat, minLon, maxLon) == Relation.CELL_INSIDE_QUERY;
     }
 
     @Override
     protected boolean cellIntersectsShape(final double minLat, final double maxLat, final double minLon, final double maxLon) {
-      return cellContains(minLat, maxLat, minLon, maxLon) || cellWithin(minLat, maxLat, minLon, maxLon)
+      return Polygon.relate(polygons, minLat, maxLat, minLon, maxLon) != Relation.CELL_OUTSIDE_QUERY;
-        || cellCrosses(minLat, maxLat, minLon, maxLon);
     }
 
     /**