Function to rescale an affine matrix

[effigies]

Given that image.header.set_zooms() does not update an affine, a useful tool would be a function that takes a given affine matrix and returns an updated matrix that encodes the desired zooms.

Adapting from @ellisdg's suggestions:

def rescale_affine(affine, target_zooms):
    ret = np.array(affine, copy=True)

    RZS = affine[:3, :3]
    zooms = np.sqrt(np.sum(RZS ** 2, axis=0))
    scale = np.divide(target_zooms, zooms)
    ret[:3, :3] = RZS * np.diag(scale)
    return ret

The correct way to fix the zooms on a file would then be:

import nibabel as nb
orig = nb.load(fname)
img = orig.__class__(np.array(orig.dataobj),
                     rescale_affine(orig.affine, target_zooms),
                     orig.header)
img.to_filename(fname)

This approach would leave the RAS coordinates of (i, j, k) = (0, 0, 0) the same. It may be desirable instead to adjust the translation so that the IJK coordinates of the RAS origin are constant. Possibly controlled by a parameter.

I'm not absolutely positive that this handles shear components correctly. I'd need to check.

Sub-issue identified in #619.

[leej3]

@effigies I think this is a nice way of dealing with the rotation and shear components:

def rescale_affine(input_affine, voxel_dims=[1, 1, 1], target_center_coords= None):
    """
    This function uses a generic approach to rescaling an affine to arbitrary
    voxel dimensions. It allows for affines with off-diagonal elements by
    decomposing the affine matrix into u,s,v (or rather the numpy equivalents)
    and applying the scaling to the scaling matrix (s).

    Parameters
    ----------
    input_affine : np.array of shape 4,4
        Result of nibabel.nifti1.Nifti1Image.affine
    voxel_dims : list
        Length in mm for x,y, and z dimensions of each voxel.
    target_center_coords: list of float
        3 numbers to specify the translation part of the affine if not using the same as the input_affine.

    Returns
    -------
    target_affine : 4x4matrix
        The resampled image.
    """
    # Initialize target_affine
    target_affine = input_affine.copy()
    # Decompose the image affine to allow scaling
    u,s,v = np.linalg.svd(target_affine[:3,:3],full_matrices=False)
    
    # Rescale the image to the appropriate voxel dimensions
    s = voxel_dims
    
    # Reconstruct the affine
    target_affine[:3,:3] = u @ np.diag(s) @ v

    # Set the translation component of the affine computed from the input
    # image affine if coordinates are specified by the user.
    if target_center_coords is not None:
        target_affine[:3,3] = target_center_coords
    return target_affine

[leej3]

@Shotgunosine helped me figure out the scaling of the translation if that is desired too. Something along the lines of:

    target_affine = img.affine.copy()
    
    # Calculate the translation part of the affine
    spatial_dimensions = (img.header['dim'] * img.header['pixdim'])[1:4]
    
    # Calculate the translation affine as a proportion of the real world
    # spatial dimensions
    image_center_as_prop = img.affine[0:3,3] / spatial_dimensions
    
    # Calculate the equivalent center coordinates in the target image
    dimensions_of_target_image = (np.array(voxel_dims) * np.array(target_shape))
    target_center_coords =  dimensions_of_target_image * image_center_as_prop 

    target_affine = rescale_affine(target_affine,voxel_dims,target_center_coords)

And then to tie it altogether to resample an image given a target_shape and voxel dimensions:

    resampled_img = image.resample_img(img, target_affine=target_affine,target_shape=target_shape)
    resampled_img.header.set_zooms((np.absolute(voxel_dims)))

[effigies]

affines.rescale_affine added in #853.

[sitandon]

@effigies : I am trying to figure what will be the new translation values in case of scaling. Lets say we have new scale and old scale. I was going thru @ellisdg's code and he has used something like this (https://github.com/ellisdg/3DUnetCNN/blob/master/unet3d/utils/affine.py)

def calculate_origin_offset(new_spacing, old_spacing):
return np.divide(np.subtract(new_spacing, old_spacing)/2, old_spacing)

Another question is leveraging other libraries like skimage and torchio. How resize/transform in skimage/ torchio approach is different from scaling affine first and then using resample_img from nilearn

https://stackoverflow.com/questions/64674612/how-to-resize-a-nifti-nii-gz-medical-image-file