Apache
/
lucene
mirror of https://github.com/apache/lucene.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

LUCENE-7139: fix bugs in geo3d's vincenty distance implementation

This commit is contained in:
Mike McCandless 2016-03-25 07:59:56 -04:00
parent 90958799c6
commit 9726903c99
3 changed files with 67 additions and 45 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
lucene/CHANGES.txt
View File

@ -211,6 +211,9 @@ Bug Fixes
* LUCENE-7111: DocValuesRangeQuery.newLongRange behaves incorrectly for
Long.MAX_VALUE and Long.MIN_VALUE (Ishan Chattopadhyaya via Steve Rowe)
* LUCENE-7139: Fix bugs in geo3d's Vincenty surface distance
implementation (Karl Wright via Mike McCandless)
Other
* LUCENE-7035: Upgrade icu4j to 56.1/unicode 8. (Robert Muir)

96
lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/PlanetModel.java
View File

@ -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 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)
sinLambda = Math.sin(lambda);
cosLambda = Math.cos(lambda);
sinSigma = Math.sqrt((cosU2*sinLambda) * (cosU2*sinLambda) +
(dCosU1SinU2 - dSinU1CosU2 * cosLambda) * (dCosU1SinU2 - dSinU1CosU2 * cosLambda));
if (sinSigma==0.0) {
return 0.0;
cosSigma = sU1 * sU2 + cU1 * cU2 * cosLambda;
}
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;
sinAlpha = dCosU1CosU2 * sinLambda / sinSigma;
cosSqAlpha = 1.0 - sinAlpha * sinAlpha;
cos2SigmaM = cosSigma - 2.0 * dSinU1SinU2 / cosSqAlpha;
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);
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);
if (iters == 0)
return 0.0;
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)));
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)));
return this.c * A * (sigma - deltaSigma);
return c * A * (sigma - deltaSigma);
}
@Override

13
lucene/spatial3d/src/test/org/apache/lucene/spatial3d/geom/GeoPointTest.java
View File

@ -67,6 +67,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);
}
}