Add CPD, LZD and make ICP Python-based

NOTICE

@@ -74,11 +74,6 @@ Copyright (c) 2012-2023 Anaconda, Inc.
 Website: https://numba.pydata.org/
 License: BSD 2-Clause.
 
-OpenCV (cv2): Computer vision library.
-Copyright (C) 2015-2023, OpenCV Foundation, all rights reserved.
-Website: https://opencv.org/
-License: Apache License 2.0
-
 scikit-image: Image processing in Python.
 Copyright (c) 2009-2022 the scikit-image team.
 Website: https://scikit-image.org/

dev-environment.yml

@@ -48,12 +48,7 @@ dependencies:
   - numpydoc
 
   - pip:
-    # Optional dependency
     - -e ./
 
-    - noisyopt
-    # "Headless" needed for opencv to install without requiring system dependencies
-    - opencv-contrib-python-headless
-
     # To run CI against latest GeoUtils
 #    - git+https://github.com/rhugonnet/geoutils.git

doc/source/api.md

@@ -259,7 +259,9 @@ To build and pass your coregistration pipeline to {func}`~xdem.DEM.coregister_3d
     coreg.VerticalShift
     coreg.NuthKaab
     coreg.DhMinimize
+    coreg.LZD
     coreg.ICP
+    coreg.CPD
 ```
 
 #### Manipulating affine transforms
@@ -268,15 +270,10 @@ To build and pass your coregistration pipeline to {func}`~xdem.DEM.coregister_3d
 .. autosummary::
     :toctree: gen_modules/
 
-    coreg.AffineCoreg.from_matrix
-    coreg.AffineCoreg.to_matrix
-    coreg.AffineCoreg.from_translations
-    coreg.AffineCoreg.to_translations
-    coreg.AffineCoreg.from_rotations
-    coreg.AffineCoreg.to_rotations
-
     coreg.apply_matrix
     coreg.invert_matrix
+    coreg.matrix_from_translations_rotations
+    coreg.translations_rotations_from_matrix
 ```
 
 ### Bias-correction

doc/source/citation.md

@@ -16,7 +16,7 @@ More details are available on each feature page!
 ### Terrain attributes
 
 ```{list-table}
-   :widths: 1 1
+   :widths: 1 2
    :header-rows: 1
    :stub-columns: 1
 
@@ -41,7 +41,7 @@ More details are available on each feature page!
 ### Coregistration
 
 ```{list-table}
-   :widths: 1 1
+   :widths: 1 2
    :header-rows: 1
    :stub-columns: 1
 
@@ -51,41 +51,41 @@ More details are available on each feature page!
      - [Nuth and Kääb (2011)](https://doi.org/10.5194/tc-5-271-2011)
    * - Dh minimization
      - N/A
+   * - Least Z-difference
+     - [Rosenholm and Torlegård (1988)](https://www.asprs.org/wp-content/uploads/pers/1988journal/oct/1988_oct_1385-1389.pdf)
    * - Iterative closest point
-     - [Besl and McKay (1992)](https://doi.org/10.1117/12.57955)
+     - [Besl and McKay (1992)](https://doi.org/10.1117/12.57955), [Chen and Medioni (1992)](https://doi.org/10.1016/0262-8856(92)90066-C)
+   * - Coherent point drift
+     - [Myronenko and Song (2010)](https://doi.org/10.1109/TPAMI.2010.46)
    * - Vertical shift
      - N/A
 ```
 
 ### Bias-correction
 
 ```{list-table}
-   :widths: 1 1
+   :widths: 1 2
    :header-rows: 1
    :stub-columns: 1
 
    * - Method
      - Reference
    * - Deramp
      - N/A
-   * - Directional bias (sum of sinuoids)
+   * - Directional bias (sinusoids)
      - [Girod et al. (2017)](https://doi.org/10.3390/rs9070704)
-   * - Directional bias (other)
-     - N/A
-   * - Terrain bias (maximum curvature)
+   * - Terrain bias (curvature)
      - [Gardelle et al. (2012)](https://doi.org/10.3189/2012JoG11J175)
    * - Terrain bias (elevation)
      - [Nuth and Kääb (2011)](https://doi.org/10.5194/tc-5-271-2011)
-   * - Terrain bias (other)
-     - N/A
    * - Vertical shift
      - N/A
 ```
 
 ### Gap-filling
 
 ```{list-table}
-   :widths: 1 1
+   :widths: 1 2
    :header-rows: 1
    :stub-columns: 1
 
@@ -101,14 +101,14 @@ More details are available on each feature page!
 ### Uncertainty analysis
 
 ```{list-table}
-   :widths: 2 1
+   :widths: 1 1
    :header-rows: 1
    :stub-columns: 1
 
    * - Method
      - Reference
-   * - R2009 (multiple correlation ranges, circular approximation)
+   * - R2009 (nested ranges, circular approx.)
      - [Rolstad et al. (2009)](http://dx.doi.org/10.3189/002214309789470950)
-   * - H2022 (heteroscedasticity, multiple correlation ranges, spatial propagation approximation)
+   * - H2022 (heterosc., nested ranges, spatial propag.)
      - [Hugonnet et al. (2022)](http://dx.doi.org/10.1109/JSTARS.2022.3188922)
 ```

doc/source/coregistration.md

@@ -111,14 +111,20 @@ Often, an `inlier_mask` has to be passed to {func}`~xdem.coreg.Coreg.fit` to iso
      - Description
      - Reference
    * - {ref}`nuthkaab`
-     - Horizontal and vertical translations
+     - Translations
      - [Nuth and Kääb (2011)](https://doi.org/10.5194/tc-5-271-2011)
    * - {ref}`dh-minimize`
-     - Horizontal and vertical translations
+     - Translations
      - N/A
+   * - {ref}`lzd`
+     - Translations and rotations
+     - [Rosenholm and Torlegård (1988)](https://www.asprs.org/wp-content/uploads/pers/1988journal/oct/1988_oct_1385-1389.pdf)
    * - {ref}`icp`
-     - Translation and rotations
-     - [Besl and McKay (1992)](https://doi.org/10.1117/12.57955)
+     - Translations and rotations
+     - [Besl and McKay (1992)](https://doi.org/10.1117/12.57955), [Chen and Medioni (1992)](https://doi.org/10.1016/0262-8856(92)90066-C)
+   * - {ref}`cpd`
+     - Translations and rotations
+     - [Myronenko and Song (2010)](https://doi.org/10.1109/TPAMI.2010.46)
    * - {ref}`vshift`
      - Vertical translation
      - N/A
@@ -305,20 +311,17 @@ _ = ax[1].set_yticklabels([])
 plt.tight_layout()
 ```
 
-(icp)=
-### Iterative closest point
+(lzd)=
+### Least Z-difference
 
-{class}`xdem.coreg.ICP`
+{class}`xdem.coreg.LZD`
 
 - **Performs:** Rigid transform transformation (3D translation + 3D rotation).
 - **Does not support weights.**
-- **Pros:** Efficient at estimating rotation and shifts simultaneously.
-- **Cons:** Poor sub-pixel accuracy for horizontal shifts, sensitive to outliers, and runs slowly with large samples.
+- **Pros:** Good at solving sub-pixels shifts by harnessing gridded nature of DEM.
+- **Cons:** Sensitive to elevation outliers.
 
-Iterative Closest Point (ICP) coregistration is an iterative point cloud registration method from [Besl and McKay (1992)](https://doi.org/10.1117/12.57955). It aims at iteratively minimizing the distance between closest neighbours by applying sequential rigid transformations. If DEMs are used as inputs, they are converted to point clouds.
-As for Nuth and Kääb (2011), the iteration stops if it reaches the maximum number of iteration limit or if the iterative transformation amplitude falls below a specified tolerance.
-
-ICP is currently based on [OpenCV's implementation](https://docs.opencv.org/4.x/dc/d9b/classcv_1_1ppf__match__3d_1_1ICP.html) (an optional dependency), which includes outlier removal arguments. This may improve results significantly on outlier-prone data, but care should still be taken, as the risk of landing in [local minima](https://en.wikipedia.org/wiki/Maxima_and_minima) increases.
+Least Z-difference (LZD) coregistration is an iterative point-grid registration method from [Rosenholm and Torlegård (1988)](https://www.asprs.org/wp-content/uploads/pers/1988journal/oct/1988_oct_1385-1389.pdf).
 
 ```{code-cell} ipython3
 :tags: [hide-cell]
@@ -327,21 +330,55 @@ ICP is currently based on [OpenCV's implementation](https://docs.opencv.org/4.x/
 :  code_prompt_hide: "Hide the code for adding a shift and rotation"
 
 # Apply a rotation of 0.2 degrees in X, 0.1 in Y and 0 in Z
-e = np.deg2rad([0.2, 0.1, 0])
+rotations = np.array([0.2, 0.1, 0])
 # Add X/Y/Z shifts in meters
 shifts = np.array([10, 20, 5])
 # Affine matrix for 3D transformation
-import pytransform3d
-matrix_rot = pytransform3d.rotations.matrix_from_euler(e, i=0, j=1, k=2, extrinsic=True)
-matrix = np.diag(np.ones(4, dtype=float))
-matrix[:3, :3] = matrix_rot
-matrix[:3, 3] = shifts
+matrix = xdem.coreg.matrix_from_translations_rotations(*shifts, *rotations)
 
-centroid = [ref_dem.bounds.left + 5000, ref_dem.bounds.top - 2000, np.median(ref_dem)]
 # We create misaligned elevation data
+centroid = [ref_dem.bounds.left + 5000, ref_dem.bounds.top - 2000, np.median(ref_dem)]
 tba_dem_shifted_rotated = xdem.coreg.apply_matrix(ref_dem, matrix, centroid=centroid)
 ```
 
+```{code-cell} ipython3
+# Define a coregistration based on LZD
+lzd = xdem.coreg.LZD()
+# Fit to data and apply
+aligned_dem = lzd.fit_and_apply(ref_dem, tba_dem_shifted_rotated)
+```
+
+```{code-cell} ipython3
+:tags: [hide-input]
+:mystnb:
+:  code_prompt_show: "Show plotting code"
+:  code_prompt_hide: "Hide plotting code"
+
+# Plot before and after
+f, ax = plt.subplots(1, 2)
+ax[0].set_title("Before LZD")
+(tba_dem_shifted_rotated - ref_dem).plot(cmap='RdYlBu', vmin=-30, vmax=30, ax=ax[0])
+ax[1].set_title("After LZD")
+(aligned_dem - ref_dem).plot(cmap='RdYlBu', vmin=-30, vmax=30, ax=ax[1], cbar_title="Elevation differences (m)")
+_ = ax[1].set_yticklabels([])
+plt.tight_layout()
+```
+
+(icp)=
+### Iterative closest point
+
+{class}`xdem.coreg.ICP`
+
+- **Performs:** Rigid transform transformation (3D translation + 3D rotation).
+- **Does not support weights.**
+- **Pros:** Efficient at estimating rotation and shifts simultaneously.
+- **Cons:** Poor sub-pixel accuracy for horizontal shifts, sensitive to outliers, and runs slowly with large samples.
+
+Iterative closest point (ICP) coregistration is an iterative point cloud registration method from [Besl and McKay (1992)](https://doi.org/10.1117/12.57955) (point-to-point) or [Chen and Medioni (1992)](https://doi.org/10.1016/0262-8856(92)90066-C) (point-to-plane).
+It aims at iteratively minimizing the distance between closest neighbours by applying sequential rigid transformations. For point-to-point, the 3D distance is used for the minimization while, for point-to-plane, the 3D distance projected on the plane normals is used instead. If DEMs are used as inputs, they are converted to point clouds.
+As for Nuth and Kääb (2011), the iteration stops if it reaches the maximum number of iteration limit or if the iterative transformation amplitude falls below a specified tolerance.
+
+
 ```{code-cell} ipython3
 # Define a coregistration based on ICP
 icp = xdem.coreg.ICP()
@@ -365,6 +402,43 @@ _ = ax[1].set_yticklabels([])
 plt.tight_layout()
 ```
 
+(cpd)=
+### Coherent point drift
+
+{class}`xdem.coreg.CPD`
+
+- **Performs:** Rigid transform transformation (3D translation + 3D rotation).
+- **Does not support weights.**
+- **Pros:** Needs fewer samples to work efficiently.
+- **Cons:** Has trouble solving for horizontal translations if X/Y scale differs from Z scale, as GMM variance is the same in all dimensions.
+
+Coherent point drift (CPD) coregistration is an iterative point cloud registration method from [Myronenko and Song (2010)](https://doi.org/10.1109/TPAMI.2010.46).
+It is a probabilistic approach, iteratively fitting Gaussian mixture model (GMM) centroids for the reference point cloud on the target point cloud by maximizing the likelihood of the probability density estimation problem.
+
+
+```{code-cell} ipython3
+# Define a coregistration based on CPD
+cpd = xdem.coreg.CPD()
+# Fit to data and apply
+aligned_dem = cpd.fit_and_apply(ref_dem, tba_dem_shifted_rotated)
+```
+
+```{code-cell} ipython3
+:tags: [hide-input]
+:mystnb:
+:  code_prompt_show: "Show plotting code"
+:  code_prompt_hide: "Hide plotting code"
+
+# Plot before and after
+f, ax = plt.subplots(1, 2)
+ax[0].set_title("Before CPD")
+(tba_dem_shifted_rotated - ref_dem).plot(cmap='RdYlBu', vmin=-30, vmax=30, ax=ax[0])
+ax[1].set_title("After CPD")
+(aligned_dem - ref_dem).plot(cmap='RdYlBu', vmin=-30, vmax=30, ax=ax[1], cbar_title="Elevation differences (m)")
+_ = ax[1].set_yticklabels([])
+plt.tight_layout()
+```
+
 ## Building coregistration pipelines
 
 ### The {class}`~xdem.coreg.CoregPipeline` object
@@ -476,6 +550,8 @@ For both **inputs** and **outputs**, four consistent categories of metadata are
 
 **4. Affine metadata (common to all affine methods)**:
 
+- An input `only_translation` to define if a coregistration should solve only for translations,
+- An input `standardize` to define if the input data should be standardized to the unit sphere before coregistration,
 - An output `matrix` that stores the estimated affine matrix,
 - An output `centroid` that stores the centroid coordinates with which to apply the affine transformation,
 - Outputs `shift_x`, `shift_y` and `shift_z` that store the easting, northing and vertical offsets, respectively.

examples/advanced/plot_norm_regional_hypso.py

@@ -58,13 +58,13 @@
 
 
 # %%
-# To test the method, we can generate a semi-random mask to assign nans to glacierized areas.
+# To test the method, we can generate a random mask to assign nans to glacierized areas.
 # Let's remove 30% of the data.
-random_nans = (xdem.misc.generate_random_field(dem_1990.shape, corr_size=5, random_state=42) > 0.7) & (
-    glacier_index_map > 0
-)
+index_nans = dem_2009.subsample(subsample=0.3, return_indices=True)
+mask_nans = dem_2009.copy(new_array=np.ones(dem_2009.shape))
+mask_nans[index_nans] = 0
 
-random_nans.plot()
+mask_nans.plot()
 
 # %%
 # The normalized hypsometric signal shows the tendency for elevation change as a function of elevation.
@@ -74,7 +74,7 @@
 # Generating it first, however, allows us to visualize and validate it.
 
 ddem = dem_2009 - dem_1990
-ddem_voided = np.where(random_nans.data, np.nan, ddem.data)
+ddem_voided = np.where(mask_nans.data, np.nan, ddem.data)
 
 signal = xdem.volume.get_regional_hypsometric_signal(
     ddem=ddem_voided,
@@ -105,7 +105,7 @@
 # %%
 # We can plot the difference between the actual and the interpolated values, to validate the method.
 
-difference = (ddem_filled - ddem)[random_nans]
+difference = (ddem_filled - ddem)[mask_nans.data]
 median = np.ma.median(difference)
 nmad = xdem.spatialstats.nmad(difference)
 

setup.cfg

@@ -49,7 +49,6 @@ include =
 
 [options.extras_require]
 opt =
-    opencv-contrib-python-headless
     pytransform3d
     noisyopt
     scikit-learn