You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@lucene.apache.org by sa...@apache.org on 2016/06/16 14:19:41 UTC
[2/2] lucene-solr:branch_5x: LUCENE-7139: fix bugs in geo3d's
vincenty distance implementation
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 <mi...@apache.org>
Authored: Fri Mar 25 07:59:56 2016 -0400
Committer: Steve Rowe <sa...@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