Add support for more parametrization in coregistration

[rhugonnet]

This PR adds support for many new parametrizations as well as uncertainty propagation for coregistration, which required to start the restructuration of the uncertainty module (previously spatialstats; see #378).

Details

This PR adds several options and optimizations to affine coregistrations, namely:

• The tolerance argument is replaced by tolerance_translation, tolerance_rotation (and tolerance_objective for CPD), and (internally) the coregistration algorithms now all compute the associated statistics at each iteration (magnitude of change in translation, rotation, objective function) and the iteration stops only when all tolerances are reached,
• A trim_residuals option is added to all methods for outlier filtering, including four sub-options: trim_central_statistic, trim_spread_statistic, trim_spread_coverage, which describe how the trimming is done (for instance median +/- 3NMAD), and trim_iterative to either perform the trimming only once, or by re-computing the trimming statistic iteratively until all outliers are removed,
• The fit_optimizer is removed in favor of fit_minimizer (more generic for optimizer supports, can vary loss function details more easily for instance) and the default consequently changes from scipy.optimize.curve_fit() to scipy.optimize.least_squares() (takes residuals = Y - f(X) as input, instead of X and Y independently). For this adaptation to work, we mirror some defaults of curve_fit(): for example, if undefined, the initial values are set as an array of ones (required argument for least_squares),
• The ICP old support for a "lsq_linearized" fit_minimizer is moved to an argument linearized=True, which corresponds the linearization (Taylor expansion) of the "point-to-plane" optimisation from Low (2004). This linearization now also supports any fit_minimizer (previously a plain NumPy least-square matrix calculation),
• The DhMinimize method is deprecated in favor of a linearized=False argument in LZD (and, to keep the old behaviour, would have to pass only_translation=True and original optimizer/loss functions, because the new LZD(linearized=False) now supports rotations as well!). They are actually the same method (both doing a "Least Z-difference" i.e. LZD), with the Rosenholm and Torlegard (1988) formulation being the linearized version; this greatly improves the consistency of naming/arguments in the package,
• A new variant of CPD named "local surface geometry" (LSG) is added through the argument lsg=True. The LSG variant uses the local normals to better constrain the optimization (like "point-to-plane" does for ICP).
• A new anisotropic argument is added to ICP and CPD to optionally perform automated weighting along the axes X/Y/Z, with strategies "xy_vs_z" and "per_axis". This is very useful to better constrain all rotations/translations for data that has very different extents in X/Y and Z (for instance a DEM that is 50km long, but with only 1km of elevation range),
• A sampling_strategy argument is added for point-point methods (ICP, CPD) which describes how a raster is converted to point cloud before coregistration: It can be set to "same_xy" (both inputs are subsampled at the same X/Y), or "independent" (both inputs are independently subsampled at different random points) or "iterative_same_xy" (both inputs are re-subsampled at the same X/Y every iteration of the algorithm by interpolating the raster input). The new "independent" sampling option required the addition of a new _subsample function logic which works on all inputs,
• An internal reproj_same_grid argument is added for preprocessing the input datasets in Coreg.fit(), it is set as True for bias-correction, but False for affine coregistrations (because point-point registration can work even if the datasets don't overlap in Z coordinate, so reprojection should not be enforced). This change adds the requirement to pass independently a ref_transform and tba_transform to subfunctions (previously a single transform valid for both grids, which were always reprojected), so this is modified for all Coreg._fit_rst_rst subfunctions.
• The iterative output metadata (last_iteration saving the last iteration number, and iteration_stats saving the iteration statistics compared to each tolerance) were derived but not saved in the final .meta until now. We now store these statistics, which can be printed in info() by users.

Finally, being finalized, this PR adds:

• An error-aware LZD with generalized least squares and total least squares implementation, requiring a new optional dependency (GPyTorch).

Resolves #596

TO-DO:

• Add tests for reproj_same_grid (and change boolean to actual value in fit()).
• Merge DhMinimize into LZD?
• Add support for LSG-CPD (uses normals) for a fair comparison of the three algorithms,
• Add support to X/Y vs Z anisotropic weighting for ICP/CPD,
• Add support for general minimizers in CPD,
• Make linearization of ICP solve with same minimizers,
• Add formatting of statistics in _iterate_method prints following Limit the number of decimal places when displaying floats with NumPy #760 (comment),

[rhugonnet]

@belletva @marinebcht The changes of this PR, specifically the addition of the internal reproj_same_grid argument in Coreg.fit() (details above), made me think of #873.

In short: Point-point coregistrations (such as ICP and CPD) work even if datasets are not reprojected to the same grid (= they are converted to points, and don't need to overlap in any way). It's actually essential not to reproject them before coregistration. So, in #873, we shouldn't prevent the user from using coregistration when datasets are not on the same grid.

Maybe, in #873, we just leave the work to Coreg() (for instance NuthKaab or LZD should fail without overlap) and catch/adapt errors to make them more digestible for CLI input/users? What do you think?