Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

c: fix "a", revert to 2-iter exact #55

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

c: fix "a", revert to 2-iter exact #55

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
7025b78
c: fix "a", revert to 2-iter exact
Artoria2e5 Jan 9, 2018
ff21c70
c: rename absX to sqrtX
Artoria2e5 Jan 12, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
70 changes: 29 additions & 41 deletions c/transform.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,43 +24,44 @@ 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 sqrtX = 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;
*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) {
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) {
Expand Down Expand Up @@ -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 = 3;
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)<threshold) && (fabs(dLng)<threshold) ) {
return;
}
if (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;
}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The I'm glad to see this method changed as it's so much faster.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Side note... It would be really nice to see consistency between versions and variable names. X_Pi,xPi,yPi,MPi.....
and consistent definitions across implementations mPi = 3000Pi/180; xPi = xmPi for example..
exact transform implemented via cajun method everywhere... or left as half distance method, and mars2wgs as cajin. Just so long as it's all consistent.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The C impl is actually the first instance where a caijun-esque method is introduced to eviltransform. Like I said in the PR, it's ditched in an "optimize" commit, the only productive part of which is fixing a "fabs" declaraction.

Copy link
Contributor Author

@Artoria2e5 Artoria2e5 Jan 12, 2018 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And yeah the whole varname thing is crazy. We probably don't even need all those, as many of them like x*pi and sqrt(fabs(x)) are really subexpression optimizations that can be done automatically by the compiler as long as you allow it to do associative math on floats.

Yeah, right, good call -- let me think really hard about putting some pragmas to tell the compilers…

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I like defining subexpressions or constants that otherwise would be recalculated multiple times later , to limit the clock cycles needed for repetitive calculations. I just like to see consistent naming.

}

Expand Down