lucene-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From sar...@apache.org
Subject [2/2] lucene-solr:branch_5x: LUCENE-7139: fix bugs in geo3d's vincenty distance implementation
Date Thu, 16 Jun 2016 14:19:41 GMT
LUCENE-7139: fix bugs in geo3d's vincenty distance implementation


Project: http://git-wip-us.apache.org/repos/asf/lucene-solr/repo
Commit: http://git-wip-us.apache.org/repos/asf/lucene-solr/commit/3691c73a
Tree: http://git-wip-us.apache.org/repos/asf/lucene-solr/tree/3691c73a
Diff: http://git-wip-us.apache.org/repos/asf/lucene-solr/diff/3691c73a

Branch: refs/heads/branch_5x
Commit: 3691c73ab767571ebe835fa04f22c3c2eed97c8a
Parents: c541d59
Author: Mike McCandless <mikemccand@apache.org>
Authored: Fri Mar 25 07:59:56 2016 -0400
Committer: Steve Rowe <sarowe@apache.org>
Committed: Thu Jun 16 10:19:27 2016 -0400

----------------------------------------------------------------------
 .../spatial/spatial4j/geo3d/GeoPointTest.java   | 13 +++
 .../org/apache/lucene/geo3d/PlanetModel.java    | 96 +++++++++++---------
 2 files changed, 64 insertions(+), 45 deletions(-)
----------------------------------------------------------------------


http://git-wip-us.apache.org/repos/asf/lucene-solr/blob/3691c73a/lucene/spatial/src/test/org/apache/lucene/spatial/spatial4j/geo3d/GeoPointTest.java
----------------------------------------------------------------------
diff --git a/lucene/spatial/src/test/org/apache/lucene/spatial/spatial4j/geo3d/GeoPointTest.java
b/lucene/spatial/src/test/org/apache/lucene/spatial/spatial4j/geo3d/GeoPointTest.java
index e652581..3587c5e 100644
--- a/lucene/spatial/src/test/org/apache/lucene/spatial/spatial4j/geo3d/GeoPointTest.java
+++ b/lucene/spatial/src/test/org/apache/lucene/spatial/spatial4j/geo3d/GeoPointTest.java
@@ -70,6 +70,19 @@ public class GeoPointTest extends LuceneTestCase {
       // Compute ellipsoid distance; it should agree for a sphere
       final double surfaceDistance = PlanetModel.SPHERE.surfaceDistance(p1,p2);
       assertEquals(arcDistance, surfaceDistance, 1e-6);
+      // Now try some WGS84 points (taken randomly and compared against a known-good implementation)
+      assertEquals(1.1444648695765323, PlanetModel.WGS84.surfaceDistance(
+        new GeoPoint(PlanetModel.WGS84, 0.038203808753702884, -0.6701260455506466),
+        new GeoPoint(PlanetModel.WGS84, -0.8453720422675458, 0.1737353153814496)), 1e-6);
+      assertEquals(1.4345148695890722, PlanetModel.WGS84.surfaceDistance(
+        new GeoPoint(PlanetModel.WGS84, 0.5220926323378574, 0.6758041581907408),
+        new GeoPoint(PlanetModel.WGS84, -0.8453720422675458, 0.1737353153814496)), 1e-6);
+      assertEquals(2.32418144616446, PlanetModel.WGS84.surfaceDistance(
+        new GeoPoint(PlanetModel.WGS84, 0.09541335760967473, 1.2091829760623236),
+        new GeoPoint(PlanetModel.WGS84, -0.8501591797459979, -2.3044806381627594)), 1e-6);
+      assertEquals(2.018421047005435, PlanetModel.WGS84.surfaceDistance(
+        new GeoPoint(PlanetModel.WGS84, 0.3402853531962009, -0.43544195327249957),
+        new GeoPoint(PlanetModel.WGS84, -0.8501591797459979, -2.3044806381627594)), 1e-6);
     }
   }
 

http://git-wip-us.apache.org/repos/asf/lucene-solr/blob/3691c73a/lucene/spatial3d/src/java/org/apache/lucene/geo3d/PlanetModel.java
----------------------------------------------------------------------
diff --git a/lucene/spatial3d/src/java/org/apache/lucene/geo3d/PlanetModel.java b/lucene/spatial3d/src/java/org/apache/lucene/geo3d/PlanetModel.java
index b113e85..572c08e 100644
--- a/lucene/spatial3d/src/java/org/apache/lucene/geo3d/PlanetModel.java
+++ b/lucene/spatial3d/src/java/org/apache/lucene/geo3d/PlanetModel.java
@@ -189,64 +189,70 @@ public class PlanetModel {
   }
   
   /** Compute surface distance between two points.
-   * @param p1 is the first point.
-   * @param p2 is the second point.
+   * @param pt1 is the first point.
+   * @param pt2 is the second point.
    * @return the adjusted angle, when multiplied by the mean earth radius, yields a surface
distance.  This will differ
    * from GeoPoint.arcDistance() only when the planet model is not a sphere. @see {@link
org.apache.lucene.geo3d.GeoPoint#arcDistance(GeoPoint)}
    */
-  public double surfaceDistance(final GeoPoint p1, final GeoPoint p2) {
-    final double latA = p1.getLatitude();
-    final double lonA = p1.getLongitude();
-    final double latB = p2.getLatitude();
-    final double lonB = p2.getLongitude();
+  public double surfaceDistance(final GeoPoint pt1, final GeoPoint pt2) {
+    final double L = pt2.getLongitude() - pt1.getLongitude();
+    final double U1 = Math.atan((1.0-flattening) * Math.tan(pt1.getLatitude()));
+    final double U2 = Math.atan((1.0-flattening) * Math.tan(pt2.getLatitude()));
+
+    final double sinU1 = Math.sin(U1);
+    final double cosU1 = Math.cos(U1);
+    final double sinU2 = Math.sin(U2);
+    final double cosU2 = Math.cos(U2);
+
+    final double dCosU1CosU2 = cosU1 * cosU2;
+    final double dCosU1SinU2 = cosU1 * sinU2;
+
+    final double dSinU1SinU2 = sinU1 * sinU2;
+    final double dSinU1CosU2 = sinU1 * cosU2;
 
-    final double L = lonB - lonA;
-    final double oF = 1.0 - this.flattening;
-    final double U1 = Math.atan(oF * Math.tan(latA));
-    final double U2 = Math.atan(oF * Math.tan(latB));
-    final double sU1 = Math.sin(U1);
-    final double cU1 = Math.cos(U1);
-    final double sU2 = Math.sin(U2);
-    final double cU2 = Math.cos(U2);
 
-    double sigma, sinSigma, cosSigma;
-    double cos2Alpha, cos2SigmaM;
-    
     double lambda = L;
-    double iters = 100;
+    double lambdaP = Math.PI * 2.0;
+    int iterLimit = 0;
+    double cosSqAlpha;
+    double sinSigma;
+    double cos2SigmaM;
+    double cosSigma;
+    double sigma;
+    double sinAlpha;
+    double C;
+    double sinLambda, cosLambda;
 
     do {
-      final double sinLambda = Math.sin(lambda);
-      final double cosLambda = Math.cos(lambda);
-      sinSigma = Math.sqrt((cU2 * sinLambda) * (cU2 * sinLambda) + (cU1 * sU2 - sU1 * cU2
* cosLambda)
-          * (cU1 * sU2 - sU1 * cU2 * cosLambda));
-      if (Math.abs(sinSigma) < Vector.MINIMUM_RESOLUTION)
-        return 0.0;
+      sinLambda = Math.sin(lambda);
+      cosLambda = Math.cos(lambda);
+      sinSigma = Math.sqrt((cosU2*sinLambda) * (cosU2*sinLambda) +
+                                    (dCosU1SinU2 - dSinU1CosU2 * cosLambda) * (dCosU1SinU2
- dSinU1CosU2 * cosLambda));
 
-      cosSigma = sU1 * sU2 + cU1 * cU2 * cosLambda;
+      if (sinSigma==0.0) {
+        return 0.0;
+      }
+      cosSigma = dSinU1SinU2 + dCosU1CosU2 * cosLambda;
       sigma = Math.atan2(sinSigma, cosSigma);
-      final double sinAlpha = cU1 * cU2 * sinLambda / sinSigma;
-      cos2Alpha = 1.0 - sinAlpha * sinAlpha;
-      cos2SigmaM = cosSigma - 2.0 * sU1 * sU2 / cos2Alpha;
-
-      final double c = this.flattening * 0.625 * cos2Alpha * (4.0 + this.flattening * (4.0
- 3.0 * cos2Alpha));
-      final double lambdaP = lambda;
-      lambda = L + (1.0 - c) * this.flattening * sinAlpha * (sigma + c * sinSigma * (cos2SigmaM
+ c * cosSigma *
-          (-1.0 + 2.0 * cos2SigmaM * cos2SigmaM)));
-      if (Math.abs(lambda - lambdaP) < Vector.MINIMUM_RESOLUTION)
-        break;
-    } while (--iters > 0);
+      sinAlpha = dCosU1CosU2 * sinLambda / sinSigma;
+      cosSqAlpha = 1.0 - sinAlpha * sinAlpha;
+      cos2SigmaM = cosSigma - 2.0 * dSinU1SinU2 / cosSqAlpha;
 
-    if (iters == 0)
-      return 0.0;
+      if (Double.isNaN(cos2SigmaM))
+        cos2SigmaM = 0.0;  // equatorial line: cosSqAlpha=0
+      C = flattening / 16.0 * cosSqAlpha * (4.0 + flattening * (4.0 - 3.0 * cosSqAlpha));
+      lambdaP = lambda;
+      lambda = L + (1.0 - C) * flattening * sinAlpha *
+        (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1.0 + 2.0 * cos2SigmaM *cos2SigmaM)));
+    } while (Math.abs(lambda-lambdaP) > Vector.MINIMUM_RESOLUTION && ++iterLimit
< 40);
 
-    final double uSq = cos2Alpha * this.squareRatio;
-    final double A = 1.0 + uSq * 0.00006103515625 * (4096.0 + uSq * (-768.0 + uSq * (320.0
- 175.0 * uSq)));
-    final double B = uSq * 0.0009765625 * (256.0 + uSq * (-128.0 + uSq * (74.0 - 47.0 * uSq)));
-    final double deltaSigma = B * sinSigma * (cos2SigmaM + B * 0.25 * (cosSigma * (-1.0 +
2.0 * cos2SigmaM * cos2SigmaM) - B * 0.16666666666666666666667 * cos2SigmaM
-            * (-3.0 + 4.0 * sinSigma * sinSigma) * (-3.0 + 4.0 * cos2SigmaM * cos2SigmaM)));
+    final double uSq = cosSqAlpha * this.squareRatio;
+    final double A = 1.0 + uSq / 16384.0 * (4096.0 + uSq * (-768.0 + uSq * (320.0 - 175.0
* uSq)));
+    final double B = uSq / 1024.0 * (256.0 + uSq * (-128.0 + uSq * (74.0 - 47.0 * uSq)));
+    final double deltaSigma = B * sinSigma * (cos2SigmaM + B / 4.0 * (cosSigma * (-1.0 +
2.0 * cos2SigmaM * cos2SigmaM)-
+                                        B / 6.0 * cos2SigmaM * (-3.0 + 4.0 * sinSigma * sinSigma)
* (-3.0 + 4.0 * cos2SigmaM * cos2SigmaM)));
 
-    return this.c * A * (sigma - deltaSigma);
+    return c * A * (sigma - deltaSigma);
   }
 
   @Override


Mime
View raw message