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.

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

Improve numerical precision of inverse spherical mercator #1105

Merged
merged 1 commit into from Aug 24, 2018

Conversation

Projects
None yet
3 participants
@kbevers
Member

kbevers commented Aug 24, 2018

Supersedes #931

@cffk

This comment has been minimized.

Contributor

cffk commented Aug 24, 2018

Much better would be to replace

pi/2 - 2*atan(exp(x))

by

-atan(sinh(x))

This is accurate everywhere, is demonstrably odd, and doesn't require branches or C99 functions.

@kbevers kbevers force-pushed the kbevers:inv-mercator-precision branch from c40f61f to 3e967d5 Aug 24, 2018

@kbevers

This comment has been minimized.

Member

kbevers commented Aug 24, 2018

This is accurate everywhere, is demonstrably odd, and doesn't require branches or C99 functions.

I knew you would have a better way to do this. Updated the code and it is indeed much simpler. Thanks for the suggestion!

@@ -57,7 +57,7 @@ static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */
static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */
LP lp = {0.0,0.0};
lp.phi = M_HALFPI - 2. * atan(exp(-xy.y / P->k0));
lp.phi = -atan(sinh(-xy.y / P->k0));

This comment has been minimized.

@rouault

rouault Aug 24, 2018

Member

lots of negation there. I guess atan(sinh(xy.y / P->k0)) would do ?

This comment has been minimized.

@cffk

cffk Aug 24, 2018

Contributor

Yes, the two negations cancel.

This comment has been minimized.

@kbevers

kbevers Aug 24, 2018

Member

You guessed right! Nice catch.

@rouault

This comment has been minimized.

Member

rouault commented Aug 24, 2018

Perhaps there should be a new test point in gie tests to demonstrate the increased accuracy ?

@kbevers kbevers force-pushed the kbevers:inv-mercator-precision branch from 3e967d5 to 8dff160 Aug 24, 2018

@kbevers

This comment has been minimized.

Member

kbevers commented Aug 24, 2018

Perhaps there should be a new test point in gie tests to demonstrate the increased accuracy ?

Yeah. The original PR didn't have on and for some reason I can't come up with one that demonstrates that this approach is better. Perhaps it isn't more precise than the original code?

@rouault

This comment has been minimized.

Member

rouault commented Aug 24, 2018

The difference between both approaches only show up at extremely small absolute values

>>> math.pi / 2 - 2 * math.atan(math.exp(-1e-15))
8.881784197001252e-16
>>> math.atan(math.sinh(1e-15))
1e-15

At 1e-10 and above, the relative error is below 1e-7
At 1e-5 and above, the relative error is below 1e-10
Perhaps the advantage of the new formula would be that it involves less floating point operations. Although it is not clear if sinh() implementations are equivalently optimized as the common exp()

@kbevers

This comment has been minimized.

Member

kbevers commented Aug 24, 2018

Yeah, and my tests has been using extremely small inputs. From the below I see a change, but not exactly an improvement. I would expect to the input to be repeated in the output in this example:

# This PR
$ echo 0 1e-15 | ./bin/proj -I -f %.9e +proj=merc +R=1
0.000000000e+00 5.729577951e-14

# PROJ 5.1
$ echo 0 1e-15 | proj -I -f %.9e +proj=merc +R=1
0.000000000e+00 5.088887490e-14
@kbevers

This comment has been minimized.

Member

kbevers commented Aug 24, 2018

Disregard my conclusion from the previous post. With the changes in this PR the precision is in fact improved. I forgot to account for the RAD2DEG conversion that proj does to the output. The output matches the input when that is taken into account.

@kbevers kbevers force-pushed the kbevers:inv-mercator-precision branch from 8dff160 to 27b2deb Aug 24, 2018

@kbevers kbevers force-pushed the kbevers:inv-mercator-precision branch from 27b2deb to 62b8164 Aug 24, 2018

@kbevers kbevers added this to the 5.2.0 milestone Aug 24, 2018

@kbevers kbevers merged commit ee57ac5 into OSGeo:master Aug 24, 2018

3 checks passed

continuous-integration/appveyor/pr AppVeyor build succeeded
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details
coverage/coveralls Coverage increased (+0.02%) to 77.82%
Details
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment