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

Speed-up +proj=cart +inv #4087

Merged
merged 4 commits into from
Mar 24, 2024
Merged

Speed-up +proj=cart +inv #4087

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
27561c1
Speed-up +proj=cart +inv
rouault Mar 14, 2024
73ebe8e
+proj=cart +inv: use P->ra
rouault Mar 14, 2024
36da998
geocentric_radius(): reduce number of additions and multiplications
rouault Mar 16, 2024
2d22cf6
geocentric_radius(): save a sqrt()
rouault Mar 16, 2024
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
74 changes: 57 additions & 17 deletions src/conversions/cart.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Copy link
Contributor

Choose a reason for hiding this comment

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

Nice trick using b/a, that is always <1, to avoid overflows and improve performance.

Copy link
Member

Choose a reason for hiding this comment

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

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)

Copy link
Member Author

Choose a reason for hiding this comment

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

Thanks. Do you know if there is a conventional short name/letter for the aspect ratio?

Copy link
Member

Choose a reason for hiding this comment

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

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.

Copy link
Member Author

Choose a reason for hiding this comment

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

after all, isn't b_div_a quite explicit ;-) ?

Copy link
Member

Choose a reason for hiding this comment

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

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!

Copy link
Member

Choose a reason for hiding this comment

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

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

Copy link
Contributor

Choose a reason for hiding this comment

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

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.

Copy link
Member

@busstoptaktik busstoptaktik Mar 17, 2024 •
edited
Loading

Choose a reason for hiding this comment

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

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 sinphi) {
/*********************************************************************
Return the geocentric radius at latitude phi, of an ellipsoid
Expand All @@ -121,8 +121,18 @@ static double geocentric_radius(double a, double b, double cosphi,
This is from WP2, but uses hypot() for potentially better
numerical robustness
***********************************************************************/
return hypot(a * a * cosphi, b * b * sinphi) /
hypot(a * cosphi, b * sinphi);
// Non-optimized version:
// const double b = a * b_div_a;
// return hypot(a * a * cosphi, b * b * sinphi) /
// hypot(a * cosphi, b * sinphi);
const double cosphi_squared = cosphi * cosphi;
const double sinphi_squared = sinphi * sinphi;
const double b_div_a_squared = b_div_a * b_div_a;
const double b_div_a_squared_mul_sinphi_squared =
b_div_a_squared * sinphi_squared;
return a * sqrt((cosphi_squared +
b_div_a_squared * b_div_a_squared_mul_sinphi_squared) /
(cosphi_squared + b_div_a_squared_mul_sinphi_squared));
}

/*********************************************************************/
Expand All @@ -147,8 +157,23 @@ static PJ_LPZ geodetic(PJ_XYZ cart, PJ *P) {
/*********************************************************************/
PJ_LPZ lpz;

// Normalize (x,y,z) to the unit sphere/ellipsoid.
#if (defined(__i386__) && !defined(__SSE__)) || defined(_M_IX86)
// i386 (actually non-SSE) code path to make following test case of
// testvarious happy
// "echo 6378137.00 -0.00 0.00 | bin/cs2cs +proj=geocent +datum=WGS84 +to
// +proj=latlong +datum=WGS84"
const double x_div_a = cart.x / P->a;
const double y_div_a = cart.y / P->a;
const double z_div_a = cart.z / P->a;
#else
const double x_div_a = cart.x * P->ra;
const double y_div_a = cart.y * P->ra;
const double z_div_a = cart.z * P->ra;
#endif

/* Perpendicular distance from point to Z-axis (HM eq. 5-28) */
const double p = hypot(cart.x, cart.y);
const double p_div_a = sqrt(x_div_a * x_div_a + y_div_a * y_div_a);

#if 0
/* HM eq. (5-37) */
Expand All @@ -158,18 +183,33 @@ static PJ_LPZ geodetic(PJ_XYZ cart, PJ *P) {
const double c = cos(theta);
const double s = sin(theta);
#else
const double y_theta = cart.z * P->a;
const double x_theta = p * P->b;
const double norm = hypot(y_theta, x_theta);
const double c = norm == 0 ? 1 : x_theta / norm;
const double s = norm == 0 ? 0 : y_theta / norm;
const double b_div_a = 1 - P->f; // = P->b / P->a
const double p_div_a_b_div_a = p_div_a * b_div_a;
const double norm =
sqrt(z_div_a * z_div_a + p_div_a_b_div_a * p_div_a_b_div_a);
double c, s;
if (norm != 0) {
const double inv_norm = 1.0 / norm;
c = p_div_a_b_div_a * inv_norm;
Copy link
Contributor

Choose a reason for hiding this comment

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

I had to write it down in paper to check that they are equivalent.

s = z_div_a * inv_norm;
} else {
c = 1;
s = 0;
}
#endif

const double y_phi = cart.z + P->e2s * P->b * s * s * s;
const double x_phi = p - P->es * P->a * c * c * c;
const double norm_phi = hypot(y_phi, x_phi);
double cosphi = norm_phi == 0 ? 1 : x_phi / norm_phi;
double sinphi = norm_phi == 0 ? 0 : y_phi / norm_phi;
const double y_phi = z_div_a + P->e2s * b_div_a * s * s * s;
const double x_phi = p_div_a - P->es * c * c * c;
const double norm_phi = sqrt(y_phi * y_phi + x_phi * x_phi);
double cosphi, sinphi;
if (norm_phi != 0) {
const double inv_norm_phi = 1.0 / norm_phi;
cosphi = x_phi * inv_norm_phi;
sinphi = y_phi * inv_norm_phi;
} else {
cosphi = 1;
sinphi = 0;
}
if (x_phi <= 0) {
// this happen on non-sphere ellipsoid when x,y,z is very close to 0
// there is no single solution to the cart->geodetic conversion in
Expand All @@ -181,18 +221,18 @@ static PJ_LPZ geodetic(PJ_XYZ cart, PJ *P) {
} else {
lpz.phi = atan(y_phi / x_phi);
}
lpz.lam = atan2(cart.y, cart.x);
lpz.lam = atan2(y_div_a, x_div_a);

if (cosphi < 1e-6) {
/* poleward of 89.99994 deg, we avoid division by zero */
/* by computing the height as the cartesian z value */
/* minus the geocentric radius of the Earth at the given */
/* latitude */
const double r = geocentric_radius(P->a, P->b, cosphi, sinphi);
const double r = geocentric_radius(P->a, b_div_a, cosphi, sinphi);
lpz.z = fabs(cart.z) - r;
} else {
const double N = normal_radius_of_curvature(P->a, P->es, sinphi);
lpz.z = p / cosphi - N;
lpz.z = P->a * p_div_a / cosphi - N;
}

return lpz;
Expand Down