From 7025b78bd3470fb29f0a4c429041ad2e324154a8 Mon Sep 17 00:00:00 2001 From: Mingye Wang Date: Tue, 9 Jan 2018 13:52:51 +0800 Subject: [PATCH 1/2] c: fix "a", revert to 2-iter exact Commits cd7a6a4 killed the 2-iter exact algo. This commit brings it back. --- c/transform.c | 64 +++++++++++++++++++++------------------------------ 1 file changed, 26 insertions(+), 38 deletions(-) diff --git a/c/transform.c b/c/transform.c index 5814ca8..5c5617f 100644 --- a/c/transform.c +++ b/c/transform.c @@ -24,26 +24,26 @@ INLINE static int outOfChina(double lat, double lng) { return 0; } -#define EARTH_R 6378137.0 +#define EARTH_R 6371000 void transform(double x, double y, double *lat, double *lng) { double xy = x * y; double absX = sqrt(fabs(x)); double xPi = x * M_PI; double yPi = y * M_PI; - double d = 20.0*sin(6.0*xPi) + 20.0*sin(2.0*xPi); + double d = 2.0*sin(6.0*xPi) + 2.0*sin(2.0*xPi); *lat = d; *lng = d; - *lat += 20.0*sin(yPi) + 40.0*sin(yPi/3.0); - *lng += 20.0*sin(xPi) + 40.0*sin(xPi/3.0); + *lat += 2.0*sin(yPi) + 4.0*sin(yPi/3.0); + *lng += 2.0*sin(xPi) + 4.0*sin(xPi/3.0); - *lat += 160.0*sin(yPi/12.0) + 320*sin(yPi/30.0); - *lng += 150.0*sin(xPi/12.0) + 300.0*sin(xPi/30.0); + *lat += 16.0*sin(yPi/12.0) + 32.0*sin(yPi/30.0); + *lng += 15.0*sin(xPi/12.0) + 30.0*sin(xPi/30.0); - *lat *= 2.0 / 3.0; - *lng *= 2.0 / 3.0; + *lat *= 20.0 / 3.0; + *lng *= 20.0 / 3.0; *lat += -100.0 + 2.0*x + 3.0*y + 0.2*y*y + 0.1*xy + 0.2*absX; *lng += 300.0 + x + 2.0*y + 0.1*x*x + 0.1*xy + 0.1*absX; @@ -53,14 +53,15 @@ static void delta(double lat, double lng, double *dLat, double *dLng) { if ((dLat == NULL) || (dLng == NULL)) { return; } + const double a = 6378245.0; const double ee = 0.00669342162296594323; transform(lng-105.0, lat-35.0, dLat, dLng); double radLat = lat / 180.0 * M_PI; double magic = sin(radLat); magic = 1 - ee*magic*magic; double sqrtMagic = sqrt(magic); - *dLat = (*dLat * 180.0) / ((EARTH_R * (1 - ee)) / (magic * sqrtMagic) * M_PI); - *dLng = (*dLng * 180.0) / (EARTH_R / sqrtMagic * cos(radLat) * M_PI); + *dLat = (*dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * M_PI); + *dLng = (*dLng * 180.0) / (a / sqrtMagic * cos(radLat) * M_PI); } void wgs2gcj(double wgsLat, double wgsLng, double *gcjLat, double *gcjLng) { @@ -94,35 +95,22 @@ void gcj2wgs(double gcjLat, double gcjLng, double *wgsLat, double *wgsLng) { } void gcj2wgs_exact(double gcjLat, double gcjLng, double *wgsLat, double *wgsLng) { - double initDelta = 0.01; - double threshold = 0.000001; - double dLat, dLng, mLat, mLng, pLat, pLng; - dLat = dLng = initDelta; - mLat = gcjLat - dLat; - mLng = gcjLng - dLng; - pLat = gcjLat + dLat; - pLng = gcjLng + dLng; + double dLat, dLng; + // n_iter=2: centimeter precision, n_iter=5: double precision + const int n_iter = 2; int i; - for (i=0; i<30; i++) { - *wgsLat = (mLat+pLat) / 2; - *wgsLng = (mLng+pLng) / 2; - double tmpLat, tmpLng; - wgs2gcj(*wgsLat, *wgsLng, &tmpLat, &tmpLng); - dLat = tmpLat - gcjLat; - dLng = tmpLng - gcjLng; - if ( (fabs(dLat) 0) { - pLat = *wgsLat; - } else { - mLat = *wgsLat; - } - if (dLng > 0) { - pLng = *wgsLng; - } else { - mLng = *wgsLng; - } + if ((wgsLat == NULL) || (wgsLng == NULL)) { + return; + } + *wgsLat = gcjLat; + *wgsLng = gcjLng; + if (outOfChina(gcjLat, gcjLng)) { + return; + } + for (i = 0; i < n_iter; i++) { + delta(*wgsLat, *wgsLng, &dLat, &dLng); + *wgsLat = gcjLat - dLat; + *wgsLng = gcjLng - dLng; } } From ff21c70adbc231a31aec58a48f5aab00dc6bd2d3 Mon Sep 17 00:00:00 2001 From: Mingye Wang Date: Fri, 12 Jan 2018 23:53:07 +0800 Subject: [PATCH 2/2] c: rename absX to sqrtX Quite frankly I would like to see these all dealt with in some -fassociative-math pragma --- c/transform.c | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/c/transform.c b/c/transform.c index 5c5617f..9de2d89 100644 --- a/c/transform.c +++ b/c/transform.c @@ -28,7 +28,7 @@ INLINE static int outOfChina(double lat, double lng) { void transform(double x, double y, double *lat, double *lng) { double xy = x * y; - double absX = sqrt(fabs(x)); + double sqrtX = sqrt(fabs(x)); double xPi = x * M_PI; double yPi = y * M_PI; double d = 2.0*sin(6.0*xPi) + 2.0*sin(2.0*xPi); @@ -45,8 +45,8 @@ void transform(double x, double y, double *lat, double *lng) { *lat *= 20.0 / 3.0; *lng *= 20.0 / 3.0; - *lat += -100.0 + 2.0*x + 3.0*y + 0.2*y*y + 0.1*xy + 0.2*absX; - *lng += 300.0 + x + 2.0*y + 0.1*x*x + 0.1*xy + 0.1*absX; + *lat += -100.0 + 2.0*x + 3.0*y + 0.2*y*y + 0.1*xy + 0.2*sqrtX; + *lng += 300.0 + x + 2.0*y + 0.1*x*x + 0.1*xy + 0.1*sqrtX; } static void delta(double lat, double lng, double *dLat, double *dLng) { @@ -97,7 +97,7 @@ void gcj2wgs(double gcjLat, double gcjLng, double *wgsLat, double *wgsLng) { void gcj2wgs_exact(double gcjLat, double gcjLng, double *wgsLat, double *wgsLng) { double dLat, dLng; // n_iter=2: centimeter precision, n_iter=5: double precision - const int n_iter = 2; + const int n_iter = 3; int i; if ((wgsLat == NULL) || (wgsLng == NULL)) { return;