diff --git a/nibabel/__init__.py b/nibabel/__init__.py index 2e4f877c5f..948cc9247f 100644 --- a/nibabel/__init__.py +++ b/nibabel/__init__.py @@ -65,6 +65,8 @@ def teardown_package(): from . import spm2analyze as spm2 from . import nifti1 as ni1 from . import ecat +from . import transform + # object imports from .fileholders import FileHolder, FileHolderError from .loadsave import load, save diff --git a/nibabel/affines.py b/nibabel/affines.py index a7d7a4e9b8..b09257c7ad 100644 --- a/nibabel/affines.py +++ b/nibabel/affines.py @@ -2,7 +2,6 @@ # vi: set ft=python sts=4 ts=4 sw=4 et: """ Utility routines for working with points and affine transforms """ - import numpy as np from functools import reduce @@ -296,3 +295,31 @@ def voxel_sizes(affine): """ top_left = affine[:-1, :-1] return np.sqrt(np.sum(top_left ** 2, axis=0)) + + +def obliquity(affine): + r""" + Estimate the *obliquity* an affine's axes represent. + + The term *obliquity* is defined here as the rotation of those axes with + respect to the cardinal axes. + This implementation is inspired by `AFNI's implementation + `_. + For further details about *obliquity*, check `AFNI's documentation + _. + + Parameters + ---------- + affine : 2D array-like + Affine transformation array. Usually shape (4, 4), but can be any 2D + array. + + Returns + ------- + angles : 1D array-like + The *obliquity* of each axis with respect to the cardinal axes, in radians. + + """ + vs = voxel_sizes(affine) + best_cosines = np.abs((affine[:-1, :-1] / vs).max(axis=1)) + return np.arccos(best_cosines) diff --git a/nibabel/tests/data/affine-LAS-itk.tfm b/nibabel/tests/data/affine-LAS-itk.tfm new file mode 120000 index 0000000000..849a5e1667 --- /dev/null +++ b/nibabel/tests/data/affine-LAS-itk.tfm @@ -0,0 +1 @@ +affine-RAS-itk.tfm \ No newline at end of file diff --git a/nibabel/tests/data/affine-LAS.afni b/nibabel/tests/data/affine-LAS.afni new file mode 120000 index 0000000000..1bcb645d8f --- /dev/null +++ b/nibabel/tests/data/affine-LAS.afni @@ -0,0 +1 @@ +affine-RAS.afni \ No newline at end of file diff --git a/nibabel/tests/data/affine-LAS.fsl b/nibabel/tests/data/affine-LAS.fsl new file mode 120000 index 0000000000..aefb123fd3 --- /dev/null +++ b/nibabel/tests/data/affine-LAS.fsl @@ -0,0 +1 @@ +affine-RAS.fsl \ No newline at end of file diff --git a/nibabel/tests/data/affine-LPS-itk.tfm b/nibabel/tests/data/affine-LPS-itk.tfm new file mode 100644 index 0000000000..b31606e275 --- /dev/null +++ b/nibabel/tests/data/affine-LPS-itk.tfm @@ -0,0 +1,5 @@ +#Insight Transform File V1.0 +#Transform 0 +Transform: MatrixOffsetTransformBase_double_3_3 +Parameters: 0.999999 -0.000999999 -0.001 0.00140494 0.621609 0.783327 -0.000161717 -0.783327 0.62161 -4 -2 -1 +FixedParameters: 0 0 0 diff --git a/nibabel/tests/data/affine-LPS.afni b/nibabel/tests/data/affine-LPS.afni new file mode 120000 index 0000000000..1bcb645d8f --- /dev/null +++ b/nibabel/tests/data/affine-LPS.afni @@ -0,0 +1 @@ +affine-RAS.afni \ No newline at end of file diff --git a/nibabel/tests/data/affine-LPS.fsl b/nibabel/tests/data/affine-LPS.fsl new file mode 100644 index 0000000000..baf012ae22 --- /dev/null +++ b/nibabel/tests/data/affine-LPS.fsl @@ -0,0 +1,4 @@ +0.999999 -0.00140494 0.000161717 -3.89014 +0.000999999 0.621609 -0.783327 105.905 +0.001 0.783327 0.62161 -34.3513 +0 0 0 1 diff --git a/nibabel/tests/data/affine-RAS-itk.tfm b/nibabel/tests/data/affine-RAS-itk.tfm new file mode 100644 index 0000000000..b31606e275 --- /dev/null +++ b/nibabel/tests/data/affine-RAS-itk.tfm @@ -0,0 +1,5 @@ +#Insight Transform File V1.0 +#Transform 0 +Transform: MatrixOffsetTransformBase_double_3_3 +Parameters: 0.999999 -0.000999999 -0.001 0.00140494 0.621609 0.783327 -0.000161717 -0.783327 0.62161 -4 -2 -1 +FixedParameters: 0 0 0 diff --git a/nibabel/tests/data/affine-RAS.afni b/nibabel/tests/data/affine-RAS.afni new file mode 100644 index 0000000000..3239906695 --- /dev/null +++ b/nibabel/tests/data/affine-RAS.afni @@ -0,0 +1,2 @@ +# 3dvolreg matrices (DICOM-to-DICOM, row-by-row): +0.999999 -0.000999999 -0.001 -4 0.00140494 0.621609 0.783327 -2 -0.000161717 -0.783327 0.62161 -1 diff --git a/nibabel/tests/data/affine-RAS.fsl b/nibabel/tests/data/affine-RAS.fsl new file mode 100644 index 0000000000..d26cf54ba0 --- /dev/null +++ b/nibabel/tests/data/affine-RAS.fsl @@ -0,0 +1,4 @@ +0.999999 -0.00140494 -0.000161717 4.14529 +0.000999999 0.621609 0.783327 -37.3811 +-0.001 -0.783327 0.62161 107.976 +0 0 0 1 diff --git a/nibabel/tests/data/affine-oblique-itk.tfm b/nibabel/tests/data/affine-oblique-itk.tfm new file mode 100644 index 0000000000..b31606e275 --- /dev/null +++ b/nibabel/tests/data/affine-oblique-itk.tfm @@ -0,0 +1,5 @@ +#Insight Transform File V1.0 +#Transform 0 +Transform: MatrixOffsetTransformBase_double_3_3 +Parameters: 0.999999 -0.000999999 -0.001 0.00140494 0.621609 0.783327 -0.000161717 -0.783327 0.62161 -4 -2 -1 +FixedParameters: 0 0 0 diff --git a/nibabel/tests/data/affine-oblique.fsl b/nibabel/tests/data/affine-oblique.fsl new file mode 100644 index 0000000000..48f811e63b --- /dev/null +++ b/nibabel/tests/data/affine-oblique.fsl @@ -0,0 +1,4 @@ +0.999998 -0.00181872 -0.0011965 4.26083 +0.00206779 0.621609 0.783325 -25.3129 +-0.000680894 -0.783326 0.621611 101.967 +0 0 0 1 diff --git a/nibabel/tests/data/someones_anatomy.nii.gz b/nibabel/tests/data/someones_anatomy.nii.gz new file mode 120000 index 0000000000..b1afffc190 --- /dev/null +++ b/nibabel/tests/data/someones_anatomy.nii.gz @@ -0,0 +1 @@ +../../../doc/source/downloads/someones_anatomy.nii.gz \ No newline at end of file diff --git a/nibabel/tests/test_affines.py b/nibabel/tests/test_affines.py index e66ed46190..13b554c5a8 100644 --- a/nibabel/tests/test_affines.py +++ b/nibabel/tests/test_affines.py @@ -7,7 +7,7 @@ from ..eulerangles import euler2mat from ..affines import (AffineError, apply_affine, append_diag, to_matvec, - from_matvec, dot_reduce, voxel_sizes) + from_matvec, dot_reduce, voxel_sizes, obliquity) from nose.tools import assert_equal, assert_raises @@ -178,3 +178,15 @@ def test_voxel_sizes(): rot_affine[:3, :3] = rotation full_aff = rot_affine.dot(aff) assert_almost_equal(voxel_sizes(full_aff), vox_sizes) + + +def test_obliquity(): + """Check the calculation of inclination of an affine axes.""" + from math import pi + aligned = np.diag([2.0, 2.0, 2.3, 1.0]) + aligned[:-1, -1] = [-10, -10, -7] + R = from_matvec(euler2mat(x=0.09, y=0.001, z=0.001), [0.0, 0.0, 0.0]) + oblique = R.dot(aligned) + assert_almost_equal(obliquity(aligned), [0.0, 0.0, 0.0]) + assert_almost_equal(obliquity(oblique) * 180 / pi, + [0.0810285, 5.1569949, 5.1569376]) diff --git a/nibabel/tests/test_transform.py b/nibabel/tests/test_transform.py new file mode 100644 index 0000000000..f8160e2ee1 --- /dev/null +++ b/nibabel/tests/test_transform.py @@ -0,0 +1,75 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +"""Tests of the transform module.""" +import os +import numpy as np +from numpy.testing import assert_array_equal, assert_almost_equal, \ + assert_array_almost_equal +import pytest + +from ..loadsave import load as loadimg +from ..nifti1 import Nifti1Image +from ..eulerangles import euler2mat +from ..affines import from_matvec +from ..volumeutils import shape_zoom_affine +from ..transform import linear as nbl +from ..testing import (assert_equal, assert_not_equal, assert_true, + assert_false, assert_raises, data_path, + suppress_warnings, assert_dt_equal) +from ..tmpdirs import InTemporaryDirectory + + +SOMEONES_ANATOMY = os.path.join(data_path, 'someones_anatomy.nii.gz') +# SOMEONES_ANATOMY = os.path.join(data_path, 'someones_anatomy.nii.gz') + + +@pytest.mark.parametrize('image_orientation', ['RAS', 'LAS', 'LPS', 'oblique']) +def test_affines_save(image_orientation): + """Check implementation of exporting affines to formats.""" + # Generate test transform + img = loadimg(SOMEONES_ANATOMY) + imgaff = img.affine + + if image_orientation == 'LAS': + newaff = imgaff.copy() + newaff[0, 0] *= -1.0 + newaff[0, 3] = imgaff.dot(np.hstack((np.array(img.shape[:3]) - 1, 1.0)))[0] + img = Nifti1Image(np.flip(img.get_fdata(), 0), newaff, img.header) + elif image_orientation == 'LPS': + newaff = imgaff.copy() + newaff[0, 0] *= -1.0 + newaff[1, 1] *= -1.0 + newaff[:2, 3] = imgaff.dot(np.hstack((np.array(img.shape[:3]) - 1, 1.0)))[:2] + img = Nifti1Image(np.flip(np.flip(img.get_fdata(), 0), 1), newaff, img.header) + elif image_orientation == 'oblique': + A = shape_zoom_affine(img.shape, img.header.get_zooms(), x_flip=False) + R = from_matvec(euler2mat(x=0.09, y=0.001, z=0.001)) + newaff = R.dot(A) + img = Nifti1Image(img.get_fdata(), newaff, img.header) + img.header.set_qform(newaff, 1) + img.header.set_sform(newaff, 1) + + T = from_matvec(euler2mat(x=0.9, y=0.001, z=0.001), [4.0, 2.0, -1.0]) + + xfm = nbl.Affine(T) + xfm.reference = img + + itk = nbl.load(os.path.join(data_path, 'affine-%s-itk.tfm' % image_orientation), + fmt='itk') + fsl = np.loadtxt(os.path.join(data_path, 'affine-%s.fsl' % image_orientation)) + afni = np.loadtxt(os.path.join(data_path, 'affine-%s.afni' % image_orientation)) + + with InTemporaryDirectory(): + xfm.to_filename('M.tfm', fmt='itk') + xfm.to_filename('M.fsl', fmt='fsl') + xfm.to_filename('M.afni', fmt='afni') + + nb_itk = nbl.load('M.tfm', fmt='itk') + nb_fsl = np.loadtxt('M.fsl') + nb_afni = np.loadtxt('M.afni') + + assert_equal(itk, nb_itk) + assert_almost_equal(fsl, nb_fsl) + assert_almost_equal(afni, nb_afni) + +# Create version not aligned to canonical diff --git a/nibabel/transform/__init__.py b/nibabel/transform/__init__.py new file mode 100644 index 0000000000..be52d50c42 --- /dev/null +++ b/nibabel/transform/__init__.py @@ -0,0 +1,23 @@ +# emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## +# +# See COPYING file distributed along with the NiBabel package for the +# copyright and license terms. +# +### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## +""" +Geometric transforms. + +.. currentmodule:: nibabel.transform + +.. autosummary:: + :toctree: ../generated + + transform +""" +from .linear import Affine +from .nonlinear import DeformationFieldTransform + + +__all__ = ['Affine', 'DeformationFieldTransform'] diff --git a/nibabel/transform/base.py b/nibabel/transform/base.py new file mode 100644 index 0000000000..ec9d4e8b6c --- /dev/null +++ b/nibabel/transform/base.py @@ -0,0 +1,199 @@ +# emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## +# +# See COPYING file distributed along with the NiBabel package for the +# copyright and license terms. +# +### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## +"""Common interface for transforms.""" +import numpy as np +import h5py + +from scipy import ndimage as ndi + +EQUALITY_TOL = 1e-5 + + +class ImageSpace(object): + """Class to represent spaces of gridded data (images).""" + + __slots__ = ['_affine', '_shape', '_ndim', '_ndindex', '_coords', '_nvox', + '_inverse'] + + def __init__(self, image): + self._affine = image.affine + self._shape = image.shape + self._ndim = len(image.shape) + self._nvox = np.prod(image.shape) # Do not access data array + self._ndindex = None + self._coords = None + self._inverse = np.linalg.inv(image.affine) + if self._ndim not in [2, 3]: + raise ValueError('Invalid image space (%d-D)' % self._ndim) + + @property + def affine(self): + return self._affine + + @property + def inverse(self): + return self._inverse + + @property + def shape(self): + return self._shape + + @property + def ndim(self): + return self._ndim + + @property + def nvox(self): + return self._nvox + + @property + def ndindex(self): + if self._ndindex is None: + indexes = tuple([np.arange(s) for s in self._shape]) + self._ndindex = np.array(np.meshgrid( + *indexes, indexing='ij')).reshape(self._ndim, self._nvox) + return self._ndindex + + @property + def ndcoords(self): + if self._coords is None: + self._coords = np.tensordot( + self._affine, + np.vstack((self.ndindex, np.ones((1, self._nvox)))), + axes=1 + )[:3, ...] + return self._coords + + def map_voxels(self, coordinates): + coordinates = np.array(coordinates) + ncoords = coordinates.shape[-1] + coordinates = np.vstack((coordinates, np.ones((1, ncoords)))) + + # Back to grid coordinates + return np.tensordot(np.linalg.inv(self._affine), + coordinates, axes=1)[:3, ...] + + def _to_hdf5(self, group): + group.attrs['Type'] = 'image' + group.attrs['ndim'] = self.ndim + group.create_dataset('affine', data=self.affine) + group.create_dataset('shape', data=self.shape) + + def __eq__(self, other): + try: + return np.allclose( + self.affine, other.affine, rtol=EQUALITY_TOL) and self.shape == other.shape + except AttributeError: + pass + return False + + +class TransformBase(object): + """Abstract image class to represent transforms.""" + + __slots__ = ['_reference'] + + def __init__(self): + self._reference = None + + def __eq__(self, other): + """Overload equals operator.""" + if not self._reference == other._reference: + return False + return np.allclose(self.matrix, other.matrix, rtol=EQUALITY_TOL) + + @property + def reference(self): + '''A reference space where data will be resampled onto''' + if self._reference is None: + raise ValueError('Reference space not set') + return self._reference + + @reference.setter + def reference(self, image): + self._reference = ImageSpace(image) + + @property + def ndim(self): + return self.reference.ndim + + def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True, + output_dtype=None): + """ + Resample the moving image in reference space. + + Parameters + ---------- + moving : `spatialimage` + The image object containing the data to be resampled in reference + space + order : int, optional + The order of the spline interpolation, default is 3. + The order has to be in the range 0-5. + mode : {'reflect', 'constant', 'nearest', 'mirror', 'wrap'}, optional + Determines how the input image is extended when the resamplings overflows + a border. Default is 'constant'. + cval : float, optional + Constant value for ``mode='constant'``. Default is 0.0. + prefilter: bool, optional + Determines if the moving image's data array is prefiltered with + a spline filter before interpolation. The default is ``True``, + which will create a temporary *float64* array of filtered values + if *order > 1*. If setting this to ``False``, the output will be + slightly blurred if *order > 1*, unless the input is prefiltered, + i.e. it is the result of calling the spline filter on the original + input. + + Returns + ------- + moved_image : `spatialimage` + The moving imaged after resampling to reference space. + + """ + moving_data = np.asanyarray(moving.dataobj) + if output_dtype is None: + output_dtype = moving_data.dtype + + moved = ndi.geometric_transform( + moving_data, + mapping=self.map_voxel, + output_shape=self.reference.shape, + output=output_dtype, + order=order, + mode=mode, + cval=cval, + prefilter=prefilter, + extra_keywords={'moving': moving}, + ) + + moved_image = moving.__class__(moved, self.reference.affine, moving.header) + moved_image.header.set_data_dtype(output_dtype) + return moved_image + + def map_point(self, coords): + """Apply y = f(x), where x is the argument `coords`.""" + raise NotImplementedError + + def map_voxel(self, index, moving=None): + """Apply ijk' = f_ijk((i, j, k)), equivalent to the above with indexes.""" + raise NotImplementedError + + def _to_hdf5(self, x5_root): + """Serialize this object into the x5 file format.""" + raise NotImplementedError + + def to_filename(self, filename, fmt='X5'): + """Store the transform in BIDS-Transforms HDF5 file format (.x5).""" + with h5py.File(filename, 'w') as out_file: + out_file.attrs['Format'] = 'X5' + out_file.attrs['Version'] = np.uint16(1) + root = out_file.create_group('/0') + self._to_hdf5(root) + + return filename diff --git a/nibabel/transform/linear.py b/nibabel/transform/linear.py new file mode 100644 index 0000000000..84017b203a --- /dev/null +++ b/nibabel/transform/linear.py @@ -0,0 +1,335 @@ +# emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## +# +# See COPYING file distributed along with the NiBabel package for the +# copyright and license terms. +# +### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## +"""Linear transforms.""" +import sys +import numpy as np +from scipy import ndimage as ndi +from pathlib import Path + +from ..loadsave import load as loadimg +from ..affines import from_matvec, voxel_sizes, obliquity +from ..volumeutils import shape_zoom_affine +from .base import TransformBase + + +LPS = np.diag([-1, -1, 1, 1]) +OBLIQUITY_THRESHOLD_DEG = 0.01 + + +class Affine(TransformBase): + """Represents linear transforms on image data.""" + + __slots__ = ['_matrix'] + + def __init__(self, matrix=None, reference=None): + """ + Initialize a linear transform. + + Parameters + ---------- + matrix : ndarray + The inverse coordinate transformation matrix **in physical + coordinates**, mapping coordinates from *reference* space + into *moving* space. + This matrix should be provided in homogeneous coordinates. + + Examples + -------- + >>> xfm = Affine([[1, 0, 0, 4], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]) + >>> xfm.matrix # doctest: +NORMALIZE_WHITESPACE + array([[1, 0, 0, 4], + [0, 1, 0, 0], + [0, 0, 1, 0], + [0, 0, 0, 1]]) + + """ + super(Affine, self).__init__() + if matrix is None: + matrix = [np.eye(4)] + + if np.ndim(matrix) == 2: + matrix = [matrix] + + self._matrix = np.array(matrix) + assert self._matrix.ndim == 3, 'affine matrix should be 3D' + assert self._matrix.shape[-2] == self._matrix.shape[-1], 'affine matrix is not square' + + if reference: + if isinstance(reference, str): + reference = loadimg(reference) + self.reference = reference + + @property + def matrix(self): + """Access the internal representation of this affine.""" + return self._matrix + + def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True, + output_dtype=None): + """ + Resample the moving image in reference space. + + Parameters + ---------- + moving : `spatialimage` + The image object containing the data to be resampled in reference + space + order : int, optional + The order of the spline interpolation, default is 3. + The order has to be in the range 0-5. + mode : {'reflect', 'constant', 'nearest', 'mirror', 'wrap'}, optional + Determines how the input image is extended when the resamplings overflows + a border. Default is 'constant'. + cval : float, optional + Constant value for ``mode='constant'``. Default is 0.0. + prefilter: bool, optional + Determines if the moving image's data array is prefiltered with + a spline filter before interpolation. The default is ``True``, + which will create a temporary *float64* array of filtered values + if *order > 1*. If setting this to ``False``, the output will be + slightly blurred if *order > 1*, unless the input is prefiltered, + i.e. it is the result of calling the spline filter on the original + input. + + Returns + ------- + moved_image : `spatialimage` + The moving imaged after resampling to reference space. + + + Examples + -------- + >>> import nibabel as nib + >>> xfm = Affine([[1, 0, 0, 4], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]) + >>> ref = nib.load('image.nii.gz') + >>> xfm.reference = ref + >>> xfm.resample(ref, order=0) + + """ + if output_dtype is None: + output_dtype = moving.header.get_data_dtype() + + try: + reference = self.reference + except ValueError: + print('Warning: no reference space defined, using moving as reference', + file=sys.stderr) + reference = moving + + nvols = 1 + if len(moving.shape) > 3: + nvols = moving.shape[3] + + movaff = moving.affine + movingdata = moving.get_data() + if nvols == 1: + movingdata = movingdata[..., np.newaxis] + + nmats = self._matrix.shape[0] + if nvols != nmats and nmats > 1: + raise ValueError("""\ +The moving image contains {0} volumes, while the transform is defined for \ +{1} volumes""".format(nvols, nmats)) + + singlemat = None + if nmats == 1: + singlemat = np.linalg.inv(movaff).dot(self._matrix[0].dot(reference.affine)) + + if singlemat is not None and nvols > nmats: + print('Warning: resampling a 4D volume with a single affine matrix', + file=sys.stderr) + + # Compose an index to index affine matrix + moved = [] + for i in range(nvols): + i2imat = singlemat if singlemat is not None else np.linalg.inv( + movaff).dot(self._matrix[i].dot(reference.affine)) + mdata = movingdata[..., i] + moved += [ndi.affine_transform( + mdata, matrix=i2imat[:-1, :-1], + offset=i2imat[:-1, -1], + output_shape=reference.shape, order=order, mode=mode, + cval=cval, prefilter=prefilter)] + + moved_image = moving.__class__( + np.squeeze(np.stack(moved, -1)), reference.affine, moving.header) + moved_image.header.set_data_dtype(output_dtype) + + return moved_image + + def map_point(self, coords, index=0, forward=True): + """Apply y = f(x), where x is the argument `coords`.""" + coords = np.array(coords) + if coords.shape[0] == self._matrix[index].shape[0] - 1: + coords = np.append(coords, [1]) + affine = self._matrix[index] if forward else np.linalg.inv(self._matrix[index]) + return affine.dot(coords)[:-1] + + def map_voxel(self, index, nindex=0, moving=None): + """Apply ijk' = f_ijk((i, j, k)), equivalent to the above with indexes.""" + try: + reference = self.reference + except ValueError: + print('Warning: no reference space defined, using moving as reference', + file=sys.stderr) + reference = moving + else: + if moving is None: + moving = reference + finally: + if reference is None: + raise ValueError('Reference and moving spaces are both undefined') + + index = np.array(index) + if index.shape[0] == self._matrix[nindex].shape[0] - 1: + index = np.append(index, [1]) + + matrix = reference.affine.dot(self._matrix[nindex].dot(np.linalg.inv(moving.affine))) + return tuple(matrix.dot(index)[:-1]) + + def _to_hdf5(self, x5_root): + """Serialize this object into the x5 file format.""" + xform = x5_root.create_dataset('Transform', data=self._matrix) + xform.attrs['Type'] = 'affine' + x5_root.create_dataset('Inverse', data=np.linalg.inv(self._matrix)) + + if self._reference: + self.reference._to_hdf5(x5_root.create_group('Reference')) + + def to_filename(self, filename, fmt='X5', moving=None): + """Store the transform in BIDS-Transforms HDF5 file format (.x5).""" + if fmt.lower() in ['itk', 'ants', 'elastix']: + with open(filename, 'w') as f: + f.write('#Insight Transform File V1.0\n') + + for i in range(self.matrix.shape[0]): + parameters = LPS.dot(self.matrix[i].dot(LPS)) + parameters = parameters[:3, :3].reshape(-1).tolist() + \ + parameters[:3, 3].tolist() + itkfmt = """\ +#Transform {0} +Transform: MatrixOffsetTransformBase_double_3_3 +Parameters: {1} +FixedParameters: 0 0 0\n""".format + f.write(itkfmt(i, ' '.join(['%g' % p for p in parameters]))) + return filename + + if fmt.lower() == 'afni': + from math import pi + + if moving and isinstance(moving, (str, bytes, Path)): + moving = loadimg(str(moving)) + + T = self.matrix.copy() + pre = LPS + post = LPS + if obliquity(self.reference.affine).min() * 180 / pi > OBLIQUITY_THRESHOLD_DEG: + print('Reference affine axes are oblique.') + M = self.reference.affine + A = shape_zoom_affine(self.reference.shape, + voxel_sizes(M), x_flip=False, y_flip=False) + pre = M.dot(np.linalg.inv(A)).dot(LPS) + + if not moving: + moving = self.reference + + if moving and obliquity(moving.affine).min() * 180 / pi > OBLIQUITY_THRESHOLD_DEG: + print('Moving affine axes are oblique.') + M2 = moving.affine + A2 = shape_zoom_affine(moving.shape, + voxel_sizes(M2), x_flip=True, y_flip=True) + post = A2.dot(np.linalg.inv(M2)) + + # swapaxes is necessary, as axis 0 encodes series of transforms + parameters = np.swapaxes(post.dot(self.matrix.copy().dot(pre)), 0, 1) + parameters = parameters[:, :3, :].reshape((T.shape[0], -1)) + np.savetxt(filename, parameters, delimiter='\t', header="""\ +3dvolreg matrices (DICOM-to-DICOM, row-by-row):""", fmt='%g') + return filename + + if fmt.lower() == 'fsl': + if not moving: + moving = self.reference + + if isinstance(moving, str): + moving = loadimg(moving) + + # Adjust for reference image offset and orientation + refswp, refspc = _fsl_aff_adapt(self.reference) + pre = self.reference.affine.dot( + np.linalg.inv(refspc).dot(np.linalg.inv(refswp))) + + # Adjust for moving image offset and orientation + movswp, movspc = _fsl_aff_adapt(moving) + post = np.linalg.inv(movswp).dot(movspc.dot(np.linalg.inv( + moving.affine))) + + # Compose FSL transform + mat = np.linalg.inv( + np.swapaxes(post.dot(self.matrix.dot(pre)), 0, 1)) + + if self.matrix.shape[0] > 1: + for i in range(self.matrix.shape[0]): + np.savetxt('%s.%03d' % (filename, i), mat[i], + delimiter=' ', fmt='%g') + else: + np.savetxt(filename, mat[0], delimiter=' ', fmt='%g') + return filename + return super(Affine, self).to_filename(filename, fmt=fmt) + + +def load(filename, fmt='X5', reference=None): + """Load a linear transform.""" + if fmt.lower() in ['itk', 'ants', 'elastix', 'nifty']: + with open(filename) as itkfile: + itkxfm = itkfile.read().splitlines() + + if 'Insight Transform File' not in itkxfm[0]: + raise ValueError( + "File %s does not look like an ITK transform text file." % filename) + + matlist = [] + for nxfm in range(len(itkxfm[1:]) // 4): + parameters = np.fromstring(itkxfm[nxfm * 4 + 3].split(':')[-1].strip(), + dtype=float, sep=' ') + offset = np.fromstring(itkxfm[nxfm * 4 + 4].split(':')[-1].strip(), + dtype=float, sep=' ') + if len(parameters) == 12: + matrix = from_matvec(parameters[:9].reshape((3, 3)), parameters[9:]) + c_neg = from_matvec(np.eye(3), offset * -1.0) + c_pos = from_matvec(np.eye(3), offset) + matlist.append(LPS.dot(c_pos.dot(matrix.dot(c_neg.dot(LPS))))) + matrix = np.stack(matlist) + # elif fmt.lower() == 'afni': + # parameters = LPS.dot(self.matrix.dot(LPS)) + # parameters = parameters[:3, :].reshape(-1).tolist() + elif fmt.lower() in ('x5', 'bids'): + raise NotImplementedError + else: + raise NotImplementedError + + if reference and isinstance(reference, str): + reference = loadimg(reference) + return Affine(matrix, reference) + + +def _fsl_aff_adapt(space): + """ + Adapt FSL affines. + + Calculates a matrix to convert from the original RAS image + coordinates to FSL's internal coordinate system of transforms + """ + aff = space.affine + zooms = list(voxel_sizes(aff)) + [1] + swp = np.eye(4) + if np.linalg.det(aff) > 0: + swp[0, 0] = -1.0 + swp[0, 3] = (space.shape[0] - 1) * zooms[0] + return swp, np.diag(zooms) diff --git a/nibabel/transform/nonlinear.py b/nibabel/transform/nonlinear.py new file mode 100644 index 0000000000..e607bd144f --- /dev/null +++ b/nibabel/transform/nonlinear.py @@ -0,0 +1,250 @@ +# emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## +# +# See COPYING file distributed along with the NiBabel package for the +# copyright and license terms. +# +### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## +"""Nonlinear transforms.""" +import numpy as np +from scipy import ndimage as ndi +# from gridbspline.maths import cubic + +from .base import ImageSpace, TransformBase +from ..funcs import four_to_three + +# vbspl = np.vectorize(cubic) + + +class DeformationFieldTransform(TransformBase): + """Represents a dense field of displacements (one vector per voxel).""" + + __slots__ = ['_field', '_moving', '_moving_space'] + __s = (slice(None), ) + + def __init__(self, field, reference=None): + """Create a dense deformation field transform.""" + super(DeformationFieldTransform, self).__init__() + self._field = field.get_data() + + ndim = self._field.ndim - 1 + if self._field.shape[:-1] != ndim: + raise ValueError( + 'Number of components of the deformation field does ' + 'not match the number of dimensions') + + self._moving = None # Where each voxel maps to + self._moving_space = None # Input space cache + + # By definition, a free deformation field has a + # displacement vector per voxel in output (reference) + # space + if reference is None: + reference = four_to_three(field)[0] + elif reference.shape[:ndim] != field.shape[:ndim]: + raise ValueError( + 'Reference ({}) and field ({}) must have the same ' + 'grid.'.format( + _pprint(reference.shape[:ndim]), + _pprint(field.shape[:ndim]))) + + self.reference = reference + + def _cache_moving(self, moving): + # Check whether input (moving) space is cached + moving_space = ImageSpace(moving) + if self._moving_space == moving_space: + return + + # Generate grid of pixel indexes (ijk) + ndim = self._field.ndim - 1 + if ndim == 2: + grid = np.meshgrid( + np.arange(self._field.shape[0]), + np.arange(self._field.shape[1]), + indexing='ij') + elif ndim == 3: + grid = np.meshgrid( + np.arange(self._field.shape[0]), + np.arange(self._field.shape[1]), + np.arange(self._field.shape[2]), + indexing='ij') + else: + raise ValueError('Wrong dimensions (%d)' % ndim) + + grid = np.array(grid) + flatgrid = grid.reshape(ndim, -1) + + # Calculate physical coords of all voxels (xyz) + flatxyz = np.tensordot( + self.reference.affine, + np.vstack((flatgrid, np.ones((1, flatgrid.shape[1])))), + axes=1 + ) + + # Add field + newxyz = flatxyz + np.vstack(( + np.moveaxis(self._field, -1, 0).reshape(ndim, -1), + np.zeros((1, flatgrid.shape[1])))) + + # Back to grid coordinates + newijk = np.tensordot(np.linalg.inv(moving.affine), + newxyz, axes=1) + + # Reshape as grid + self._moving = np.moveaxis( + newijk[0:3, :].reshape((ndim, ) + self._field.shape[:-1]), + 0, -1) + + self._moving_space = moving_space + + def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True, + output_dtype=None): + """ + Resample the `moving` image applying the deformation field. + + Examples + -------- + >>> import numpy as np + >>> import nibabel as nb + >>> ref = nb.load('t1_weighted.nii.gz') + >>> field = np.zeros(tuple(list(ref.shape) + [3])) + >>> field[..., 0] = 4.0 + >>> fieldimg = nb.Nifti1Image(field, ref.affine, ref.header) + >>> xfm = nb.transform.DeformationFieldTransform(fieldimg) + >>> new = xfm.resample(ref) + >>> new.to_filename('deffield.nii.gz') + + """ + self._cache_moving(moving) + return super(DeformationFieldTransform, self).resample( + moving, order=order, mode=mode, cval=cval, prefilter=prefilter) + + def map_voxel(self, index, moving=None): + """Apply ijk' = f_ijk((i, j, k)), equivalent to the above with indexes.""" + return tuple(self._moving[index + self.__s]) + + def map_coordinates(self, coordinates, order=3, mode='mirror', cval=0.0, + prefilter=True): + """Apply y = f(x), where x is the argument `coords`.""" + coordinates = np.array(coordinates) + # Extract shapes and dimensions, then flatten + ndim = coordinates.shape[-1] + output_shape = coordinates.shape[:-1] + flatcoord = np.moveaxis(coordinates, -1, 0).reshape(ndim, -1) + + # Convert coordinates to voxel indices + ijk = np.tensordot( + np.linalg.inv(self.reference.affine), + np.vstack((flatcoord, np.ones((1, flatcoord.shape[1])))), + axes=1) + deltas = ndi.map_coordinates( + self._field, + ijk, + order=order, + mode=mode, + cval=cval, + prefilter=prefilter) + + deltas = np.moveaxis(deltas[0:3, :].reshape((ndim, ) + output_shape), + 0, -1) + + return coordinates + deltas + + +class BSplineFieldTransform(TransformBase): + """Represent a nonlinear transform parameterized by BSpline basis.""" + + __slots__ = ['_coeffs', '_knots', '_refknots', '_order', '_moving'] + __s = (slice(None), ) + + def __init__(self, reference, coefficients, order=3): + """Create a smooth deformation field using B-Spline basis.""" + super(BSplineFieldTransform, self).__init__() + self._order = order + self.reference = reference + + if coefficients.shape[-1] != self.ndim: + raise ValueError( + 'Number of components of the coefficients does ' + 'not match the number of dimensions') + + self._coeffs = coefficients.get_data() + self._knots = ImageSpace(four_to_three(coefficients)[0]) + self._cache_moving() + + def _cache_moving(self): + self._moving = np.zeros((self.reference.shape) + (3, ), + dtype='float32') + ijk = np.moveaxis(self.reference.ndindex, 0, -1).reshape(-1, self.ndim) + xyz = np.moveaxis(self.reference.ndcoords, 0, -1).reshape(-1, self.ndim) + print(np.shape(xyz)) + + for i in range(np.shape(xyz)[0]): + print(i, xyz[i, :]) + self._moving[tuple(ijk[i]) + self.__s] = self._interp_transform(xyz[i, :]) + + def _interp_transform(self, coords): + # Calculate position in the grid of control points + knots_ijk = self._knots.inverse.dot(np.hstack((coords, 1)))[:3] + neighbors = [] + offset = 0.0 if self._order & 1 else 0.5 + # Calculate neighbors along each dimension + for dim in range(self.ndim): + first = int(np.floor(knots_ijk[dim] + offset) - self._order // 2) + neighbors.append(list(range(first, first + self._order + 1))) + + # Get indexes of the neighborings clique + ndindex = np.moveaxis( + np.array(np.meshgrid(*neighbors, indexing='ij')), 0, -1).reshape( + -1, self.ndim) + + # Calculate the tensor B-spline weights of each neighbor + # weights = np.prod(vbspl(ndindex - knots_ijk), axis=-1) + ndindex = [tuple(v) for v in ndindex] + + # Retrieve coefficients and deal with boundary conditions + zero = np.zeros(self.ndim) + shape = np.array(self._knots.shape) + coeffs = [] + for ijk in ndindex: + offbounds = (zero > ijk) | (shape <= ijk) + coeffs.append( + self._coeffs[ijk] if not np.any(offbounds) + else [0.0] * self.ndim) + + coords[:3] += weights.dot(np.array(coeffs, dtype=float)) + return self.reference.inverse.dot(np.hstack((coords, 1)))[:3] + + def map_voxel(self, index, moving=None): + """Apply ijk' = f_ijk((i, j, k)), equivalent to the above with indexes.""" + return tuple(self._moving[index + self.__s]) + + def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True, + output_dtype=None): + """ + Resample the `moving` image applying the deformation field. + + Examples + -------- + >>> import numpy as np + >>> import nibabel as nb + >>> ref = nb.load('t1_weighted.nii.gz') + >>> coeffs = np.zeros((6, 6, 6, 3)) + >>> coeffs[2, 2, 2, ...] = [10.0, -20.0, 0] + >>> aff = ref.affine + >>> aff[:3, :3] = aff.dot(np.eye(3) * np.array(ref.header.get_zooms()[:3] / 6.0) + >>> coeffsimg = nb.Nifti1Image(coeffs, ref.affine, ref.header) + >>> xfm = nb.transform.BSplineFieldTransform(ref, coeffsimg) + >>> new = xfm.resample(ref) + >>> new.to_filename('deffield.nii.gz') + + """ + self._cache_moving() + return super(BSplineFieldTransform, self).resample( + moving, order=order, mode=mode, cval=cval, prefilter=prefilter) + + +def _pprint(inlist): + return 'x'.join(['%d' % s for s in inlist]) diff --git a/nibabel/volumeutils.py b/nibabel/volumeutils.py index 2cc083ecb6..a7aa6f4e7b 100644 --- a/nibabel/volumeutils.py +++ b/nibabel/volumeutils.py @@ -1456,7 +1456,7 @@ def finite_range(arr, check_nan=False): return np.nanmin(mins), np.nanmax(maxes) -def shape_zoom_affine(shape, zooms, x_flip=True): +def shape_zoom_affine(shape, zooms, x_flip=True, y_flip=False): ''' Get affine implied by given shape and zooms We get the translations from the center of the image (implied by @@ -1471,6 +1471,9 @@ def shape_zoom_affine(shape, zooms, x_flip=True): x_flip : {True, False} whether to flip the X row of the affine. Corresponds to radiological storage on disk. + y_flip : {False, True} + whether to flip the Y row of the affine. Corresponds to + DICOM storage on disk when x_flip is also True. Returns ------- @@ -1510,6 +1513,9 @@ def shape_zoom_affine(shape, zooms, x_flip=True): zooms = full_zooms if x_flip: zooms[0] *= -1 + + if y_flip: + zooms[1] *= -1 # Get translations from center of image origin = (shape - 1) / 2.0 aff = np.eye(4)