bug fix in nibabel.affines.rescale_affine

[romainVala]

Hello

I compared the affine matrix after resampling, with different soft (mrtrix, sitk from torchio)
and I notice a small difference in the translation part of the affine when using nibabel function: nib.affines.rescale_affine()

I could find the same result if I modify nibabel code like that

Instead of

    centroid = nib.affines.apply_affine(affine, (shape -1) // 2)
    t_out = centroid - rzs_out @ ((new_shape -1) // 2)

I chage to

    centroid = nib.affines.apply_affine(affine, (shape -1) / 2)
    t_out = centroid - rzs_out @ ((new_shape -1) / 2)

it make sense to me to stay with float value for the centroid coordinate.

Cheers
Romain

[effigies]

I wouldn't call it a bug, but a design decision. The goal was to preserve the RAS coordinate of the central voxel, as opposed to the central point, which could be in the middle of a voxel or on an edge. If I can recall my thinking, it was that interpolation would be unnecessary as long as voxel sizes remained the same, while using the central point would depend on whether the size of each dimension in the input and output images had matching parity.

Perhaps it would make sense to make this a flag?

[romainVala]

Hi
thanks for the explanation, I agree this is not a bug.
I always have difficulty to fully understand, when it goes to the exact position of the volume. Since we are talking about a half (or a quarter) of the voxel size, this is not easy to visually check.

I can understand that there are different convention, but I thing there is only one exact solution for position of the volume so that the resampled volume overlayed on the original volume will show the same position. (right ?)

Matching the central point looks to me like this exact solution ... but I may be wrong

[effigies]

If we resample an image from one space (affine + shape) to another, what we will do is take the $(i, j, k)$ coordinates in the target volume, convert to RAS using the target affine and then find the $(i, j, k)$ coordinates in the source image using the inverse source affine.

$$ \left[\array{i \\ j \\ k \\ 1}\right]_s = \mathbf A_s^{-1} \mathbf A_t \left[\array{i \\ j \\ k \\ 1}\right]_t $$

When the resulting values are not integral, then we need to interpolate (or extrapolate, in cases where the target lies outside the original field of view).

When you plot, you are doing the same thing: Choosing a collection of points to visualize, either by selecting a plane in a reference image or in RAS space, then mapping the points into each of the voxel arrays you want to visualize, and interpolating as needed. So if the resampling is done correctly, viewing one image and its resampled image, we should see identical values, up to interpolation error.

If you are simply overlaying voxels, then you'll see a shift because you aren't resampling into a common RAS space. Most viewers I can think of (mricron, fsleyes, freeview) do this correctly. If you have a plotting function that uses nearest-neighbor interpolation in order to plot two images with translations that are off by a half voxel, it will probably appear as a shift.

[effigies]

Perhaps better to say, if we resample from one space to another, we map into RAS, transform from the target RAS to the source RAS, and then map into the source IJK. If two images are in the same RAS space, then the transform is the identity, so all the affines are doing is telling us the difference in how their voxel arrays are sampling that single space.

[romainVala]

Well actually, I am still not convince that this is only a convention issue.
Let's take a simple example (a case where the resample grid indeed match the original one, with only a zoom factor 2)
imagine a 4x4x4 volume with 1mm voxel size that we resample to a voxel size of 2mm

here is a reproducible code

import nibabel as nib
import numpy as np

def rescale_affine_corrected(affine, shape, zooms, new_shape=None):
    """copy (modify from nibable
    """
    shape = np.asarray(shape)
    new_shape = np.array(new_shape if new_shape is not None else shape)

    s = nib.affines.voxel_sizes(affine)
    rzs_out = affine[:3, :3] * zooms / s

    # Using xyz = A @ ijk, determine translation
    centroid = nib.affines.apply_affine(affine, (shape -1) / 2)
    t_out = centroid - rzs_out @ ((new_shape -1) / 2)
    return nib.affines.from_matvec(rzs_out, t_out)


aa =np.arange(64).reshape((4,4,4))

affine = np.eye(4)
img = nib.Nifti1Image(aa.astype(np.int16), affine)
nib.save(img,'test.nii')

#manual resample, 
aa_resample = aa.reshape((2,2,2,8)).mean(axis=3)

#get the resample affine two ways :

affine_nibabel_resample = nib.affines.rescale_affine(affine, (4,4,4) , 2, (2,2,2) )
affine_corrected = rescale_affine_corrected(affine, (4,4,4) , 2, (2,2,2) )

#now save all volumes
img_orig = nib.Nifti1Image(aa.astype(np.int16), affine)
nib.save(img,'test_1iso.nii')

img_r1 = nib.Nifti1Image(aa_resample.astype(np.int16), affine_nibabel_resample)
nib.save(img_r1,'test_2iso_nibabel.nii')

img_r2 = nib.Nifti1Image(aa_resample.astype(np.int16), affine_corrected)
nib.save(img_r2,'test_2iso_corrected.nii')


print(affine_nibabel_resample)
print(affine_corrected)

[[2. 0. 0. 1.]
 [0. 2. 0. 1.]
 [0. 0. 2. 1.]
 [0. 0. 0. 1.]]

[[2.  0.  0.  0.5]
 [0.  2.  0.  0.5]
 [0.  0.  2.  0.5]
 [0.  0.  0.  1. ]]

now viewing original volume with the resample on in overlay. Using actual nibabel rescale_affine function.

This is not correct, there is a shift, which is not the case for the corrected affine: superposition is exactly matching original volume

Note that using nibable resample, make the same error in the affine (kind of logical) but I also get an unexpected shape ... but this is an other issue (or a concequence of this one .. ? )

from nibabel.processing import resample_to_output
r_img = resample_to_output(img_orig, voxel_sizes=2)
print r_img.shape

Out[262]: (3, 3, 3)

[romainVala]

@effigies did I miss something ?

[effigies]

img_r1 = nib.Nifti1Image(aa_resample.astype(np.int16), affine_nibabel_resample)
nib.save(img_r1,'test_2iso_nibabel.nii')

img_r2 = nib.Nifti1Image(aa_resample.astype(np.int16), affine_corrected)
nib.save(img_r2,'test_2iso_corrected.nii')

You can't "resample" by saving the same array with two different affines. The test should use some kind of resampler, such as in nibabel.processing.

Try the following:

import nibabel as nib
import nibabel.processing
import numpy as np
import scipy.ndimage


def rescale_affine_corrected(affine, shape, zooms, new_shape=None):
    """copy (modify from nibable"""
    shape = np.asarray(shape)
    new_shape = np.array(new_shape if new_shape is not None else shape)

    s = nib.affines.voxel_sizes(affine)
    rzs_out = affine[:3, :3] * zooms / s

    # Using xyz = A @ ijk, determine translation
    centroid = nib.affines.apply_affine(affine, (shape - 1) / 2)
    t_out = centroid - rzs_out @ ((new_shape - 1) / 2)
    return nib.affines.from_matvec(rzs_out, t_out)


rng = np.random.default_rng()

# Generate a smooth random field to have some features to look at
noise = rng.random((20, 20, 20), dtype='f4')
smoothed = scipy.ndimage.gaussian_filter(noise, 4)
scaled = ((smoothed - smoothed.mean()) / smoothed.std()) * 20 + 50

affine = np.eye(4)
affine[:3, 3] = -10  # Center image
orig = nib.Nifti1Image(scaled, affine)
orig.to_filename("orig.nii")

# If input shape is even, target shape should be odd to show a difference
int_aff = nib.affines.rescale_affine(affine, (20, 20, 20), 1, (9, 9, 9))
float_aff = rescale_affine_corrected(affine, (20, 20, 20), 1, (9, 9, 9))

int_resampled = nib.processing.resample_from_to(orig, ((9, 9, 9), int_aff))
float_resampled = nib.processing.resample_from_to(orig, ((9, 9, 9), float_aff))
int_resampled.to_filename("int_resampled.nii")
float_resampled.to_filename("float_resampled.nii")

In the following images, you should see that both image arrays line up with the original, but with the integer one, the voxel edges line up and the values are in fact the same (since I used the same voxel size), while the float one is off by half.

Either way, with a valid resampling tool, the target affine does not affect the validity of the resampling.

Int resampled

Float resampled:

[romainVala]

You can't "resample" by saving the same array with two different affines. The test should use some kind of resampler, such as in nibabel.processing.

I do not really understand why. The whole point of this function is to compute the correct affine, for a given resampling. (independently of the exact resampling done). In my example I do the resampling just after the comment #manual resample,

I came to this function for a specific use case, where I want to compute the resampling myself, so I need to compute the correct affine. (the reason why is related to this discussion #1264).

So the question of interest here, is not about resampling, but how to compute the affine for a given resampling

In the example you chose, the question is then: where should be defined the new grid (size 9^3), compare to the old one of size 20^3
Actually there are many possibility, to position on 9^3 grid on the 20^3 (with same resolution). This example looks more like a cropping problem, ... and there are many possible crop ...! So my understanding, is that there is also an implicit constrain in order to be able to find only one solution : the two grid should share the same center ...
If you agree with this implicit constrain then the corrected solution is the correction I proposed ... ( The actual implemented solution give a shifted grid (I do no see why do you want an assymetry with 6 voxel on the left, and 5 voxel on the right ... ??? .)
Without this constrain, then there are many possible solution ...

On the opposite, for the example I choose, there is an explicit unique solution. In my case I start from an even grid, and ask a zoom factor of 2. The resampled grid is then exactly aligned with the original one. Each voxel on the new grid, is the average of the 2^3 voxel of the original grid. ( we did not agreed on the exact resampling on #1264) but again, we do not care on the exact resampling here, just the grid position

This is why I choose this example, since the true solution is easy to find: the two grid should share the exact same border

[effigies]

I do not really understand why. The whole point of this function is to compute the correct affine, for a given resampling. (independently of the exact resampling done). In my example I do the resampling just after the comment #manual resample,

The point of the function is to provide an affine with the target shape and zooms while preserving the rotations and shears, which can then be used with a generic resampling function. There is not a unique resampling for which to find the affine.

As I suggested in the first post, we could parameterize this function to use // or / based on whether you are trying to preserve the image center or the RAS coordinate of a central voxel (which will be slightly off the image center for dimensions with even numbers of voxels). Or we could provide a second function.

[romainVala]

Hi
my apology for the confusion. ... I do agree with you.
This function aims to define an affine for a given target and zooms, let me just add : while preserving either RAS coordiante of a central voxel, or image center.

Our disagreement remains, to the fact that I still believe there is a unique solution if one want to preserve the alignment between the resample and the original image
( then one needs to preserve the image center )

by the way
what is the central voxel for even number of voxel ?
For instance if there are 4 voxel in the volume, is the central voxel the second or the third voxel ?

If you still disagree, then it is fine, (I may be wrong) and I'll be happy with the parametrized solution

many thanks for your time

[effigies]

I agree there's a unique solution if you're trying to preserve the field of view, modulo the voxel size. It's just that that wasn't the goal of the function when written. In the context of image conformation (#853) the goal was explicitly to set a target FoV, including voxel size and shape. The decision then was to attempt to minimize interpolation when possible, since conforming is usually the first step in a processing chain. And now changing the default behavior could change somebody's results, hence my reluctance to adjust this without either making a new function or a parameter.

And the central voxel is just rounded toward the origin from the image center, so for even dimensions.

[romainVala]

Great, we did agreed at the end

thanks for pointing the thread ,I know better understand the different goals

Then yes, the best is to add an option, and guide the user, with a proper documentation ()

[romainVala]

I forget to ask if this "issue' is relate to the reported behaviour of resample (Cf the end of the code example)

Note that using nibable resample, make the same error in the affine (kind of logical) but I also get an unexpected shape ... but this is an other issue (or a concequence of this one .. ? )

[effigies]

No, I believe that results from this line:

nibabel/nibabel/spaces.py

Line 91 in 3b1c7b3

out_shape = np.ceil((out_mx - out_mn) / out_vox) + 1

I'm guessing it's a fudge factor intended to ensure that there was no cropping at all. Possibly the effect was slightly different in Python 2, when integer division happened if both arguments were integers.