Source code for spinalcordtoolbox.resampling

#########################################################################################
#
# Resample data using nibabel.
#
# ---------------------------------------------------------------------------------------
# Copyright (c) 2014 Polytechnique Montreal <www.neuro.polymtl.ca>
# Authors: Julien Cohen-Adad, Sara Dupont
#
# About the license: see the file LICENSE.TXT
#########################################################################################

# TODO: remove resample_file (not needed)

import logging

import numpy as np
import nibabel as nib
from nibabel.processing import resample_from_to

from spinalcordtoolbox.image import Image, add_suffix
from spinalcordtoolbox.utils import display_viewer_syntax

logger = logging.getLogger(__name__)


[docs]def resample_nib(image, new_size=None, new_size_type=None, image_dest=None, interpolation='linear', mode='nearest'): """ Resample a nibabel or Image object based on a specified resampling factor. Can deal with 2d, 3d or 4d image objects. :param image: nibabel or Image image. :param new_size: list of float: Resampling factor, final dimension or resolution, depending on new_size_type. :param new_size_type: {'vox', 'factor', 'mm'}: Feature used for resampling. Examples: new_size=[128, 128, 90], new_size_type='vox' --> Resampling to a dimension of 128x128x90 voxels new_size=[2, 2, 2], new_size_type='factor' --> 2x isotropic upsampling new_size=[1, 1, 5], new_size_type='mm' --> Resampling to a resolution of 1x1x5 mm :param image_dest: Destination image to resample the input image to. In this case, new_size and new_size_type are ignored :param interpolation: {'nn', 'linear', 'spline'}. The interpolation type :param mode: Outside values are filled with 0 ('constant') or nearest value ('nearest'). :return: The resampled nibabel or Image image (depending on the input object type). """ # set interpolation method dict_interp = {'nn': 0, 'linear': 1, 'spline': 2} # If input is an Image object, create nibabel object from it if type(image) == nib.nifti1.Nifti1Image: img = image elif type(image) == Image: img = nib.nifti1.Nifti1Image(image.data, image.hdr.get_best_affine()) else: raise Exception(TypeError) if image_dest is None: # Get dimensions of data p = img.header.get_zooms() shape = img.header.get_data_shape() if img.ndim == 4: new_size += ['1'] # needed because the code below is general, i.e., does not assume 3d input and uses img.shape # compute new shape based on specific resampling method if new_size_type == 'vox': shape_r = tuple([int(new_size[i]) for i in range(img.ndim)]) elif new_size_type == 'factor': if len(new_size) == 1: # isotropic resampling new_size = tuple([new_size[0] for i in range(img.ndim)]) # compute new shape as: shape_r = shape * f shape_r = tuple([int(np.round(shape[i] * float(new_size[i]))) for i in range(img.ndim)]) elif new_size_type == 'mm': if len(new_size) == 1: # isotropic resampling new_size = tuple([new_size[0] for i in range(img.ndim)]) # compute new shape as: shape_r = shape * (p_r / p) shape_r = tuple([int(np.round(shape[i] * float(p[i]) / float(new_size[i]))) for i in range(img.ndim)]) else: raise ValueError("'new_size_type' is not recognized.") # Generate 3d affine transformation: R affine = img.affine[:4, :4] affine[3, :] = np.array([0, 0, 0, 1]) # satisfy to nifti convention. Otherwise it grabs the temporal logger.debug('Affine matrix: \n' + str(affine)) R = np.eye(4) for i in range(3): try: R[i, i] = img.shape[i] / float(shape_r[i]) except ZeroDivisionError: raise ZeroDivisionError("Destination size is zero for dimension {}. You are trying to resample to an " "unrealistic dimension. Check your NIFTI pixdim values to make sure they are " "not corrupted.".format(i)) affine_r = np.dot(affine, R) reference = (shape_r, affine_r) # If reference is provided else: if type(image_dest) == nib.nifti1.Nifti1Image: reference = image_dest elif type(image_dest) == Image: reference = nib.nifti1.Nifti1Image(image_dest.data, image_dest.hdr.get_best_affine()) else: raise Exception(TypeError) if img.ndim == 3: # we use mode 'nearest' to overcome issue #2453 img_r = resample_from_to( img, to_vox_map=reference, order=dict_interp[interpolation], mode=mode, cval=0.0, out_class=None) elif img.ndim == 4: # TODO: Cover img_dest with 4D volumes # Import here instead of top of the file because this is an isolated case and nibabel takes time to import data4d = np.zeros(shape_r) # Loop across 4th dimension and resample each 3d volume for it in range(img.shape[3]): # Create dummy 3d nibabel image nii_tmp = nib.nifti1.Nifti1Image(img.get_data()[..., it], affine) img3d_r = resample_from_to( nii_tmp, to_vox_map=(shape_r[:-1], affine_r), order=dict_interp[interpolation], mode=mode, cval=0.0, out_class=None) data4d[..., it] = img3d_r.get_data() # Create 4d nibabel Image img_r = nib.nifti1.Nifti1Image(data4d, affine_r) # Convert back to proper type if type(image) == nib.nifti1.Nifti1Image: return img_r elif type(image) == Image: return Image(img_r.get_data(), hdr=img_r.header, orientation=image.orientation, dim=img_r.header.get_data_shape())
[docs]def resample_file(fname_data, fname_out, new_size, new_size_type, interpolation, verbose, fname_ref=None): """This function will resample the specified input image file to the target size. Can deal with 2d, 3d or 4d image objects. :param fname_data: The input image filename. :param fname_out: The output image filename. :param new_size: The target size, i.e. 0.25x0.25 :param new_size_type: Unit of resample (mm, vox, factor) :param interpolation: The interpolation type :param verbose: verbosity level :param fname_ref: Reference image to resample input image to """ # Load data logger.info('load data...') nii = nib.load(fname_data) if fname_ref is not None: nii_ref = nib.load(fname_ref) else: nii_ref = None nii_r = resample_nib(nii, new_size.split('x'), new_size_type, image_dest=nii_ref, interpolation=interpolation) # build output file name if fname_out == '': fname_out = add_suffix(fname_data, '_r') else: fname_out = fname_out # save data nib.save(nii_r, fname_out) # to view results display_viewer_syntax([fname_out], verbose=verbose) return nii_r