[BUG] Results for default level 1 (translation) vary between runs

[schlegelp]

Describe the bug
I'm not sure if this is a bug or expected but I noticed that the results for the first level of the default recipe vary between runs.

To Reproduce
This uses the example data:

wget -nv https://github.com/andreasmang/nirep/raw/refs/heads/master/nifti/na01.nii.gz
wget -nv https://github.com/andreasmang/nirep/raw/refs/heads/master/nifti/na15.nii.gz

>>> import warpfield
>>> fixed = warpfield.load_data('na01.nii.gz')[0]
>>> moving = warpfield.load_data('na15.nii.gz')[0]
>>> recipe = warpfield.Recipe()
>>> recipe.levels = recipe.levels[:1]
>>> for i in range(3):
>>>   _, warpmap, _ = warpfield.register_volumes(fixed, moving, this_recipe, verbose=False)
>>>   print(f"Run {i + 1}:\n", warpmap.warp_field.flatten())
Run 1:
 [  1.625 -30.75   28.375]
Run 2:
 [-0.125 -3.    -2.25 ]
Run 3:
 [ 0.75  0.   -0.75]

I had a quick hunt for potential culprits (random initialisations, uninitialized memory, desynchronized streams, etc.) but came up empty. For what it's worth, the deformation field seems to stabilize after the second level - so maybe it doesn't matter?

Environment:
Tested on Google Colab

• Python 3.12
• warpfield 0.2.2
• T4 GPU, CUDA 12.8.r12.8

[bjudkewitz]

Interesting! I can see it too (maybe not so surprising, since I also tested on colab). What's also strange: setting fixed and moving to be identical does not result in a near-0 shift and the estimated shift still jumps. The jitter seems to get less with increasing repeats (which for this affine step does't make any sense to me, but maybe I forgot sth important).
I don't have immediate ideas to explain it. But if this only occurs for single-level recipes, while even the first step works for multi-level recipes, then it might have to be an error that only occurs for low-res warpfields, which would point to a potential issue with the extrapolation in _warp_volume_kernel.

[bjudkewitz]

Having just written that, I doubt linear extrapolation works as expected for a size-1 dimension.

[bjudkewitz]

Indeed for a size-1 axis, the line float v1 = arr[ravel3d(shape, xn, yn, zn)]; of the following section reads of bounds memory. This can also explain the unpredictable behavior you observed.

warpfield/src/warpfield/warp.py

Lines 51 to 61 in 338810a

	if (mode == EXTRAP_LINEAR) {
	int xc = max(0, min(shape[0]-1, xi));
	int xn = (xi < 0) ? xc + 1 : (xi >= shape[0]) ? xc - 1 : xc;
	int yc = max(0, min(shape[1]-1, yi));
	int yn = (yi < 0) ? yc + 1 : (yi >= shape[1]) ? yc - 1 : yc;
	int zc = max(0, min(shape[2]-1, zi));
	int zn = (zi < 0) ? zc + 1 : (zi >= shape[2]) ? zc - 1 : zc;
	float v0 = arr[ravel3d(shape, xc, yc, zc)];
	float v1 = arr[ravel3d(shape, xn, yn, zn)];
	return v0 + (v0 - v1); // linear extrapolation
	}

Even though a size-1 axis seems unlike to happen in normal use (basically only in single-level recipes with 1x1x1 warp map blocks) we could guard against this error by clamping xn/yn/zn within bounds, e.g.:

int xn = (xi < 0) ? min(xc + 1, shape[0]-1) : (xi >= shape[0]) ? max(xc - 1, 0) : xc;