Note
Go to the end to download the full example code
Iterative closest point coregistration#
Iterative Closest Point (ICP) is a registration methods accounting for both rotation and translation.
It is used primarily to correct rotations, as it performs worse than Nuth and Kääb (2011) for sub-pixel shifts. Fortunately, xDEM provides the best of two worlds by allowing a combination of the two methods in a pipeline, demonstrated below!
Reference: Besl and McKay (1992).
import matplotlib.pyplot as plt
import numpy as np
import xdem
We load a DEM and crop it to a single mountain on Svalbard, called Battfjellet. Its aspects vary in every direction, and is therefore a good candidate for coregistration exercises.
dem = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"))
subset_extent = [523000, 8660000, 529000, 8665000]
dem = dem.crop(subset_extent)
Let’s plot a hillshade of the mountain for context.
xdem.terrain.hillshade(dem).plot(cmap="gray")
To try the effects of rotation, we can artificially rotate the DEM using a transformation matrix. Here, a rotation of just one degree is attempted. But keep in mind: the window is 6 km wide; 1 degree of rotation at the center equals to a 52 m vertical difference at the edges!
rotation = np.deg2rad(1)
rotation_matrix = np.array(
[
[np.cos(rotation), 0, np.sin(rotation), 0],
[0, 1, 0, 0],
[-np.sin(rotation), 0, np.cos(rotation), 0],
[0, 0, 0, 1],
]
)
# This will apply the matrix along the center of the DEM
rotated_dem = xdem.coreg.apply_matrix(dem, matrix=rotation_matrix)
We can plot the difference between the original and rotated DEM. It is now artificially tilting from east down to the west.
diff_before = dem - rotated_dem
diff_before.plot(cmap="coolwarm_r", vmin=-20, vmax=20)
plt.show()
As previously mentioned, NuthKaab
works well on sub-pixel scale but does not handle rotation.
ICP
works with rotation but lacks the sub-pixel accuracy.
Luckily, these can be combined!
Any xdem.coreg.Coreg
subclass can be added with another, making a xdem.coreg.CoregPipeline
.
With a pipeline, each step is run sequentially, potentially leading to a better result.
Let’s try all three approaches: ICP
, NuthKaab
and ICP + NuthKaab
.
approaches = [
(xdem.coreg.ICP(), "ICP"),
(xdem.coreg.NuthKaab(), "NuthKaab"),
(xdem.coreg.ICP() + xdem.coreg.NuthKaab(), "ICP + NuthKaab"),
]
plt.figure(figsize=(6, 12))
for i, (approach, name) in enumerate(approaches):
approach.fit(
reference_elev=dem,
to_be_aligned_elev=rotated_dem,
)
corrected_dem = approach.apply(elev=rotated_dem)
diff = dem - corrected_dem
ax = plt.subplot(3, 1, i + 1)
plt.title(name)
diff.plot(cmap="coolwarm_r", vmin=-20, vmax=20, ax=ax)
plt.tight_layout()
plt.show()
The results show what we expected:
ICP
alone handled the rotational offset, but left a horizontal offset as it is not sub-pixel accurate (in this case, the resolution is 20x20m).NuthKaab
barely helped at all, since the offset is purely rotational.ICP + NuthKaab
first handled the rotation, then fit the reference with sub-pixel accuracy.
The last result is an almost identical raster that was offset but then corrected back to its original position!
Total running time of the script: (0 minutes 3.366 seconds)