Speed-up +proj=cart +inv #4087
Speed-up +proj=cart +inv #4087
Conversation
|
@wblangdon Assuming this PR pass the tests, can you give it a try on your side to measure the performance? |
rouault
force-pushed
the
faster_geodetic
branch
4 times, most recently
from
4d72380 to
923f6f8
Compare
|
On my system, with a RelWithDebInfo build (-O2), |
|
Passes (64 bit) tests here. Speed up on geodetic 37% (GCC 10.2.1 -O3, i7-4790 CPU @ 3.60GHz) |
rouault
force-pushed
the
faster_geodetic
branch
from
923f6f8 to
73ebe8e
Compare
src/conversions/cart.cpp
Outdated
| #else | ||
| return a * | ||
| sqrt(cosphi * cosphi + | ||
| (b_div_a * b_div_a) * (b_div_a * b_div_a) * (sinphi * sinphi)) / |
There was a problem hiding this comment.
Those parenthesis are not needed. Are for readability or helping the compiler?
There was a problem hiding this comment.
That might help the compiler. Otherwise if we do "a * a * a * a", I suspect he's supposed to evaluate it as "((a * a) * a) * a" given the left-to-right evaluation order mandated by C/C++. Grouping by pair will save one multiplication (since the optimizer should be smart enough to evaluate once b_div_a * b_div_a
There was a problem hiding this comment.
actually I've just added a commit to further optimize things
| @@ -112,7 +112,7 @@ static double normal_radius_of_curvature(double a, double es, double sinphi) { | |||
| } | |||
|
|
|||
| /*********************************************************************/ | |||
| static double geocentric_radius(double a, double b, double cosphi, | |||
| static double geocentric_radius(double a, double b_div_a, double cosphi, | |||
There was a problem hiding this comment.
Nice trick using b/a, that is always <1, to avoid overflows and improve performance.
There was a problem hiding this comment.
Note that b/a is also called the aspect ratio of the ellipsoid, which may be more descriptive than b_div_a. I have noted another pair of identities here. There's no reference to literature, but I remember adding the name and the referred ellipsoid method right after having seen it mentioned in an absolutely reliable source (perhaps Rapp, OSU)
There was a problem hiding this comment.
Thanks. Do you know if there is a conventional short name/letter for the aspect ratio?
There was a problem hiding this comment.
Thanks. Do you know if there is a conventional short name/letter for the aspect ratio?
I do not remember one, and as I was too stupid to write down the reference in the Rust Geodesy Bibliography it will be hard to find again - but this paper uses alpha (and swaps the meaning of a and b!), while Rapp uses alpha for the angular eccentricity, and hence b/a = arccos(alpha). So it does not look like there is any firmly established convention.
Also, I see from this Wikipedia page, that it is also common to define the aspect ratio as the larger parameter divided by the smaller, so even in that aspect (!) there does not seem to be any clear convention.
There was a problem hiding this comment.
after all, isn't b_div_a quite explicit ;-) ?
There was a problem hiding this comment.
after all, isn't b_div_a quite explicit ;-) ?
It can probably harder get any more explicit :-) so given the lack of common nomenclature, it is obviously the better choice!
There was a problem hiding this comment.
Nice trick using
b/a, that is always<1, to avoid overflows and improve performance.
... although I think a few of Jupiter's moons are actually prolate, so this will not hold if we at some time accidentally blindly import IAU ellipsoids
There was a problem hiding this comment.
The algorithm is still stable. If it were thaaat eccentric, then using an ellipsoid model wouldn't be a good idea ;)
I'm curious what is the physics behind a spinning prolate moon.
There was a problem hiding this comment.
The algorithm is still stable.
Ah, yes, I now see that it is, since we may safely assume that b never equals the imaginary unit times a :-)
I'm curious what is the physics behind a spinning prolate moon.
Physically, I believe it's as stable as an oblate, as long as the spin axis and the principal axes of inertia are reasonably aligned. Apparently it just happens to spin around the longest axis, for whichever reason $DEITY came up with during the creation of the solar system (and, in the case of Saturn, whichever perturbations the poor moon suffers from the debris of the ring it shepherds)
| double c, s; | ||
| if (norm != 0) { | ||
| const double inv_norm = 1.0 / norm; | ||
| c = p_div_a_b_div_a * inv_norm; |
There was a problem hiding this comment.
I had to write it down in paper to check that they are equivalent.
rouault
force-pushed
the
faster_geodetic
branch
from
ef3502a to
36da998
Compare
|
:-) I am out of the office next week, but do not hesitate to let me know if I can help when I get back, eg timings |
yes, you're feedback on timings would be helpful since my own measurements show only a modest 13% improvement |
|
sure:-) Please do not delay waiting for me, as I say I will be out-of-action for a while... |
|
geodetic() speed up cart.cpp_2d22cf6 relative PROJ-master 37% (GCC 10.2.1 -O3 i7-4790 3.60GHz) |
Alternative implementation to #4067