Computational performance of coregistration

[rhugonnet]

Following the presentation of @marinebcht yesterday, I tried to look more into the code today. 🕵️
@belletva @adebardo @adehecq

First, the cause is certainly not the early preprocess step during fit(), because this step is exactly the same for all methods.

That leaves the underlying registration algorithm functions, which are composed of 2 steps for each: A _preprocess_subsample step (yes, there are 2 internal preprocess 😅, I'm actually renaming the second into just _subsample for clarity in #759), and the _iterate_method step, which has different inputs/outputs for each method but relies on the same core function.

To explain the large difference in memory usage between LZD/ICP/CPD and NK/DM/VS, a likely culprit could be the respective preprocess-subsample step in each method, as there are 2 slightly different versions of it:

• The _preprocess_rst_pts_subsample function (example here) is used only in LZD/ICP/CPD, and also in VS (VerticalShift),
• The _preprocess_rst_pts_subsample_interpolator (example here that is used in NK/DM).

However, if I remember correctly, VerticalShift did not have large RAM usage, like NK/DM? If so, then this is probably not the right culprit, as it uses the same faction...
Another factor could be that the _preprocess_subsample step can accept auxiliary variables, which might create more memory usage. But, while auxiliary variables exist in the case of LZD/ICP, they don't in the case of CPD. So probably not that either...

The last step is _iterate_method, which is used by all iterative methods (NK, LZD, ICP, CPD). The only difference here between NK and LZD/ICP/CPD is the fact that NK only applies a translation to the data during iteration (normally, no additional memory usage once the interpolator is created), while LZD/ICP/CPD all apply a translation+rotation. Then they perform a similar optimization.
I manually checked the apply_matrix subfunction utilized here, and memory usage is very small for the subsample size we have (500,000):

import xdem
import numpy as np
test_pc = np.ones((3, 500000), dtype=np.float32)
matrix = xdem.coreg.matrix_from_translations_rotations(2, 3, 4, 5, 6, 7)
new_pc = xdem.coreg.base._apply_matrix_pts_mat(test_pc, matrix=matrix)

All of this points towards the fact that the culprit is likely not the _iterate_method, but _preprocess_subsample, but it's unclear why VerticalShift wouldn't be affected...

[rhugonnet]

Actually, there is a difference between _preprocess_subsample applied to VS or to ICP/LZD/CPD: The return_coords argument is used for ICP/LZD/CPD, and not for VS.
This could explain the high-memory usage.

Maybe something in the underlying function here triggers higher memory usage for return_coords=True or a point cloud input...
But the gridded coordinates estimation just creates 2 more matrices the size of the array. A 5000x5000 float32 array should be 100MB. So, it should be +200MB, not +12GB... 🤔

[rhugonnet]

Conclusion: I really don't know! 😅
It would be great to solve this. Ideally, point cloud methods like ICP should actually use less memory than grid-based method (like NK/LZD) as, once we've taken a point subsample of the grid, which is not particularly demanding, that should be the end of large memory usage...

[adehecq]

Nice investigation! When you mentioned return_coords this made me think of meshgrid which is used when providing the coordinates as 2D arrays, as called in geoutils here. If I'm correct, the default behavior of _coords is grid=True so meshgrid is being called each time. My experience with it is that it is very memory consuming, so worth checking.
Is it possible to use the 1D coordinate arrays instead of 2D here?

[rhugonnet]

Turns out after investigation using the profiling module:
The big memory usage is the apply() part, using either griddata or our custom apply. For both, we could chunk them and linearize the code to apply it bit by bit (less efficient in CPU, more efficient in RAM).

[adehecq]

Adding some notes after investigating the problem of a slow LZD for ~large DEM (8000 x 8000 pixels). I found two possible ways to improve the computation time:

• in base _apply_matrix_rst there is a step to convert from raster to EPC that takes a long time here. Inside this geoutils function, there is a comprehension list that becomes enormous for large DEMs at this place. Might be possible to replace with these few lines instead (to check that results are equal):

x_coords, y_coords = source_raster.coords(force_offset=force_pixel_offset)
x_coords_2 = x_coords[indices[0], indices[1]]
y_coords_2 = y_coords[indices[0], indices[1]]

• a few lines later in the same function, the conversion from array to GeoDataFrame and then to EPC is really slow. We might be able to speed up _apply_matrix_rst by using arrays instead of EPC here.