Add BiasCorr classes and rename previous coreg.BiasCorr in coreg.VerticalShift

[rhugonnet]

Summary

This PR adds non-rigid co-registration "bias-correction" methods through the class BiasCorr and subclasses located in biascorr.py! 😄

In details, this PR:

• Refactors the old vertical shift method coreg.BiasCorr into coreg.VerticalShift.
• Renames the old Coreg parent class into Rigid, as all of its methods (to_matrix, etc) and its subclasses are only relevant for rigid transformations,
• Creates a new "elder" parent class named Coreg that contains only the generic methods to any processing step (e.g., residuals()) and is subclassed by Rigid for rigid transforms, BiasCorr for non-rigid transforms, and CoregPipeline and BlockwiseCoreg to be able to use those with any subclass of Coreg (can now do a pipeline with a Coreg method and a BiasCorr method, for example, and any other class we choose to add below Coreg),
• Rewrite the content of Coreg.fit() into a separate _preprocess function, and uses subsample_raster from GeoUtils to consistently subsample the inputs,
• Streamlines the optimization pipelines in fit.py to be consistent in the way they return optimized parameters and accept user kwargs,
• Adds a fully modular framework based on the choice "fit_or_bin" to run any BiasCorr class: the user can either pass any fit_func (e.g., np.polyval) and fit_optimizer (e.g., scipy.optimize.curve_fit) functions to determine a parametric solution to their bias, or any bin_sizes (e.g., 10), bin_statistic (e.g., np.nanmedian) and bin_apply_method (e.g., "linear" interpolation) to determine an empirical solution to their bias, or to do both (recommended option)!
• Defines robust and case-specific default parameters of "fit_or_bin" for subclasses Deramp, DirectionalBias and TerrainBias to allow user to easily run the classes without customization,
• Fixes several bugs that were in fit.py and coreg.py,
• Adds test_biascorr.py with many tests using synthetic biases,
• Adds a basic documentation page (to be re-worked once coreg.py is updated, see thoughts below),
• Adds a basic gallery example (also to be re-worked once coreg.py is updated).

TO-DO:

• Rename Deramp into Tilt in coreg.py, and fix degree to 1 (pointing to Deramp for degree > 1 ?),
• Add the option for users to do bin_and_fit to simplify big data problems,
• Write documentation page for BiasCorr (to streamline later with updated Coreg?),
• Add gallery examples for BiasCorr (to streamline later with updated Coreg?).

Thoughts on Coreg moving forward

We should probably copy the same logic in Rigid classes:

• Being able to pass any optimizer or binning parameters to resolve a minimization problem parametrically or empirically,
• Reflect these changes and other recent changes in the documentation!

But this should probably live in a separate PR! 😅

Old draft text below

Draft PR: will need PR #151 to be merged to become functional for testing

Mostly trying to organize things logically in biascorr.py right now. It's a bit complex! @erikmannerfelt @adehecq Tell me what you think

Reminder of how biascorr.py differs from coreg.py: integrates all transformations that are not "affine". You can pass any "bias_fun" function to fit the bias along the given dimension. Right now: robust_polynomial or robust_sumsin in 1D.

Structural draft: DirectionalBias and TerrainBias

DirectionalBias can be used to correct in 1D along any angle, notably including along-track and across-track corrections.

TerrainBias can be used to correct in 1D along any variable, notably terrain (elevation, curvature). Because those are inter-connected, I thought it'd be good to have them all in the same module. But we could expand this to any kind of variable really.

But in practice, those are part of a same category of a 1D bias correction with an external variable. It's the exact same operations, except that DirectionalBias is able to derive the binned direction from the np.ndarray directly (no need to pass it to the function) and could possibly find the along/across angle automatically from the shape of the data (would be an amazing functionality 🤯), and TerrainBias can derive the terrain attribute of interest directly from the provided DEMs so we wouldn't need the extra lines to pass it to the function (except if the user wants it).

So... maybe we should make those subclasses of a BiasCorr1D itself subclass or Coreg? So that we don't copy/paste the input/outputs?

Previous BiasCorr -> VerticalShift and previous ZScaleCorr -> TerrainBias

Refactored BiasCorr into VerticalShift in all functions and documentations, and also comments.
Removed ZScaleCorr that'll be replaced by TerrainBias (or the like)

EDIT 09/09/2021:

Regarding the structure, I'm thinking that this is the most logical approach:

• Bias1D:
  including "convenience classes" (that have methods/1D sampling pre-parametrized):
  -> DirectionalBias
  -> TerrainBias

• Bias2D:
  including convenience class:
  -> Deramp

• BiasND (for 3+ dimensions)

Each method can work with a different method of bias prediction:

• Parametric:
  -> Polynomial fit (available for 1D, 2D)
  -> Sum of sinusoid fit (1D only)
• Numeric approximation:
  -> Binning + bilinear interpolation (1D, 2D, ND)
• (TO BE IMPLEMENTED) Non-parametric:
  -> Gaussian Process fit

TODO to finalize PR:

• Homogenize functions in Add robust polynomial, sum of sinusoids fitting #151 to facilitate coefficient retrieval/apply functions
• Structure in: BiasCorr1D, BiasCorr2D, BiasCorrND
• Create convenience subclasses DirectionalBias, TerrainBias, Deramp, Tilt (that last one in coreg.py)
• Add tests on the implementation

[erikmannerfelt]

Very nice!! I think some existing Coreg classes should be in biascorr.py as well. What do you think, @adehecq?

[erikmannerfelt]

Oh, and there's no tests/test_biascorr.py. Should this not be added?

[adehecq]

Seems to be in good shape! I guess you're still waiting to merge #151 before moving further?

[rhugonnet]

Merged with #346

[erikmannerfelt]

Hi @rhugonnet,

This is really the right way forward and I'm grateful that you've spent so much time improving this! While I do have questions on the current setup, I think splitting the Coreg class into rigid/nonrigid methods may be very useful.

1. Could the name Rigid be RigidCoreg or similar? Calling it only Rigid to me sounds like it's just a transform, but it's in reality a co-registration class which estimates a rigid transform.
2. Tilt does not result in a rigid transform, last I checked. Since there's no horizontal transformation, applying it is akin to rotating and scaling at the same time. I can elaborate if that was confusing but we've already discussed this at length in one of the PRs. I think a name like Deramp (like it was) or Detrend is more suitable, as it's not just a rotation.
3. I think RigidCoreg and NonRigidCoreg would make more sense instead of Rigid and BiasCorr. I think it would make it easier to understand how they work for new people getting in to the code. Technically, a rigid transformation corrects a 3D bias, right? So technically, all Coreg classes correct biases, no?
4. If Filter will be subclassed from Coreg in the future, should Coreg change name? Something like just Corr, Correction, or Filter actually, would make more sense, and rigid /non-rigid coreg is subclassed from it. I personally don't think it makes sense to subclass a filter from a co-registration object, but just changing the name of the latter would solve that confusion. Bonus idea, the Corr or whatever class could be in its separate module to make the structure even simpler and to reduce on the already very lengthy coreg.py module.
5. Small thing but coreg.CoregDict is hanging in the air in the inheritance diagram.
6. GradientDescending is missing from the inheritance diagram!

Adding a fit_and_apply() is a great idea I think. It's similar to how OpenCV's feature detection algorithms inherit, detect(), compute() and detectAndCompute() from cv::Feature2D. See ORB for example.

I'll do a code review as well soon (hopefully today). These were just ideas from reading the documentation and your comments!

Again, I'm really excited for this improvement!

[erikmannerfelt]

Is this the largest PR to date? :P

The coreg.py diff has become messed up since things have moved up and down. Therefore it was hard for me to evaluate the last 1/3 of the changes. I guess I'll just have to trust you! If tests pass, however, it can't be that wrong ;)

[erikmannerfelt]

Here and below it's called deramping, while the class itself is now called TerrainBias. It also brings up what I mentioned before about the horizontal scaling, making it non-rigid.

[erikmannerfelt]

Oh yeah. What if instead of Rigid(-Coreg), we call it Affine(-Coreg)? Then, any process that can be explained as an affine matrix would fit into this. It would solve the problem with Tilt not being rigid, but being affine.

[rhugonnet]

Thanks a lot @erikmannerfelt for the thorough read and comments! All accounted for 🙂

We just need to decide on names! I agree with all of your replies 1. to 6. above!

More specifically:

• For the parent Coreg class: I like Corr or Correction for affine/non-affine, but not sure it'd work so well with Filter... maybe Processing or ProcStep but it's a bit generic. We could also keep the name Coreg but make it CoregStep or CoregProc to include Filter? The Filter subclass could also be called CoregFilter, as it will be specific to a co-registration pipeline and a dDEM.
• For Rigid: I like AffineCoreg proposed by @erikmannerfelt, as affine transformations can always be written in a 4x4 matrix, Tilt would fit in it, and it encompasses rigid methods! :)

Looking into the computer vision/DEM literature, it looks like "alignment", "registration" or "co-registration" are synonyms largely used for all kind of transformations (affine and not affine, see https://en.wikipedia.org/wiki/Image_registration). However, nonaffine methods generally imply a certain type of "spatial deformation" is used (spline, quadric, etc...), which is not the case for all of our current BiasCorr classes (terrain) and potentially for other specific bias corrections (radar penetration ?).
I agree that "bias corrections" is a bit generic, it technically applies to all methods... even though "bias" is more often used to refer to what is left uncorrected after affine transformation for DEMs. In our IEEE paper we called those "specific biases" (which didn't mean much except "anything but affine").

So, in short, I see two "base name" options for the parent classes:

1. Data-centered (we refer to the nature of the inputs = two DEMs):
DiffDEMStep, DiffDEMPipeline, DiffDEMBlockwise (or dh or dDEM)

2. Objective-centered (we refer to the pipeline objective that is co-registration):
CoregStep, CoregPipeline, CoregBlockwise

And for the structure, could be:
CoregStep (not public?):

AffineCoreg

ICP
NuthKaab
VerticalShift
Tilt

NonAffineCoreg or BiasCorr

CoherentPointDrift
Deramp
DirectionalBias
TerrainBias
BiasND

CoregFilter

GlobalFilter
KernelFilter

or we could differentiate non-affine registration methods and "bias correction" methods, considering non-affine registration methods are only using the spatial dimension?

CoregStep (not public?):

AffineCoreg

ICP
NuthKaab
VerticalShift
Tilt

NonAffineCoreg

CoherentPointDrift
Deramp
Directional

BiasCorr

TerrainBias
RadarPenetration
BiasND

CoregFilter

GlobalFilter
KernelFilter

What do you think @erikmannerfelt @adehecq?

[adehecq]

Hey,
Wow, that feels so weird to look back at comments that date back as far as 2 years ago !
Thanks a lot for getting back to this @rhugonnet and the huge effort of the past weeks !! 🙌 I think it's heading in a very nice direction. Let me make some generic comments first so we can agree on the structure/names.

1. Regarding the class structure and naming (it's difficult to separate both).

• I think we should have a Pipeline (or ProcPipeline) and a ProcStep living in a separate submodule, that are designed to be generic and be parent classes for coreg, filters and why not interpolation or volume calculation. These should probably be hidden classes, because the user has no reason to modify their behavior.
• I would have a coreg.py module with a CoregAffine and CoregNonAffine subclasses. Note: shouldn't "Coreg" lie before "Affine/NonAffine" because the main characteristic is coreg, then the subtype? I would not have a 3rd BiasCorr subclass as it makes things a bit more complex, but I'm not strongly opposed to it if you're both in favor. I think with a nice diagram, the differences between them should be quite clear.
• I would have a separate filters.py module, with the two suggested FilterGlobal and FilterKernel (again switching the items for consistency) and more if appropriate.
  For me the diagram would look something like this (not sure about how to represent the link between ProcStep and Pipeline and also BlockWise with the other coreg methods). This also made me think. Probably the BlockWise approach does not make sense for filters, because it is the same idea that the FilterKernel, so maybe it should stay in coreg.py?
  

2. Regarding fit/apply, I totally agree that in most cases both steps should be run at the same time. I prefer run rather than fit_and_apply which does not read so well, but if this is consistent with what is done in other libraries like OpenCV, I'm open to using the latest.

3. I had a quick look at the documentation and I really like it! I like that you have the bullet points briefly describing each method at the top. I had a few minor comments on this (e.g. Deramping is also very frequently used for SAR-derived DEMs), but maybe it's easier if I go through the documentation in detail once we have the final terminology/structure.

I'm happy to have a more detailed look at the code later, but again, given that it's a very large PR and will take a significant amount of time, I'd rather do it on a close-to-final version. Plus, @erikmannerfelt already made a first round and I don't think I can compete in term of Python rigorousness 😄

[erikmannerfelt]

I like @adehecq's suggestion on the structure! Depending on the new file length (since the Coreg functionality is now in a proc_step.py file or something) they could live in different python files maybe:

-- xdem
  | - proc_step.py
  | filters/
      | - __init__.py
      | - global_filters.py
      | - kernel_filters.py
  | coreg/
      | - __init__.py
      | - base.py  # Maybe not needed, but `Coreg` could live here if we need extra inheritance from `ProcStep`
      | - affine_coreg.py
      | - non_affine_coreg.py
  | ...

where the new __init__.py would import all classes in the coreg/ and filter/ submodules. Blockwise and Pipeline could live in proc_step.py.

Thoughts?

I'm fine with either run() or fit_and_apply(). If you think run() is better, then lets take that.

@rhugonnet, when this is figured out, I think we should finally merge!!

[rhugonnet]

Thanks both for the great inputs! 😄

I initially agreed with everything and tried to move forward, but then realized there's an issue with this structure 😅: our processing steps and pipelines as they stand are not generic, but specific to two DEMs later differenced, or a DEM and point data, and expect these inputs. So, for instance, we wouldn't be able to apply the same GlobalFilter for a co-registration pipeline, or for filtering a single DEM...
And, if we tried to do so, it could rapidly turn into a mess in terms of consistency of method arguments in the classes, or knowing what object is considered by the filter (or processing step in general).

So we'd need to clearly define to what objects the "steps" and "pipelines" refer to. As of now, we just have things that work on dh, so maybe we should just keep it that way and have dh-centered classes in the pipeline?
I'm not sure it would be very useful to use a filter class on a single DEM anyway, as there isn't many other processing steps to combine them with right now...
Actually, I think the justification we have for building pipelines is specific to co-registration, because of the need to fit a model, and then apply it to one or several outputs independently, and being able to intelligently combine several transformations through associated metadata, and access all of that metadata in the class. I don't think that this is really needed for other processing steps (filtering, terrain, interpolation, volume change)?
In the case of filters for DEMs, maybe it's better to have them as a class method DEM.filter() and a corresponding generic function in filters/? (like the ones we have in terrain.py, we should wrap terrain methods in the DEM class too!)

In short:

• For allowing filters in coregistration pipelines: we could have a specific convenience CoregFilter class that lives in coreg/ and would wrap function-based functionalities in filters/, and could apply them to the ref_dem, or the tba_dem, or both, or the dh?
• For the file structure: building on what @erikmannerfelt laid out, it would keep all the same, and just remove the proc_step.py, add coreg_filters.py in coreg/ and imply that the filters/ module is function-based:

-- xdem
  | filters/
      | - __init__.py
      | - global_filters.py # Only functions
      | - kernel_filters.py # Only functions
  | coreg/
      | - __init__.py
      | - base.py  # To have `CoregStep`, `CoregPipeline` and `BlockwiseCoreg`
      | - affine_coreg.py
      | - non_affine_coreg.py
      | - coreg_filters.py # Classes wrapping filters for coreg
  | ...

• For the naming, as everything would be related to co-registration, we could keep Coreg everywhere. For the order @adehecq mentioned, the adjective seem to usually come first in English (e.g., in pyproj: VerticalCRS or GeographicalCRS are subtypes of CRS), so AffineCoreg and NonAffineCoreg sound better to me:

All inherited from CoregBase:
-- BlockwiseCoreg
-- CoregPipeline
-- CoregStep
  | - AffineCoreg
      | - ICP
      | - NuthKaab
  | - NonAffineCoreg
      | - CPD
      | - Deramp
  | - CoregFilter
      | - MagicFilter
      | - TheBestFilterInTheWorld

What do you think? I hope I didn't miss anything.

If you are both fine with the naming, you can just send a 👍 and I'll start making the changes related to this PR. I agree @erikmannerfelt, we should merge this soon!

[adehecq]

Let's do that then !
Indeed, the Pipeline makes most sense for coregistration and filtering during coregistration. I would argue that it could be interesting for other aspects too, but with some difficult differences in the inputs as you highlighted. Wouldn't it be cool to have a pipeline like:
coreg.NuthKaab() + filters.stdfilter() + volume.hypsometric_filling + volume.volume_change()
but I think we can live with regular line by line processing pipelines 😉
I believe we could get around the issue in the number of inputs by providing a list of rasters rather than 1 or 2 rasters, but it would indeed be a hassle to change all of that and may make things overly complex.

Ok for the order. It also sounded strange to me, but it made more sense to me if we had for example coreg and filters living in the same submodule, because then it would be easier to identify each category (e.g. if we had xdem.AffineCoreg, xdem.GlobalFilter, xdem.NonAffineCoreg etc it would be puzzling). Because they live in two different submodules, it's fine.

[rhugonnet]

Re-structuration done according to comments!

Two points:

• I kept biascorr.py for non-affine/bias-correction methods, because we barely have any common "non-affine" registration methods (Coherent Point Drift), rather they're all corrections. The documentation that will come in the next PR will make the definitions clear.
• I added a pipelines.py module in xdem/coreg/ for the dem_coregistration function (which has to be in another file than base.py, affine.py or biascorr.py, otherwise there are circular import issues), and where we'll also be able to define future "pre-defined" pipelines! 😄

After more than 2 years... merging! 🥳