Bending energy loss is not scale invariant

[ebrahimebrahim]

Bending energy loss is intended to be applied to a vector field to tell us something about the regularity of the image deformation that the vector field defines when interpreted as a displacement map. I think such a measure should not depend on the resolution of the images that are involved, and so it should also not depend on the spatial resolution used to describe the vector field.

The problem: Currently, BendingEnergyLoss does vary with image resolution, scaling roughly as the inverse square of the image resolution.

Why it's a problem: This makes it tricky to use BendingEnergyLoss in practice. Often it is the regularization objective for training a registration model. The weight to put on a regularization objective like this is typically a very sensitive hyperparameter. Having to retune it every time one uses a different image resolution (e.g. in a fully convolutional network where different resolutions can be used) is prohibitively inconvenient. One can always correct for the way in which it scales, but this makes it overly technical to use properly.

Demo: I made a demo notebook here while trying to work out the issue for myself. There's an estimate of the scaling law, and a plot that shows that the estimate agrees with the empirical result. There's no need to look at this notebook to understand the issue or make comments below-- I'm just linking it here as a supplement.

Example: To show an example from the notebook, the following two vector fields are the "same" vector field, but at different resolutions:

The bending loss should not be so different for these two resolutions.

Clarification: I do believe the current implementation is correctly following the bending energy formula. What I'm suggesting is that we should divide out the image size and work with a "unitless" version of that formula instead. This would be more useful.

Clarification 2: Perhaps a better term is "resolution-invariant." When I say "scale-invariant" I'm talking about image scale, not the scale of the displacement vectors alone. A good way to think about it is in terms of units. The components of a displacement vector field might be described in meters, and the underlying spatial locations at which those vectors live are also described in meters. If one converts both to centimeters, then that shouldn't change the bending energy loss. And it wouldn't, if we understand bending energy loss to have units. But normally when dealing with loss functions we do not think about units. It should be made unitless.

Previous discussion: I see that something similar to this issue was previously discussed, but I don't understand how it was resolved. Tagging people involved in that discssion: @kate-sann5100 @YipengHu @wyli

Potential fix: I can contribute a fix-- I suggest one at the bottom of the demo notebook that I linked. However I'd first like to hear if there's agreement about this... Should bending energy loss be scale invariant?

[ebrahimebrahim]

Another thought: maybe there should be a parameter normalize=True in the constructor of BendingEnergyLoss, which provides as a default behavior the scale-invariance I am describing, and which can be turned off to get the "original" bending energy.

[wyli]

Thanks for the detailed analysis, I think it'd be nice to have a normalize=True option to make it resolution agnostic. The demo notebook could be a monai tutorial on this topic.

For the previous discussions, I can't remember the reasoning, perhaps there's a miscommunication between the input displacement field and the normalised relative coordinates (range [-1, 1] for torch grid_sample).

[kate-sann5100]

@ebrahimebrahim Thank you for the nicely written notebook which clearly explains the need of a resolution-invariant implementation of bending energy.

The concern mentioned by @YipengHu in the previous discussion was that the normalisation of ddf to [-1, 1] in the resample process (to fit the pytorch function grid_sample) may affect the bending energy calculation, which is different from your mentioned problem (if I understood your demo correctly). After off-line discussion, we resolved the concern as the computation of bending energy is based on the unnormalised ddf.

The suggested fix appears sensible to me and I agree with @wyli that have a normalize option would be helpful. It would be great if you would like to contribute the fix. Otherwise, I am happy to do it.

[ebrahimebrahim]

was that the normalisation of ddf to [-1, 1] in the resample process (to fit the pytorch function grid_sample) may affect the bending energy calculation, which is different from your mentioned problem (if I understood your demo correctly). After off-line discussion, we resolved the concern as the computation of bending energy is based on the unnormalised ddf.

I see, that makes sense now.

The suggested fix appears sensible to me and I agree with @wyli that have a normalize option would be helpful. It would be great if you would like to contribute the fix. Otherwise, I am happy to do it.

Great, I'll work on cleaning up my fix and contributing it!

Thanks for the comments @wyli and @kate-sann5100