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.

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

@kbevers
Copy link
Member

@kbevers kbevers commented Aug 24, 2018

Supersedes #931

@cffk
Copy link
Contributor

@cffk 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
Copy link
Member Author

@kbevers 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
Author Member

You guessed right! Nice catch.

@rouault
Copy link
Member

@rouault 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
Copy link
Member Author

@kbevers 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
Copy link
Member

@rouault 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
Copy link
Member Author

@kbevers 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
Copy link
Member Author

@kbevers 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
4 checks passed
4 checks passed
@travis-ci
Travis CI - Pull Request Build Passed
Details
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
@kbevers kbevers deleted the kbevers:inv-mercator-precision branch Feb 6, 2019
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Linked issues

Successfully merging this pull request may close these issues.

None yet

3 participants