divide by zero error in DiffeomorphicRegistration of small image volumes

[grlee77]

I came across the following traceback while using SyN DiffeomorphicRegistration registration today. The error occured because the input image size was too small for the default level_iters, but one would be very hard pressed to determine that from the following ZeroDivisionError traceback:

  File "/home/lee8rx/anaconda/envs/py35mkl/lib/python3.5/site-packages/dipy/align/imwarp.py", line 1461, in optimize
    self._optimize()

  File "/home/lee8rx/anaconda/envs/py35mkl/lib/python3.5/site-packages/dipy/align/imwarp.py", line 1387, in _optimize
    derivative = self._iterate()

  File "/home/lee8rx/anaconda/envs/py35mkl/lib/python3.5/site-packages/dipy/align/imwarp.py", line 1258, in _iterate
    der = self._get_energy_derivative()

  File "/home/lee8rx/anaconda/envs/py35mkl/lib/python3.5/site-packages/dipy/align/imwarp.py", line 1350, in _get_energy_derivative
    y = [v / ss for v in y]

  File "/home/lee8rx/anaconda/envs/py35mkl/lib/python3.5/site-packages/dipy/align/imwarp.py", line 1350, in <listcomp>
    y = [v / ss for v in y]

ZeroDivisionError: float division by zero

It turns out that this occurs because the cross-correlation metric can't compute anything meaningful if any of the image dimensions are < 2*CCMetric.radius + 1. The current default CCMetric radius is 4 and by default iterations start at an image downsampling factor of 4, so any inputs that have a dimension < 36 will cause this error to occur unless fewer downsampling levels or a smaller radius are manually specified.

The underlying issue is that compute_cc_forward_step_3d or compute_cc_forward_step_2d will return only zeros if any of the image dimensions are < 2*radius + 1.
This is to be expected from the description in the docstring: "The returned vector field will be zero along a boundary of width radius voxels"

A relatively straightforward solution would be to add a dimension check either when calling the CCMetric or when initializing the DiffeoMorphicRegistration (or AffineRegistration) object to raise a more informative error. If I find time I may make a PR, but it could be a while.

[omarocegueda]

Nice catch! thanks for reporting this @grlee77!, I like your proposed fix, I'll implement it asap! =)

[gimmepizza]

Hello, @omarocegueda are you working on this? I want to start contributing and this is an interesting issue to start with. Also, are you all convinced with the straightforward solution proposed by @grlee77?

[grlee77]

Hi @gimmepizza. It has been some time since this issue was updated, so I think it is safe to say no one else is currently working on it. I would happy to help review a PR .

[skoudoro]

it seems to be fixed by #1540, Can we close this issue @grlee77?

[grlee77]

#1540 should avoid the divide by zero error, but not the reason I encountered cc=0 in the first place. Will the change in #1540 be enough that registration would now finish for a dataset with a dimension smaller than 2*CCMetric.radius + 1? The small size was what caused cc=0 in the first place in my use case (see initial comment in this issue). A solution to that case is to reduce the number of levels in the multi-resolution pyramid, for which I had suggested raising a more helpful error message.

I don't have time to test it prior to ISMRM next week, but can try to take a look afterwards. If anyone wants to test in the meantime, just try running on a dataset where one of the dimensions has been cropped to something small like 16.