Excessive padding when resampling zooms

[mgxd]

When working with RegridToZooms (aka resample_by_spacing), I was going from even -> odd slice lengths when doubling the zooms of an image. After a little digging, it seems to me we should be rounding (maybe to the 3rd decimal spot) before calling np.ceil to avoid any extra slices..

Relevant code section:

niworkflows/niworkflows/utils/images.py

Lines 227 to 232 in 47d6034

	card = nb.affines.from_matvec(np.diag(pre_zooms))
	extent = card[:3, :3].dot(np.array(in_file.shape[:3]))
	card[:3, 3] = -0.5 * extent
	# Cover the FoV with the new grid
	new_size = np.ceil(extent / zooms).astype(int)

Test case:

ipdb> card
array([[0.8, 0. , 0. , 0. ],
       [0. , 0.8, 0. , 0. ],
       [0. , 0. , 0.8, 0. ],
       [0. , 0. , 0. , 1. ]], dtype=float32)
ipdb> pre_zooms
array([0.8, 0.8, 0.8], dtype=float32)
ipdb> in_file.shape[:3]
(208, 300, 320)
ipdb> extent
array([166.40000248, 240.00000358, 256.00000381])
ipdb> zooms
array([1.6, 1.6, 1.6])
ipdb> extent / zooms
array([104.00000155, 150.00000224, 160.00000238])
ipdb> new_shape
array([105, 151, 161])
# I would expect a shape of [104, 150, 160] here

Does this make sense, or is it fine to just include the extra slice? cc @oesteban

[effigies]

That seems reasonable to me.

Just to derail a little bit, this function worries me. It's a complicated, one-off procedure that implements a lot of things that should be composable operations. It seems like it would be better implemented as:

def resample_to_zooms(img, target_zooms, **kwargs):
    shape, affine = rezoom(img.affine, target_zooms)
    new_img = resample_img(img, shape, affine, **kwargs)
    # Possibly do some header fiddling
    return new_img

def resample_img(img, shape, affine, order=3):
    ...

def rezoom(affine, target_zooms):
    ...

Then we could even parameterize rezoom to determine what happens in these edge cases, rather than opaquely adjusting a call to np.ceil(). Or make it clear in the docstring that resample_to_zoom will always round up and that the caller should check the expected shape and possibly use img.slicer[:-1, :-1, :-1] to get the desired shape.

[mgxd]

agreed - it will also simplify testing

[oesteban]

+1