LUCENE-7195: Clockwise/counterclockwise detection was rotating coordinates in the wrong direction.

lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoPolygonFactory.java

@@ -116,6 +116,7 @@ public class GeoPolygonFactory {
       final Boolean isPoleInside = isInsidePolygon(pole, pointList);
       if (isPoleInside != null) {
         // Legal pole
+        //System.out.println("Pole = "+pole+"; isInside="+isPoleInside+"; pointList = "+pointList);
         return makeGeoPolygon(planetModel, pointList, holes, pole, isPoleInside);
       }
       // If pole choice was illegal, try another one
@@ -177,24 +178,17 @@ public class GeoPolygonFactory {
    */
   private static Boolean isInsidePolygon(final GeoPoint point, final List<GeoPoint> polyPoints) {
     // First, compute sine and cosine of pole point latitude and longitude
-    final double norm = 1.0 / point.magnitude();
-    final double xyDenom = Math.sqrt(point.x * point.x + point.y * point.y);
-    final double sinLatitude = point.z * norm;
-    final double cosLatitude = xyDenom * norm;
-    final double sinLongitude;
-    final double cosLongitude;
-    if (Math.abs(xyDenom) < Vector.MINIMUM_RESOLUTION) {
-      sinLongitude = 0.0;
-      cosLongitude = 1.0;
-    } else {
-      final double xyNorm = 1.0 / xyDenom;
-      sinLongitude = point.y * xyNorm;
-      cosLongitude = point.x * xyNorm;
-    }
+    final double latitude = point.getLatitude();
+    final double longitude = point.getLongitude();
+    final double sinLatitude = Math.sin(latitude);
+    final double cosLatitude = Math.cos(latitude);
+    final double sinLongitude = Math.sin(longitude);
+    final double cosLongitude = Math.cos(longitude);
     
     // Now, compute the incremental arc distance around the points of the polygon
     double arcDistance = 0.0;
     Double prevAngle = null;
+    //System.out.println("Computing angles:");
     for (final GeoPoint polyPoint : polyPoints) {
       final Double angle = computeAngle(polyPoint, sinLatitude, cosLatitude, sinLongitude, cosLongitude);
       if (angle == null) {
@@ -215,6 +209,7 @@ public class GeoPolygonFactory {
         }
         //System.out.println(" angle delta = "+angleDelta);
         arcDistance += angleDelta;
+        //System.out.println(" For point "+polyPoint+" angle is "+angle+"; delta is "+angleDelta+"; arcDistance is "+arcDistance);
       }
       prevAngle = angle;
     }
@@ -237,7 +232,9 @@ public class GeoPolygonFactory {
       }
       //System.out.println(" angle delta = "+angleDelta);
       arcDistance += angleDelta;
+      //System.out.println(" For point "+polyPoints.get(0)+" angle is "+lastAngle+"; delta is "+angleDelta+"; arcDistance is "+arcDistance);
     }
+
     // Clockwise == inside == negative
     //System.out.println("Arcdistance = "+arcDistance);
     if (Math.abs(arcDistance) < Vector.MINIMUM_RESOLUTION) {
@@ -266,21 +263,23 @@ public class GeoPolygonFactory {
     // We need to rotate the point in question into the coordinate frame specified by
     // the lat and lon trig functions.
     // To do this we need to do two rotations on it.  First rotation is in x/y.  Second rotation is in x/z.
+    // And we rotate in the negative direction.
     // So:
-    // x1 = x0 cos az - y0 sin az
-    // y1 = x0 sin az + y0 cos az
+    // x1 = x0 cos az + y0 sin az
+    // y1 = - x0 sin az + y0 cos az
     // z1 = z0
-    // x2 = x1 cos al - z1 sin al
+    // x2 = x1 cos al + z1 sin al
     // y2 = y1
-    // z2 = x1 sin al + z1 cos al
+    // z2 = - x1 sin al + z1 cos al
       
-    final double x1 = point.x * cosLongitude - point.y * sinLongitude;
-    final double y1 = point.x * sinLongitude + point.y * cosLongitude;
+    final double x1 = point.x * cosLongitude + point.y * sinLongitude;
+    final double y1 = - point.x * sinLongitude + point.y * cosLongitude;
     final double z1 = point.z;
-    //final double x2 = x1 * cosLatitude - z1 * sinLatitude;
-    final double y2 = y1;
-    final double z2 = x1 * sinLatitude + z1 * cosLatitude;
       
+    // final double x2 = x1 * cosLatitude + z1 * sinLatitude;
+    final double y2 = y1;
+    final double z2 = - x1 * sinLatitude + z1 * cosLatitude;
+    
     // Now we should be looking down the X axis; the original point has rotated coordinates (N, 0, 0).
     // So we can just compute the angle using y2 and z2.  (If Math.sqrt(y2*y2 + z2 * z2) is 0.0, then the point is on the pole and we need another one).
     if (Math.sqrt(y2*y2 + z2*z2) < Vector.MINIMUM_RESOLUTION) {

lucene/spatial3d/src/test/org/apache/lucene/spatial3d/TestGeo3DPoint.java

@@ -517,6 +517,16 @@ public class TestGeo3DPoint extends LuceneTestCase {
     verify(lats, lons);
   }
 
+  public void testPolygonOrdering() {
+    final double[] lats = new double[] {
+      51.204382859999996, 50.89947531437482, 50.8093624806861,50.8093624806861, 50.89947531437482, 51.204382859999996, 51.51015366140113, 51.59953838204167, 51.59953838204167, 51.51015366140113, 51.204382859999996};
+    final double[] lons = new double[] {
+      0.8747711978759765, 0.6509219832137298, 0.35960265165247807, 0.10290284834752167, -0.18841648321373008, -0.41226569787597667, -0.18960465285650027, 0.10285893781346236, 0.35964656218653757, 0.6521101528565002, 0.8747711978759765};
+    final Query q = Geo3DPoint.newPolygonQuery("point", new Polygon(lats, lons));
+    //System.out.println(q);
+    assertTrue(!q.toString().contains("GeoConcavePolygon"));
+  }
+  
   private static final double MEAN_EARTH_RADIUS_METERS = PlanetModel.WGS84_MEAN;
   
   private static Query random3DQuery(final String field) {
@@ -982,26 +992,26 @@ public class TestGeo3DPoint extends LuceneTestCase {
     // x1 = x0 cos T - y0 sin T
     // y1 = x0 sin T + y0 cos T
     // We're in essence undoing the following transformation (from GeoPolygonFactory):
-    // x1 = x0 cos az - y0 sin az
-    // y1 = x0 sin az + y0 cos az
+    // x1 = x0 cos az + y0 sin az
+    // y1 = - x0 sin az + y0 cos az
     // z1 = z0
-    // x2 = x1 cos al - z1 sin al
+    // x2 = x1 cos al + z1 sin al
     // y2 = y1
-    // z2 = x1 sin al + z1 cos al
+    // z2 = - x1 sin al + z1 cos al
     // So, we reverse the order of the transformations, AND we transform backwards.
     // Transforming backwards means using these identities: sin(-angle) = -sin(angle), cos(-angle) = cos(angle)
     // So:
-    // x1 = x0 cos al + z0 sin al
+    // x1 = x0 cos al - z0 sin al
     // y1 = y0
-    // z1 = - x0 sin al + z0 cos al
-    // x2 = x1 cos az + y1 sin az
-    // y2 = - x1 sin az + y1 cos az
+    // z1 = x0 sin al + z0 cos al
+    // x2 = x1 cos az - y1 sin az
+    // y2 = x1 sin az + y1 cos az
     // z2 = z1
-    final double x1 = x * cosLatitude + z * sinLatitude;
+    final double x1 = x * cosLatitude - z * sinLatitude;
     final double y1 = y;
-    final double z1 = - x * sinLatitude + z * cosLatitude;
-    final double x2 = x1 * cosLongitude + y1 * sinLongitude;
-    final double y2 = - x1 * sinLongitude + y1 * cosLongitude;
+    final double z1 = x * sinLatitude + z * cosLatitude;
+    final double x2 = x1 * cosLongitude - y1 * sinLongitude;
+    final double y2 = x1 * sinLongitude + y1 * cosLongitude;
     final double z2 = z1;
 
     // Scale final (x,y,z) to land on planet surface