Optimization: replace iterative Aitoff inverse with closed-form inverse

[robertkleffner]

• AI (Copilot or something similar) supported my development of this PR. See our policy about AI tool use. Use of AI tools must be indicated.
• Tests added
• Added clear title that can be used to generate release notes
• Fully documented, including updating docs/source/*.rst for new API

Justification

In a side project I've been working on, I uncovered the (already-known) existence of an exact closed-form inverse for the Aitoff projection (in part thanks to an AI agent I was having help build the project). It's fairly straightforward to derive: squaring and adding the forward equations eliminates φ and λ from the right-hand side with some trig identities, leaving the identity (x/2)² + y² = α². Then α in terms of x and y, and thus derivations of φ and λ, follow directly:

α = √((x/2)² + y²)
φ = asin(y·sin α / α)
λ = 2·atan2((x/2)·sin α, α·cos α)

This replaces the nested Newton–Raphson inverse in use since 2015. The closed form matches the iterative result to machine precision, is noticeably faster, needs no special case at the poles, and returns longitudes in [−π, π] by construction.

Winkel Tripel shares this file and unfortunately has no closed-form inverse I'm aware of, despite my attempts. I also attempted seeding the Winkel Tripel iteration with the closed-form inverse, but this ended up hampering performance a bit since the iteration converges so quickly in most cases. So, its Newton–Raphson solver and verification loop are retained, and its results are unchanged.

Correctness

• Every existing Aitoff and Winkel Tripel gie case passes unchanged.
• Added Aitoff round-trip tests over the interior, poles, and anti-meridian, as well as explicit interior inverse tests and a domain-rejection case.

Performance

Measured on an M2 MacBook Air with test/benchmark/bench_proj_trans, inverting via a pipeline:

bench_proj_trans -l 5000000 --noise-x 2000000 --noise-y 1000000 \
  -p "+proj=pipeline +step +inv +proj=aitoff +R=6400000" 5000000 2500000

~1.8 → ~11 M coordinates/s = ~6× speedup end-to-end through proj_trans.

References

• Inverse formulas: J. P. Snyder (1987). Map Projections — A Working Manual.
  USGS Professional Paper 1395, pp. 195–197. https://doi.org/10.3133/pp1395
• Prior implementation of the same inverse (via a pre- and post-scaled
  azimuthal equidistant inverse): Justin Kunimune, Map-Projections,
  https://github.com/jkunimune/Map-Projections

[cffk]

1. Your blog post covering the inverse Aitoff projection is slightly
   overblown. As you note, the Aitoff projection is a pre- and post-scaled
   version of the azimuthal equidistant projection. So the inverse
   projection is the inverse of the azimuthal equidistant projection
   similarly pre- and post-scaled (and this is how Kunimune's library
   implements it). So your formulas are precisely given by Snyder (1987),
   p. 196, with R = 1, phi1 = 0, lam0 = 0, substituting

  x -> x/2
  lam -> 2*lam

2. In your blog post you say "Thus, an easy guard is just to check that
   alpha is ‘close’ to zero; we choose a value of ‘close’ that works with
   the 32-bit floating point numbers used in WebGPU." It's more
   straighforward to make the test "equal to zero", since
   alpha / sin(alpha) is computed accurately for any non-zero small
   alpha.

3. You restrict lam to [-pi, pi] and thus the sphere is projected to
   an ellipse with semiaxes pi and pi/2. However projection is
   perfectly well-behaved at lam = +/- pi, so you might as well accept
   lam in [-2*pi, 2*pi] and the (now double covered) sphere projects to
   an ellipse with semiaxes [2*pi, pi].

4. The resulting inverse projection is

  double rho = hypot(x/2, y),  // prefer hypot(x, y) to sqrt(x^2 + y^2)
    srho = sin(rho), crho = cos(rho);
  if (rho > M_PI || srho < 0)
    // Need to test srho because with long doubles, we have sin(M_PI) < 0.
    return false;
  // Arg to asin is guaranteed in [-1, 1], so don't need aasin
  phi = rho == 0 ? y : asin((y / rho) * srho);
  lam = 2 * atan2(x/2 * srho, rho * crho);

5. Since noone else will be looking at this code for the next decade, it
   would be good if you could also fix the forward projection, while you're
   about it. The use of acos means that it is shockingly inaccurate: the
   point phi = 1e-8 rad, lam = 2e-8 rad (about 0.1 m from the origin)
   projects to [0, 0]. Instead use atan2:

  double clam2=cos(lam/2), slam2=sin(lam/2),
    cphi=cos(phi), sphi=sin(phi),
    crho = cphi * clam2,
    srho = hypot(cphi * slam2, sphi),
    A = srho == 0 ? 1 : atan2(srho, crho) / srho;
  x = 2 * cphi * slam2 * A;
  y = sphi * A;

6. Add tests to check accuracy near the origin and the new bounding
   ellipse.

In summary, this is a worthwhile addition (the use of an iterative solution
was silly in this case). However, it's worth tightening up the numerics.

[robertkleffner]

@cffk Thanks, this is exactly the sort of review I was looking for.

1. Snyder is now cited. I have his other book (Flattening the Earth), looks like I need to read the rest of his work.
2. Exact test adopted, the original guard was an artifact of the WebGPU version I'm using where I'm restricted to f32. Unnecessary for doubles.
3. Great observation, domain extension included.
4. Adopted in full, although I kept the name alpha for uniformity with forward.
5. Adopted in full. The caveat is that the forward is about 30% slower now, but in absolute terms it's small and a worthy price to pay for substantially improved accuracy in key regions of both Aitoff and Winkel Tripel.

I will amend the blog post to properly cite Snyder as well.

[rouault]

1. the Aitoff projection is a pre- and post-scaled
   version of the azimuthal equidistant projection.

shouldn't we then merge aitoff.cpp implementation in aeqd.cpp with 2 new modes in the pj_aeqd_ns::Mode enumeration ? That way we could even get almost for free an ellipsoidal formulation for Aitoff.

[cffk]

Since noone else will be looking at this code for the next decade, it
would be good if you could also fix the forward projection, while you're
about it.  The use of `acos` means that it is shockingly inaccurate: the
point `phi = 1e-8 rad, lam = 2e-8 rad` (about 0.1 m from the origin)
projects to `[0, 0]`.

I'm wrong about the errors in the original code for the forward
projection. Even though alp = acos(cos(phi)*cos(lambda/2)) is
inaccurate near the origin, the error cancels when alp/sin(alp) is
computed. Nevertheless, the replacement code using atan2 is still
needed to maintain accuracy near phi = 0 and lambda = 2*pi.

[cffk]

1. the Aitoff projection is a pre- and post-scaled
   version of the azimuthal equidistant projection.

shouldn't we then merge aitoff.cpp implementation in aeqd.cpp with 2 new modes in the pj_aeqd_ns::Mode enumeration ? That way we could even get almost for free an ellipsoidal formulation for Aitoff.

In principle, this is simple enough to do.

However, I see some downsides:

1. aitoff.cpp already has modes AITOFF and WINKEL_TRIPEL. So someone
   would mesh these with the aeqd modes.

2. aeqd_s_forward needs some cleanup to maintain accuracy for points far
   from the origin.

3. Even though an ellipsoidal generalization of Aitoff is possible,
   perhaps it's better to restrict such "compromise" projections to the
   sphere. After all they are not used for making precise measurements.

So I recommend holding off on combining aitoff and aeqd for now.

[robertkleffner]

2. aeqd_s_forward needs some cleanup to maintain accuracy for points far
   from the origin.

I think I was able to address this as well and have a commit on my local branch for it, touching 2 additional files. However, I'm open to keeping it out of this PR so it can focus on the Aitoff portion, and putting it in a follow-on PR if that order of operations is preferred.

[kbevers]

However, I'm open to keeping it out of this PR so it can focus on the Aitoff portion, and putting it in a follow-on PR if that order of operations is preferred.

It's fine to include it here, just put it in a separate commit.

[robertkleffner]

Done, added the aeqd accuracy improvements alongside verifying tests. It does slow throughput by about 30%.