Non-rigid registration in simpleITK: please help

Beginner Questions
, ,
1

Hi, everyone,
I am trying to register (non-rigidly) two volumes of microscopy images, a head of C.elegans worm imaged with lightsheet microscope from two different directions. I managed to register them rigidly, but anything I try with non-rigid (FFD, Demons) either runs forever or produces highly distorted results. Would anyone be able to look at my code on github and give me some tips? I am quite desperate to make this work. If someone helps me out, I’ll be happy to have him/her as a co-author :wink:

2

Hello @nvladimus,

Looked at the code and I’ve got some questions:

  1. The text says “original images are highly anisotropic in Z vs XY. Resampled in simpleITK for isotropic resolution.” For registration you don’t need to do this, actually you shouldn’t do this. It is only relevant for visualization. As you increase the number of pixels (image size), you are also increasing the number of samples that are used in registration as these are specified as a percentage of the image size (slower registration).
  2. The global transformation you are using is rigid (that code isn’t there). Have you tried an affine transformation? If yes, why did you prefer rigid.
  3. As I’m not familiar with the modality, have you tried SetMetricAsCorrelation instead of SetMetricAsMeanSquares? Mean squares metric assumes the intensities in the images don’t change. Correlation accommodates for linear changes, possibly more robust in your case.

Please try these changes and then we can proceed to the deformable step. The FFD is recommended, but using the Demons metric in the registration framework isn’t that great. The standalone Demons algorithms are faster.

Generally speaking, to tune a registration using a Jupyter notebook:

  1. Attach observers so that you plot the similarity metric (you will remove these once you’re done with the tuning as the plotting often takes longer than the registration).
  2. Use a small number of iterations (<=100 or <=50 depending on how long it takes) so that you can modify settings and quickly observe the convergence behavior.
  3. At the end of registration, print the termination condition to ensure that it is reaching the max number of iterations - if not you can decrease the number of iterations.

If you have reasonable settings you should see the metric going down, otherwise you’ll see it bumping along and likely terminating very quickly. Once the registration looks like it is going in the right direction and you are terminating due to max-iterations you can increase that and let it run.

2 Likes
3

Hi, @zivy!
Thanks for looking into this. After your feedback, I managed to get some FFD and Demons registration working. Very excited about this progress!
To answer your questions,

  1. Anisotropy in this dataset is weird, because original stacks are rotated 45 deg around Y axis. To preserve information, stacks are therefore up-sampled to isotropic resolution before rotation.
  2. I tried affine (see the updated notebook) - the three optimizers I tried crash, I don’t know why. In general, the worm may wiggle between the acquisition time points, and affine would not be enough, some non-linear transform is needed.
  3. Good point, I used SetMetricAsCorrelation this time, and things started to work!

Do you have any tips why affine transformation may crash? Any further tips how to improve the registration?
Many thanks.

4

Hello @nvladimus,

Happy to hear you are making progress.

The suggestion to use affine was not intended to replace FFD, just make it easier to deal with the deformation. The affine takes care of the coarser part and FFD the finer.

MeanSquares is wonderful for very constrained image acquisition situations. In most real world cases even when the intensities are supposed to remain the same they often don’t and that breaks it in terms of convergence.

The registration should never crash! Throwing an exception due to a problem is expected, but crashing is a problem in our toolkit. A quick look at the affine registration code, I see you didn’t call the SetMetricSamplingStrategy and SetMetricSamplingPercentage. Try calling these and let us know if that made a difference (if these prevent the crash then there is an issue with our default settings).

As I’m not familiar with your data, all I see is some gray blurry stuff that looks reasonably aligned (less magenta and green). How good is good enough? If you can point out what isn’t aligned then possibly we can think of how to improve the results.

1 Like
5

Hi, @zivy,
I tried what you suggested, but introducing SetMetricSamplingStrategy and SetMetricSamplingPercentage did not change anything, affine registration still fails with status Optimizer's stopping condition, LBFGSBOptimizerv4: Failure. Not sure what to try next - any suggestions?

The FFD did a reasonable job, actually a bit too much: the volumes don’t overlap fully, but FFD squeezes view1 to overlap more than it should. Is there any way to prevent this? Would using a mask help against such over-fitting?
Thank you!

6

Hi @nvladimus,

Can you share the data you are using for registration? Or a similar dataset that recreates the problems that you experience with the affine registration?

I can try to fiddle a bit with the affine registration. That has to work. I’d also like to recreate the crashing behavior so that we can figure out how to fix that. A failure of the optimizer isn’t great behavior but it is acceptable, crashing is just not acceptable.

Finally, with respect to FFD overfitting, the only way we have to mitigate that is with hard constraints, the SetOptimizerAsLBFGSB has lowerBound and upperBound scalar parameters that enforce hard limits on the parameter values. Set these based on the parameter values you observe as the output of the current overfitting registration.

7

I would not recommend using an LBFGS optimizer for the low number of parameters in a Affine transform. Consider using a variant of a gradient decent such as the GradientDescentLineSearch.

1 Like
8

@zivy, the repo I posted has the data, file total size 72M. Please, have a look!

@zivy, @blowekamp so, the SetOptimizerAsGradientDescentLineSearch() and SetOptimizerAsLBFGSB() have upper bound limit argument that applies to ALL parameters of the transform matrix? Is this how it works?
Thanks!

[update: repo link is corrected]

9

OK, so the only optimizer that has bounds constraints is the LBFGSB and it is best suited for the FFD transform optimization as it has a large number of parameters (the L stands for limited-memory so efficient memory-wise approximation of the BFGSB algorithm). As @blowekamp commented, it is not appropriate for the low dimensional transformation types (affine, rigid,…). The same bounds are for all parameters so this is also most appropriate for the FFD transformation, as the parameters are the displacements of the sparse grid of control points [dx0,dy0,dz0,dx1,dy1,dz1…].

The link you used for the repo is to a U of Glasgow zoom meeting service - I guess everything is supposed to be done via zoom nowadays, including image registration :wink: . I did go to this repo, so will take a look at the affine registration shortly.

1 Like
10

Hello @nvladimus,

This pull request solves your problem using a three step registration workflow. I kept it as simple as possible (no multi-resolution pyramids). Some things worth mentioning with respect to the provided solution:

  1. For the global transformations, the call to SetOptimizerScalesFromPhysicalShift is critical as the parameters are not commensurate (combinations of angles in radians, translation in mm/um, and scaling which is unitless) a change of one unit in a parameter has a different effect. This function tells the optimizer to account for this difference.
  2. The FFD is terminated early (low value of numberOfIterations) - change as you see fit.
  3. I did not place constraints on the FFD, you can add these in as discussed above.
  4. The plotting of the similarity metric slows the registration considerably, once your happy with the settings comment the AddCommand lines.
  5. You are working with SimpleITK 1.2.4 or lower, we are very close to releasing SimpleITK 2.0. I suggest you upgrade to the current release candidate. This will require a minor change in your util.py file resampleSliceFilter.Execute becomes sitk.Resample.

Hopefully with this PR and post you’ll be able to move forward with your image analysis.

1 Like
11

Hi, @zivy
Sorry for the wrong repo link - it’s my failure at multitasking :roll_eyes:
Many thanks for looking at the data, your PR and feedback!