Skip to content

feat: implement to_wgs84_extended#85

Merged
jonhoo merged 8 commits into
helsing-ai:mainfrom
jpteb:feat/closed-wgs84-conv
Apr 7, 2026
Merged

feat: implement to_wgs84_extended#85
jonhoo merged 8 commits into
helsing-ai:mainfrom
jpteb:feat/closed-wgs84-conv

Conversation

@jpteb

@jpteb jpteb commented Mar 8, 2026

Copy link
Copy Markdown
Contributor

Implement a closed-form ECEF to WGS84 algorithm based on the following paper: https://link.springer.com/article/10.1007/s00190-024-01821-w by Quan&Zhang (2024).

This was proposed in #65, as this algorithm is newer and does not have the same height limitations as the current implementation. The paper states accurate calculations between -6.33e6 m and 1e10 m, much better than the current -10000 m to 50000 m. A trade off is a ~57% runtime increase on my machine (47.4ns vs 74.6ns). This PR also introduces a benchmark to test both implementations against each other. There is an optimization in the paper (Appendix 1 Step 2) that is currently not implemented as it relies on a threshold describing the distance of the point to the astroid. This threshold is only defined as "near" by the paper. Testing a naive threshold improved the performance by ~13.5% which brought the increase down to ~38.9% (47.4ns vs 65.8ns). As I do not know at which threshold the optimization can be used, it is not used as of now.

Below is a violin plot created by Criterion:

violin

@jpteb jpteb force-pushed the feat/closed-wgs84-conv branch from 6539462 to 09a1b09 Compare March 8, 2026 00:04
Implement a closed-form ECEF to WGS84 algorithm based on the following
paper: https://link.springer.com/article/10.1007/s00190-024-01821-w by
Quan&Zhang (2024). This was proposed in helsing-ai#65, as this algorithm is newer
and does not have the same height limitations as the current
implementation. The paper states accurate calculations between -6.33e6m
and 1e10m, much better than the current -10_000m to 50_000m. A trade off
is a ~57% runtime increase on my machine (47.4ns vs 74.6ns). This commit
also introduces a benchmark to test both implementations against each
other. There is an optimization in the paper (Appendix 1 Step 2) that is
currently not implemented as it relies on a threshold describing the
distance of the point to the astroid. This threshold is only defined as
"near" by the paper. Testing a naive threshold improved the performance
by ~13.5% which brought the increase down to ~38.9% (47.4ns vs 65.8ns).
As I do not know at which threshold the optimization can be used, it is
not used as of now.
@jpteb jpteb force-pushed the feat/closed-wgs84-conv branch from 09a1b09 to 31fcaba Compare March 8, 2026 14:11

@jonhoo jonhoo left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Thanks for taking this on!

Comment thread src/geodetic.rs Outdated
Comment thread src/geodetic.rs Outdated
Comment thread src/geodetic.rs Outdated
Comment thread src/geodetic.rs Outdated
Comment thread src/geodetic.rs
Comment thread src/geodetic.rs Outdated
Comment on lines +441 to +442
let h =
(w * i + n - SEMI_MAJOR_AXIS * FloatMath::sqrt(FloatMath::powi(i, 2) + n_c)) / s;

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Since we know (?) that e^2 != 0 for the Earth, I believe we can use the "accelerated" formula for h here.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Oh, you are totally right...
Even if this library should get ported to other bodies in the future, it should still hold in reality.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Well, this is a little awkward. With the "optimization" the performance regresses on my machine. 🤔
I'll put both versions in the commit for now and I'd appreciate, if you could test both versions on your machine and see if the behavior is the same for you.

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Hah, that's odd. And presumably they produce the same results?

Comment thread src/geodetic.rs Outdated
Comment thread benches/ecef_to_geodetic.rs Outdated
Comment thread src/geodetic.rs Outdated
Comment thread src/geodetic.rs Outdated

@jpteb jpteb left a comment

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Thank you for the quick review!

The resolved comments will be fixed in the following commit. Otherwise there are still some questions open.

Comment thread benches/ecef_to_geodetic.rs Outdated
Comment thread src/geodetic.rs Outdated
#[doc(alias = "l")]
const L: f64 = (SEMI_MAJOR_AXIS * SEMI_MAJOR_AXIS) * (ECCENTRICITY_SQ * ECCENTRICITY_SQ);
#[doc(alias = "6*l")]
const SIX_L: f64 = 6.0 * L;

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

You're right, that "optimization" is too eager

Comment thread src/geodetic.rs Outdated
Comment thread src/geodetic.rs Outdated
Comment thread src/geodetic.rs Outdated
Comment thread src/geodetic.rs Outdated
#[doc(alias = "l")]
const L: f64 = (SEMI_MAJOR_AXIS * SEMI_MAJOR_AXIS) * (ECCENTRICITY_SQ * ECCENTRICITY_SQ);
#[doc(alias = "6*l")]
const SIX_L: f64 = 6.0 * L;

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Also a style question, I see that you used n. for integral values, would you prefer, if I changed my code to do the same?

Comment thread src/geodetic.rs Outdated
Comment thread src/geodetic.rs Outdated
Comment thread src/geodetic.rs Outdated
Comment thread src/geodetic.rs Outdated
Comment on lines +441 to +442
let h =
(w * i + n - SEMI_MAJOR_AXIS * FloatMath::sqrt(FloatMath::powi(i, 2) + n_c)) / s;

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Well, this is a little awkward. With the "optimization" the performance regresses on my machine. 🤔
I'll put both versions in the commit for now and I'd appreciate, if you could test both versions on your machine and see if the behavior is the same for you.

@jpteb jpteb left a comment

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

I think this should be it.
Comments I wasn't sure about are unresolved, so you can decide on them 👍🏼

Comment thread src/geodetic.rs Outdated
Comment thread src/geodetic.rs
Comment thread src/geodetic.rs
Comment thread src/geodetic.rs Outdated
Comment thread src/geodetic.rs
Comment thread src/geodetic.rs Outdated
Comment thread src/geodetic.rs
@jpteb

jpteb commented Mar 12, 2026

Copy link
Copy Markdown
Contributor Author

I have received an answer from the author of the paper!
And they confirm, that W is missing for the calculation of I, so the current code is correct, and I'll update the comment. They are also working on getting the journal to update the paper to correct this error.

Secondly, on the calculation of t, it is up to us, to decide what we want to use as a close enough threshold. They said the statement is mainly to avoid numerical instabilities if the denominator is equal to zero.

Thinking about this more, in practice we do not care about the values close to the astroid at all. The astroid extends 42700m on the semi-major axis and 42841m on the semi-minor axis from the origin. So it is a pretty small area inside or close to Earth's core.
I propose, that we skip the whole if-else-branch of the calculation of t and use the single cubic root solution always. This would reduce the applicable height from -6.33e6m to -6.31e6m on the semi-minor axis. But I think this is an acceptable trade-off.

See Figure 1 of https://link.springer.com/article/10.1007/s00190-010-0419-x
for reference of the astroid.

@jonhoo jonhoo left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

This looks great now, thank you! And exciting that you heard back from the authors 🎉

Comment thread src/geodetic.rs Outdated
Comment thread src/geodetic.rs Outdated
@jonhoo

jonhoo commented Apr 7, 2026

Copy link
Copy Markdown
Collaborator

After thinking some more about this, I actually think I'd like to make to_wgs84_extended be private (ie, remove pub), and instead have to_wgs84 call to_wgs84_extended when can_convert_to_wgs84 returns false. I'll fiddle with this directly in your branch and then merge so you don't have to deal with it!

@jonhoo jonhoo merged commit 5bb0d6c into helsing-ai:main Apr 7, 2026
16 checks passed
@jonhoo

jonhoo commented Apr 7, 2026

Copy link
Copy Markdown
Collaborator

Releasing in #94 🎉

@jpteb

jpteb commented Apr 12, 2026

Copy link
Copy Markdown
Contributor Author

It was a pleasure! Thank you for the great review, finishing up the work and merging it!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants