Comparing Nuth–Kääb co-registration and direct NMAD-based optimization

[sgascoin]

I examined how Nuth–Kääb co-registration algorithm compares with a more direct strategy that minimizes NMAD using the Nelder–Mead algorithm. Until now I have mostly used NMAD optimization in my case studies, as I found it more straighforward to implement but I see that everyone uses Nuth-Kääb. I evaluated both methods on a synthetic DEM with a known subpixel translation and vertical offset, and on a Copernicus GLO-30 sample with the same translation plus controlled perturbations: a 1° rotation, planar elevation tilt, combined rotation and tilt, and sinusoidal jitter. See the full code in this notebook.

Case	Nuth–Kääb NMAD (m)	NMAD-optimized NMAD (m)	Nuth–Kääb norm error (px)	NMAD-optimized norm error (px)
Synthetic translation	0.072	0.072	0.002	0.010
GLO-30 translation	1.767	1.741	0.001	0.035
GLO-30 rotation only	13.284	12.900	0.207	0.889
GLO-30 planar tilt	6.830	6.796	0.048	0.149
GLO-30 rotation + tilt	13.585	13.358	0.287	0.708
GLO-30 sinusoidal jitter	2.383	2.357	0.013	0.051

For a clean translation, both methods recover the known offset closely. Under rotation and tilt, the residuals remain large, which is expected because both approaches are limited to translation. The direct NMAD optimization yields slightly lower NMAD values (as expected since it minimizes it), but the Nuth–Kääb formulation generally recovers the horizontal shift more accurately (lower norm error). Computation time is also higher for Nelder–Mead, although the difference is not substantial because both algorithms evaluate only a random subset of up to 1e5 pixels to maintain speed.

In conclusion, both methods perform equally well in the tested cases.

Illustration 1: synthetic case

Illustration 2: real-world DEM + synthetic jitter

Nuth-Kääb

NMAD optimization

Residual histograms

[rhugonnet]

@sgascoin Thanks a lot for sharing! This is great to know for everyone here, and also very reassuring as I found similar results on my side.
(I've been working on a coregistration study for a couple years now, to be submitted soon. By the way, given that you've worked on this quite a bit, and that I tested the direct minimization approach in good part because of your initial contribution here, I'd be happy to invite you onboard it if you like! 😉)

To delve into more details into what you present above: After digging a lot into it these past years, I found that it comes down to the fact that the underlying maths for these two methods are pretty much exactly the same (when using the same objective function, so not technically when using the NMAD, would need an L2 loss instead, but close enough). So it makes perfect sense they would agree so well!
The only differences are: One method is explicit (Nuth and Kääb) by using a linearization of the offset (Taylor first-order; which has little effect because convergence towards 0 shift makes the approximation negligible bit by bit), while the other is implicit (Minimization) as it applies that offset by re-evaluating the sub-pixel grid instead.
This is why they are about to be merged together in #759 (the only difference between the two becoming linearized=True/False and loss_function= parameters).

Going deeper, I found that the maths of LZD and Nuth and Kääb are also exactly the same, one is just expressed in polar coordinates (Nuth and Kääb), while the other is expressed in cartesian while accounting for rotations (LZD). So Nuth and Kääb essentially re-invented the same maths, but only for translations (no rotations), 23 years later. Rosenholm & Torlegard had it all figured out in 1988! 🥲

I did find some slight difference in performance there. It seems the conversion in trigonometrical functions of NK reduces the overall accuracy (but not much), so LZD is very slightly better. I couldn't find another reason for the cause.

So based on the above, the impact of these differences is (if I'm not missing anything):

• LZD/NK are naturally faster (factor of around 10 I think) because they can rely on linear solvers. But we actually didn't have this support internally until the recent addition of @fnands here: Linearize fits in Nuth & Kääb, 1D/2D polynomials #940.
• For large datasets, this also means that we could technically scale the regression on billions of points by using out-of-memory linear regression (such as https://docs.dask.org/en/latest/generated/dask.array.linalg.lstsq.html), which is otherwise impossible for non-linear problems. But that might be overkill, as we usually converge pretty well with a large subsample (the main current mechanic to deal with large data; soon completely out-of-memory with the Xarray/Pandas accessors that are live on the main branch right now).
• And it's also easier to analytically propagate uncertainty through it, which is what I worked on a lot the last couple years! 😄

We've been waiting to highlight this in the documentation until the paper is out as a preprint. Its objective is partly to clarify these things, but mostly to do a wider comparison with point registration, and constrain uncertainty propagation/robustness through many different terrain settings (a bit in the same direction as your added jitter above, but not only numerical, also analytical).