From 3cfff0d14f5d36950a2661bb262c52ffabc66bac Mon Sep 17 00:00:00 2001 From: oesteban Date: Thu, 19 Sep 2019 21:31:47 -0700 Subject: [PATCH 01/36] ENH: Add an utility to calculate obliquity of affines This implementation mimics the implementation of AFNI. --- nibabel/affines.py | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/nibabel/affines.py b/nibabel/affines.py index a7d7a4e9b8..03ab8b80e8 100644 --- a/nibabel/affines.py +++ b/nibabel/affines.py @@ -3,6 +3,7 @@ """ Utility routines for working with points and affine transforms """ +from math import acos, pi as PI import numpy as np from functools import reduce @@ -296,3 +297,24 @@ 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, in degrees from plumb. + + This implementation is inspired by `AFNI's implementation + `_ + + >>> affine = np.array([ + ... [2.74999725e+00, -2.74999817e-03, 2.74999954e-03, -7.69847980e+01], + ... [2.98603540e-03, 2.73886840e+00, -2.47165887e-01, -8.36692043e+01], + ... [-2.49170222e-03, 2.47168626e-01, 2.73886865e+00, -8.34056848e+01], + ... [ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00]]) + >>> abs(5.15699 - obliquity(affine)) < 0.0001 + True + + """ + vs = voxel_sizes(affine) + fig_merit = (affine[:-1, :-1] / vs[np.newaxis, ...]).max(axis=1).min() + return abs(acos(fig_merit) * 180 / PI) From c92d5609f37ac2ec18f94606b8bb3075759e27ea Mon Sep 17 00:00:00 2001 From: oesteban Date: Fri, 20 Sep 2019 11:26:42 -0700 Subject: [PATCH 02/36] enh(tests): add not-oblique test, move tests to ``test_affines.py`` Addressing one of @matthew-brett's comments. --- nibabel/affines.py | 8 -------- nibabel/tests/test_affines.py | 12 +++++++++++- 2 files changed, 11 insertions(+), 9 deletions(-) diff --git a/nibabel/affines.py b/nibabel/affines.py index 03ab8b80e8..bc800fd318 100644 --- a/nibabel/affines.py +++ b/nibabel/affines.py @@ -306,14 +306,6 @@ def obliquity(affine): This implementation is inspired by `AFNI's implementation `_ - >>> affine = np.array([ - ... [2.74999725e+00, -2.74999817e-03, 2.74999954e-03, -7.69847980e+01], - ... [2.98603540e-03, 2.73886840e+00, -2.47165887e-01, -8.36692043e+01], - ... [-2.49170222e-03, 2.47168626e-01, 2.73886865e+00, -8.34056848e+01], - ... [ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00]]) - >>> abs(5.15699 - obliquity(affine)) < 0.0001 - True - """ vs = voxel_sizes(affine) fig_merit = (affine[:-1, :-1] / vs[np.newaxis, ...]).max(axis=1).min() diff --git a/nibabel/tests/test_affines.py b/nibabel/tests/test_affines.py index e66ed46190..143c80312b 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,13 @@ 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.""" + 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) + assert_almost_equal(obliquity(oblique), 5.1569948883) From da8da564adfccee6e87ec3695b9ce0aa55fa1be9 Mon Sep 17 00:00:00 2001 From: oesteban Date: Sat, 21 Sep 2019 17:46:12 -0700 Subject: [PATCH 03/36] enh: return radians unless degrees=True --- nibabel/affines.py | 7 +++++-- nibabel/tests/test_affines.py | 2 +- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/nibabel/affines.py b/nibabel/affines.py index bc800fd318..462798b6f1 100644 --- a/nibabel/affines.py +++ b/nibabel/affines.py @@ -299,7 +299,7 @@ def voxel_sizes(affine): return np.sqrt(np.sum(top_left ** 2, axis=0)) -def obliquity(affine): +def obliquity(affine, degrees=False): r""" Estimate the obliquity an affine's axes represent, in degrees from plumb. @@ -309,4 +309,7 @@ def obliquity(affine): """ vs = voxel_sizes(affine) fig_merit = (affine[:-1, :-1] / vs[np.newaxis, ...]).max(axis=1).min() - return abs(acos(fig_merit) * 180 / PI) + radians = abs(acos(fig_merit)) + if not degrees: + return radians + return radians * 180 / PI diff --git a/nibabel/tests/test_affines.py b/nibabel/tests/test_affines.py index 143c80312b..db4fbe0630 100644 --- a/nibabel/tests/test_affines.py +++ b/nibabel/tests/test_affines.py @@ -187,4 +187,4 @@ def test_obliquity(): 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) - assert_almost_equal(obliquity(oblique), 5.1569948883) + assert_almost_equal(obliquity(oblique, degrees=True), 5.1569948883) From 253a256d18d75664001ef3d9872c25952675b7f4 Mon Sep 17 00:00:00 2001 From: oesteban Date: Mon, 23 Sep 2019 09:31:42 -0700 Subject: [PATCH 04/36] enh: address @matthew-brett's comments --- nibabel/affines.py | 26 +++++++++++++++++--------- nibabel/tests/test_affines.py | 6 ++++-- 2 files changed, 21 insertions(+), 11 deletions(-) diff --git a/nibabel/affines.py b/nibabel/affines.py index 462798b6f1..25c9a4cdbd 100644 --- a/nibabel/affines.py +++ b/nibabel/affines.py @@ -2,8 +2,6 @@ # vi: set ft=python sts=4 ts=4 sw=4 et: """ Utility routines for working with points and affine transforms """ - -from math import acos, pi as PI import numpy as np from functools import reduce @@ -299,17 +297,27 @@ def voxel_sizes(affine): return np.sqrt(np.sum(top_left ** 2, axis=0)) -def obliquity(affine, degrees=False): +def obliquity(affine): r""" - Estimate the obliquity an affine's axes represent, in degrees from plumb. + 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 `_ + 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) - fig_merit = (affine[:-1, :-1] / vs[np.newaxis, ...]).max(axis=1).min() - radians = abs(acos(fig_merit)) - if not degrees: - return radians - return radians * 180 / PI + best_cosines = np.abs((affine[:-1, :-1] / vs).max(axis=1)) + return np.arccos(best_cosines) diff --git a/nibabel/tests/test_affines.py b/nibabel/tests/test_affines.py index db4fbe0630..13b554c5a8 100644 --- a/nibabel/tests/test_affines.py +++ b/nibabel/tests/test_affines.py @@ -182,9 +182,11 @@ def test_voxel_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) - assert_almost_equal(obliquity(oblique, degrees=True), 5.1569948883) + assert_almost_equal(obliquity(aligned), [0.0, 0.0, 0.0]) + assert_almost_equal(obliquity(oblique) * 180 / pi, + [0.0810285, 5.1569949, 5.1569376]) From 04584b6d8e0e50928ee7c5ff09de61471d6e5acb Mon Sep 17 00:00:00 2001 From: oesteban Date: Tue, 24 Sep 2019 10:20:19 -0700 Subject: [PATCH 05/36] doc: add link to AFNI's documentation about *obliquity* [skip ci] --- nibabel/affines.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/nibabel/affines.py b/nibabel/affines.py index 25c9a4cdbd..b09257c7ad 100644 --- a/nibabel/affines.py +++ b/nibabel/affines.py @@ -304,7 +304,9 @@ def obliquity(affine): 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 ---------- From 70d9620e7ed037e6d2d8ab9b43dc84ccc400b950 Mon Sep 17 00:00:00 2001 From: oesteban Date: Tue, 31 Jul 2018 20:02:46 -0700 Subject: [PATCH 06/36] Initial commit --- nibabel/transform/__init__.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) create mode 100644 nibabel/transform/__init__.py diff --git a/nibabel/transform/__init__.py b/nibabel/transform/__init__.py new file mode 100644 index 0000000000..8e34b518a2 --- /dev/null +++ b/nibabel/transform/__init__.py @@ -0,0 +1,17 @@ +# 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 36a17649041692156d07f3e6fa3760319da880e5 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Wed, 1 Aug 2018 08:02:52 -0700 Subject: [PATCH 07/36] add TransformBase --- nibabel/transform/base.py | 85 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 85 insertions(+) create mode 100644 nibabel/transform/base.py diff --git a/nibabel/transform/base.py b/nibabel/transform/base.py new file mode 100644 index 0000000000..02ec22e97b --- /dev/null +++ b/nibabel/transform/base.py @@ -0,0 +1,85 @@ +# 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 ''' +from __future__ import division, print_function, absolute_import +import numpy as np + + +class ImageSpace(object): + ''' + Abstract class to represent spaces of gridded data (images) + ''' + __slots__ = ['_affine', '_shape'] + + def __init__(self, image): + self._affine = image.affine + self._shape = image.shape + + @property + def affine(self): + return self._affine + + @property + def shape(self): + return self._shape + + +class TransformBase(object): + ''' + Abstract image class to represent transforms + ''' + __slots__ = ['_reference'] + + def __init__(self): + self._reference = None + + @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) + + def resample(self, moving, order=3, dtype=None): + '''A virtual method to resample the moving image in the reference space''' + raise NotImplementedError + + def transform(self, coords, vox=False): + '''Find the coordinates in moving space corresponding to the + input reference coordinates''' + raise NotImplementedError + + +class Affine(TransformBase): + '''Linear transforms''' + __slots__ = ['_affine'] + + def __init__(self, affine): + self._affine = np.array(affine) + assert self._affine.ndim == 2, 'affine matrix should be 2D' + assert self._affine.shape[0] == self._affine.shape[1], 'affine matrix is not square' + super(Affine, self).__init__() + + def resample(self, moving, order=3, dtype=None): + + if dtype is None: + dtype = moving.get_data_dtype() + + # Compose an index to index affine matrix + xfm = self.reference.affine.dot(self._affine.dot(np.linalg.inv(moving.affine))) + + def transform(self, coords): + return self._affine.dot(coords)[:-1] + + + From 2a700c9295ab22e927f04212ac31f377b46e89e7 Mon Sep 17 00:00:00 2001 From: oesteban Date: Wed, 1 Aug 2018 10:21:55 -0700 Subject: [PATCH 08/36] add new module --- nibabel/__init__.py | 2 ++ nibabel/transform/__init__.py | 2 ++ 2 files changed, 4 insertions(+) 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/transform/__init__.py b/nibabel/transform/__init__.py index 8e34b518a2..824b84f357 100644 --- a/nibabel/transform/__init__.py +++ b/nibabel/transform/__init__.py @@ -15,3 +15,5 @@ transform """ +from __future__ import absolute_import +from .base import ImageSpace, Affine From fda2bfe3cc524274fec44ed65aa463307d0c96c4 Mon Sep 17 00:00:00 2001 From: oesteban Date: Wed, 1 Aug 2018 10:23:07 -0700 Subject: [PATCH 09/36] add documentation, core functionality --- nibabel/transform/base.py | 150 +++++++++++++++++++++++++++++++++----- 1 file changed, 131 insertions(+), 19 deletions(-) diff --git a/nibabel/transform/base.py b/nibabel/transform/base.py index 02ec22e97b..63b0ccdb13 100644 --- a/nibabel/transform/base.py +++ b/nibabel/transform/base.py @@ -9,6 +9,7 @@ ''' Common interface for transforms ''' from __future__ import division, print_function, absolute_import import numpy as np +from scipy import ndimage as ndi class ImageSpace(object): @@ -54,32 +55,143 @@ def resample(self, moving, order=3, dtype=None): '''A virtual method to resample the moving image in the reference space''' raise NotImplementedError - def transform(self, coords, vox=False): + def map_point(self, coords): '''Find the coordinates in moving space corresponding to the input reference coordinates''' raise NotImplementedError + def map_voxel(self, index, moving=None): + '''Find the voxel indices in the moving image corresponding to the + input reference voxel''' + raise NotImplementedError -class Affine(TransformBase): - '''Linear transforms''' - __slots__ = ['_affine'] - def __init__(self, affine): - self._affine = np.array(affine) - assert self._affine.ndim == 2, 'affine matrix should be 2D' - assert self._affine.shape[0] == self._affine.shape[1], 'affine matrix is not square' +class Affine(TransformBase): + '''Represents linear transforms on image data''' + __slots__ = ['_matrix'] + + def __init__(self, matrix): + '''Initialize a 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]]) + + ''' + self._matrix = np.array(matrix) + assert self._matrix.ndim == 2, 'affine matrix should be 2D' + assert self._matrix.shape[0] == self._matrix.shape[1], 'affine matrix is not square' super(Affine, self).__init__() - def resample(self, moving, order=3, dtype=None): - - if dtype is None: - dtype = moving.get_data_dtype() + @property + def matrix(self): + 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') + reference = moving # Compose an index to index affine matrix - xfm = self.reference.affine.dot(self._affine.dot(np.linalg.inv(moving.affine))) - - def transform(self, coords): - return self._affine.dot(coords)[:-1] - - - + matrix = reference.affine.dot(self._matrix.dot(np.linalg.inv(moving.affine))) + mdata = moving.get_data() + moved = ndi.affine_transform( + mdata, matrix=matrix[:mdata.ndim, :mdata.ndim], + offset=matrix[:mdata.ndim, mdata.ndim], + output_shape=reference.shape, order=order, mode=mode, + cval=cval, prefilter=prefilter) + + img = moving.__class__(moved, moving.affine, moving.header) + img.header.set_data_dtype(output_dtype) + + return img + + def map_point(self, coords, forward=True): + coords = np.array(coords) + if coords.shape[0] == self._matrix.shape[0] - 1: + coords = np.append(coords, [1]) + affine = self._matrix if forward else np.linalg.inv(self._matrix) + return affine.dot(coords)[:-1] + + def map_voxel(self, index, moving=None): + try: + reference = self.reference + except ValueError: + print('Warning: no reference space defined, using moving as reference') + 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.shape[0] - 1: + index = np.append(index, [1]) + + matrix = reference.affine.dot(self._matrix.dot(np.linalg.inv(moving.affine))) + return matrix.dot(index)[:-1] From 3ef5adc04fb4a8843343b8bd5d7157044410690a Mon Sep 17 00:00:00 2001 From: oesteban Date: Wed, 1 Aug 2018 16:46:54 -0700 Subject: [PATCH 10/36] add a deformation field transform --- nibabel/transform/__init__.py | 1 + nibabel/transform/base.py | 63 ++++++++++++++++++++++++++++++---- nibabel/transform/nonlinear.py | 55 +++++++++++++++++++++++++++++ 3 files changed, 113 insertions(+), 6 deletions(-) create mode 100644 nibabel/transform/nonlinear.py diff --git a/nibabel/transform/__init__.py b/nibabel/transform/__init__.py index 824b84f357..96e8c98204 100644 --- a/nibabel/transform/__init__.py +++ b/nibabel/transform/__init__.py @@ -17,3 +17,4 @@ """ from __future__ import absolute_import from .base import ImageSpace, Affine +from .nonlinear import DeformationFieldTransform diff --git a/nibabel/transform/base.py b/nibabel/transform/base.py index 63b0ccdb13..6e9cada95b 100644 --- a/nibabel/transform/base.py +++ b/nibabel/transform/base.py @@ -51,9 +51,60 @@ def reference(self): def reference(self, image): self._reference = ImageSpace(image) - def resample(self, moving, order=3, dtype=None): - '''A virtual method to resample the moving image in the reference space''' - raise NotImplementedError + 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. + + ''' + + if output_dtype is None: + output_dtype = moving.header.get_data_dtype() + + moving_data = moving.get_data() + 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): '''Find the coordinates in moving space corresponding to the @@ -164,10 +215,10 @@ def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True, output_shape=reference.shape, order=order, mode=mode, cval=cval, prefilter=prefilter) - img = moving.__class__(moved, moving.affine, moving.header) - img.header.set_data_dtype(output_dtype) + moved_image = moving.__class__(moved, reference.affine, moving.header) + moved_image.header.set_data_dtype(output_dtype) - return img + return moved_image def map_point(self, coords, forward=True): coords = np.array(coords) diff --git a/nibabel/transform/nonlinear.py b/nibabel/transform/nonlinear.py new file mode 100644 index 0000000000..ff9f6322ce --- /dev/null +++ b/nibabel/transform/nonlinear.py @@ -0,0 +1,55 @@ +# 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 ''' +from __future__ import division, print_function, absolute_import +import numpy as np + +from .base import TransformBase +from ..funcs import four_to_three + + +class DeformationFieldTransform(TransformBase): + '''Represents linear transforms on image data''' + __slots__ = ['_field'] + + def __init__(self, field): + self.reference = four_to_three(field)[0] + # Set the field components in front + self._field = field.get_data() + + def map_voxel(self, index, moving=None): + ''' + + 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') + + ''' + vox2ras = self.reference.affine + ras2vox = np.linalg.inv(vox2ras if moving is None + else moving.affine) + + index = np.array(index) + if index.shape[0] == vox2ras.shape[0] - 1: + index = np.append(index, [1.0]) + + idx = tuple([int(i) for i in index[:-1]]) + point = vox2ras.dot(index)[:-1] + delta = np.squeeze(self._field[idx + (slice(None), )]) + newindex = ras2vox.dot(np.append(point + delta, [1]))[:-1] + return tuple(newindex) From 189a662ac696c8bd2f84bdf3868232b9a76df26e Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Wed, 1 Aug 2018 20:52:51 -0700 Subject: [PATCH 11/36] optimize deformation field --- nibabel/transform/nonlinear.py | 50 ++++++++++++++++++++++++---------- 1 file changed, 36 insertions(+), 14 deletions(-) diff --git a/nibabel/transform/nonlinear.py b/nibabel/transform/nonlinear.py index ff9f6322ce..33874cdf29 100644 --- a/nibabel/transform/nonlinear.py +++ b/nibabel/transform/nonlinear.py @@ -16,12 +16,40 @@ class DeformationFieldTransform(TransformBase): '''Represents linear transforms on image data''' - __slots__ = ['_field'] + __slots__ = ['_field', '_moving'] def __init__(self, field): + super(DeformationFieldTransform, self).__init__() self.reference = four_to_three(field)[0] - # Set the field components in front + # Cache the new indexes self._field = field.get_data() + 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) + flatxyz = np.tensordot( + self.reference.affine, + np.vstack((flatgrid, np.ones((1, ndim)))), + axes=1 + ) + + delta = np.moveaxis( + flatxyz[0:3, :].reshape((ndim, ) + self._field.shape[:-1]), + 0, -1) + self._moving = self._field + delta def map_voxel(self, index, moving=None): ''' @@ -40,16 +68,10 @@ def map_voxel(self, index, moving=None): >>> new.to_filename('deffield.nii.gz') ''' - vox2ras = self.reference.affine - ras2vox = np.linalg.inv(vox2ras if moving is None - else moving.affine) - - index = np.array(index) - if index.shape[0] == vox2ras.shape[0] - 1: - index = np.append(index, [1.0]) - - idx = tuple([int(i) for i in index[:-1]]) - point = vox2ras.dot(index)[:-1] - delta = np.squeeze(self._field[idx + (slice(None), )]) - newindex = ras2vox.dot(np.append(point + delta, [1]))[:-1] + ras2vox = np.linalg.inv( + self.reference.affine if moving is None + else moving.affine) + + point = self._moving[index + (slice(None), )] + newindex = ras2vox.dot(np.append(point, [1]))[:-1] return tuple(newindex) From 2619ae40484765fe83756f7f6c78923b5a6078a1 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Wed, 1 Aug 2018 22:40:45 -0700 Subject: [PATCH 12/36] pre-cache transformed indexes --- nibabel/transform/nonlinear.py | 23 ++++++++++++++++++----- 1 file changed, 18 insertions(+), 5 deletions(-) diff --git a/nibabel/transform/nonlinear.py b/nibabel/transform/nonlinear.py index 33874cdf29..bf861a1275 100644 --- a/nibabel/transform/nonlinear.py +++ b/nibabel/transform/nonlinear.py @@ -21,8 +21,11 @@ class DeformationFieldTransform(TransformBase): def __init__(self, field): super(DeformationFieldTransform, self).__init__() self.reference = four_to_three(field)[0] - # Cache the new indexes self._field = field.get_data() + self._moving = None + + def _cache_moving(self, moving): + # Cache the new indexes ndim = self._field.ndim - 1 if ndim == 2: grid = np.meshgrid( @@ -42,14 +45,24 @@ def __init__(self, field): flatgrid = grid.reshape(ndim, -1) flatxyz = np.tensordot( self.reference.affine, - np.vstack((flatgrid, np.ones((1, ndim)))), + np.vstack((flatgrid, np.ones((1, flatgrid.shape[1])))), axes=1 ) - delta = np.moveaxis( - flatxyz[0:3, :].reshape((ndim, ) + self._field.shape[:-1]), + newxyz = flatxyz + np.vstack((self._field.reshape(ndim , -1), + np.zeros((1, flatgrid.shape[1])))) + newijk = np.tensordot(np.linalg.inv(moving.affine), + newxyz, axes=1) + + self._moving = np.moveaxis( + newijk[0:3, :].reshape((ndim, ) + self._field.shape[:-1]), 0, -1) - self._moving = self._field + delta + + def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True, + output_dtype=None): + 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): ''' From 4e7c19b21bb545dbf644b5460cdaa1d40ef74813 Mon Sep 17 00:00:00 2001 From: oesteban Date: Wed, 1 Aug 2018 22:52:26 -0700 Subject: [PATCH 13/36] used cached field --- nibabel/transform/nonlinear.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/nibabel/transform/nonlinear.py b/nibabel/transform/nonlinear.py index bf861a1275..be34d8189f 100644 --- a/nibabel/transform/nonlinear.py +++ b/nibabel/transform/nonlinear.py @@ -81,10 +81,5 @@ def map_voxel(self, index, moving=None): >>> new.to_filename('deffield.nii.gz') ''' - ras2vox = np.linalg.inv( - self.reference.affine if moving is None - else moving.affine) + return tuple(self._moving[index + (slice(None), )]) - point = self._moving[index + (slice(None), )] - newindex = ras2vox.dot(np.append(point, [1]))[:-1] - return tuple(newindex) From 23d70025545120bead57a4e3fc058e470dc29541 Mon Sep 17 00:00:00 2001 From: oesteban Date: Thu, 2 Aug 2018 17:16:06 -0700 Subject: [PATCH 14/36] add caching of deltas in voxel coordinates --- nibabel/transform/base.py | 9 ++++++- nibabel/transform/nonlinear.py | 46 +++++++++++++++++++++++++--------- 2 files changed, 42 insertions(+), 13 deletions(-) diff --git a/nibabel/transform/base.py b/nibabel/transform/base.py index 6e9cada95b..52b1e25141 100644 --- a/nibabel/transform/base.py +++ b/nibabel/transform/base.py @@ -30,6 +30,13 @@ def affine(self): def shape(self): return self._shape + def __eq__(self, other): + try: + return np.allclose(self.affine, other.affine) and self.shape == other.shape + except AttributeError: + pass + return False + class TransformBase(object): ''' @@ -245,4 +252,4 @@ def map_voxel(self, index, moving=None): index = np.append(index, [1]) matrix = reference.affine.dot(self._matrix.dot(np.linalg.inv(moving.affine))) - return matrix.dot(index)[:-1] + return tuple(matrix.dot(index)[:-1]) diff --git a/nibabel/transform/nonlinear.py b/nibabel/transform/nonlinear.py index be34d8189f..70eefeeea7 100644 --- a/nibabel/transform/nonlinear.py +++ b/nibabel/transform/nonlinear.py @@ -9,23 +9,37 @@ ''' Common interface for transforms ''' from __future__ import division, print_function, absolute_import import numpy as np +from scipy import ndimage as ndi -from .base import TransformBase +from .base import ImageSpace, TransformBase from ..funcs import four_to_three class DeformationFieldTransform(TransformBase): '''Represents linear transforms on image data''' - __slots__ = ['_field', '_moving'] + __slots__ = ['_field', '_moving', '_moving_space'] + __s = (slice(None), ) def __init__(self, field): + ''' + Create a dense deformation field transform + ''' super(DeformationFieldTransform, self).__init__() + # By definition, a free deformation field has a + # displacement vector per voxel in output (reference) + # space self.reference = four_to_three(field)[0] self._field = field.get_data() - self._moving = None + self._moving = None # Where each voxel maps to + self._moving_space = None # Input space cache def _cache_moving(self, moving): - # Cache the new indexes + # 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( @@ -43,28 +57,32 @@ def _cache_moving(self, moving): 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 ) - newxyz = flatxyz + np.vstack((self._field.reshape(ndim , -1), - np.zeros((1, flatgrid.shape[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): - 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): ''' Examples @@ -81,5 +99,9 @@ def map_voxel(self, index, moving=None): >>> new.to_filename('deffield.nii.gz') ''' - return tuple(self._moving[index + (slice(None), )]) + 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): + return tuple(self._moving[index + self.__s]) From 2acd9f39c4500181ef5c52c238bcd880d0cd1f29 Mon Sep 17 00:00:00 2001 From: oesteban Date: Fri, 3 Aug 2018 11:58:45 -0700 Subject: [PATCH 15/36] add bspline cython extension --- .gitignore | 1 + nibabel/maths/__init__.py | 17 ++++++++++++++ nibabel/maths/bspline.pyx | 43 ++++++++++++++++++++++++++++++++++ nibabel/transform/nonlinear.py | 28 ++++++++++++++++++++++ 4 files changed, 89 insertions(+) create mode 100644 nibabel/maths/__init__.py create mode 100644 nibabel/maths/bspline.pyx diff --git a/.gitignore b/.gitignore index df018f0ead..bbd0a11b43 100644 --- a/.gitignore +++ b/.gitignore @@ -84,3 +84,4 @@ Thumbs.db doc/source/reference venv/ .buildbot.patch +nibabel/maths/bspline.c diff --git a/nibabel/maths/__init__.py b/nibabel/maths/__init__.py new file mode 100644 index 0000000000..33bd530321 --- /dev/null +++ b/nibabel/maths/__init__.py @@ -0,0 +1,17 @@ +# 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.maths + +.. autosummary:: + :toctree: ../generated + + maths +""" diff --git a/nibabel/maths/bspline.pyx b/nibabel/maths/bspline.pyx new file mode 100644 index 0000000000..e357ca6dfb --- /dev/null +++ b/nibabel/maths/bspline.pyx @@ -0,0 +1,43 @@ +# 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. +# +### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## +''' Cython math extension ''' +from __future__ import division +import numpy as np +cimport numpy as np +cimport cython +from libc.math cimport fabs + + +cdef double c_cubic(double x) nogil: + cdef: + double x_t = fabs(x) + + if x_t >= 2.0: + return(0.0) + if x_t <= 1.0: + return(2.0 / 3.0 - x_t**2 + 0.5 * x_t**3) + elif x_t <= 2.0: + return((2 - x_t)**3 / 6.0) + + +def cubic(double x): + """ + Evaluate the univariate cubic bspline at x + + Pure python implementation: :: + + def bspl(x): + if x >= 2.0: + return 0.0 + if x <= 1.0: + return 2.0 / 3.0 - x**2 + 0.5 * x**3 + elif x <= 2.0: + return (2 - x)**3 / 6.0 + """ + return(c_cubic(x)) diff --git a/nibabel/transform/nonlinear.py b/nibabel/transform/nonlinear.py index 70eefeeea7..ecf7ad0a3a 100644 --- a/nibabel/transform/nonlinear.py +++ b/nibabel/transform/nonlinear.py @@ -105,3 +105,31 @@ def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True, def map_voxel(self, index, moving=None): return tuple(self._moving[index + self.__s]) + + def map_coordinates(self, coordinates, order=3, mode='constant', cval=0.0, + prefilter=True): + 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) + + print(deltas) + + deltas = np.moveaxis(deltas[0:3, :].reshape((ndim, ) + output_shape), + 0, -1) + + return coordinates + deltas From 385ceff957d1e6b1430d1ff6769f358461223e8d Mon Sep 17 00:00:00 2001 From: oesteban Date: Fri, 3 Aug 2018 18:02:55 -0700 Subject: [PATCH 16/36] a smarter ImageSpace --- nibabel/transform/base.py | 47 +++++++++++++++++++++++++++++++++++---- 1 file changed, 43 insertions(+), 4 deletions(-) diff --git a/nibabel/transform/base.py b/nibabel/transform/base.py index 52b1e25141..2027741d85 100644 --- a/nibabel/transform/base.py +++ b/nibabel/transform/base.py @@ -13,14 +13,18 @@ class ImageSpace(object): - ''' - Abstract class to represent spaces of gridded data (images) - ''' - __slots__ = ['_affine', '_shape'] + '''Class to represent spaces of gridded data (images)''' + __slots__ = ['_affine', '_shape', '_ndim', '_ndindex', '_coords', '_nvox'] def __init__(self, image): self._affine = image.affine self._shape = image.shape + self._ndim = len(image.shape) + self._nvox = np.prod(image.shape) + self._ndindex = None + self._coords = None + if self._ndim not in [2, 3]: + raise ValueError('Invalid image space (%d-D)' % self._ndim) @property def affine(self): @@ -30,6 +34,41 @@ def affine(self): 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 __eq__(self, other): try: return np.allclose(self.affine, other.affine) and self.shape == other.shape From fa030a42f982eab348be125a62cddfe478f82e20 Mon Sep 17 00:00:00 2001 From: oesteban Date: Fri, 10 Aug 2018 16:36:39 -0700 Subject: [PATCH 17/36] add comment --- nibabel/transform/base.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nibabel/transform/base.py b/nibabel/transform/base.py index 2027741d85..bf4940bc75 100644 --- a/nibabel/transform/base.py +++ b/nibabel/transform/base.py @@ -20,7 +20,7 @@ def __init__(self, image): self._affine = image.affine self._shape = image.shape self._ndim = len(image.shape) - self._nvox = np.prod(image.shape) + self._nvox = np.prod(image.shape) # Do not access data array self._ndindex = None self._coords = None if self._ndim not in [2, 3]: From fd418fe100c6c13eee6e6e98bf44bb3be25b2cd2 Mon Sep 17 00:00:00 2001 From: oesteban Date: Fri, 10 Aug 2018 16:43:07 -0700 Subject: [PATCH 18/36] remove cython module --- nibabel/maths/__init__.py | 17 ---------------- nibabel/maths/bspline.pyx | 43 --------------------------------------- 2 files changed, 60 deletions(-) delete mode 100644 nibabel/maths/__init__.py delete mode 100644 nibabel/maths/bspline.pyx diff --git a/nibabel/maths/__init__.py b/nibabel/maths/__init__.py deleted file mode 100644 index 33bd530321..0000000000 --- a/nibabel/maths/__init__.py +++ /dev/null @@ -1,17 +0,0 @@ -# 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.maths - -.. autosummary:: - :toctree: ../generated - - maths -""" diff --git a/nibabel/maths/bspline.pyx b/nibabel/maths/bspline.pyx deleted file mode 100644 index e357ca6dfb..0000000000 --- a/nibabel/maths/bspline.pyx +++ /dev/null @@ -1,43 +0,0 @@ -# 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. -# -### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## -''' Cython math extension ''' -from __future__ import division -import numpy as np -cimport numpy as np -cimport cython -from libc.math cimport fabs - - -cdef double c_cubic(double x) nogil: - cdef: - double x_t = fabs(x) - - if x_t >= 2.0: - return(0.0) - if x_t <= 1.0: - return(2.0 / 3.0 - x_t**2 + 0.5 * x_t**3) - elif x_t <= 2.0: - return((2 - x_t)**3 / 6.0) - - -def cubic(double x): - """ - Evaluate the univariate cubic bspline at x - - Pure python implementation: :: - - def bspl(x): - if x >= 2.0: - return 0.0 - if x <= 1.0: - return 2.0 / 3.0 - x**2 + 0.5 * x**3 - elif x <= 2.0: - return (2 - x)**3 / 6.0 - """ - return(c_cubic(x)) From 13d722fc431b140125e4121798f9c349686a4b21 Mon Sep 17 00:00:00 2001 From: oesteban Date: Fri, 10 Aug 2018 17:23:44 -0700 Subject: [PATCH 19/36] cleanup --- nibabel/transform/nonlinear.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/nibabel/transform/nonlinear.py b/nibabel/transform/nonlinear.py index ecf7ad0a3a..7dca0e47ca 100644 --- a/nibabel/transform/nonlinear.py +++ b/nibabel/transform/nonlinear.py @@ -10,6 +10,7 @@ from __future__ import division, print_function, absolute_import import numpy as np from scipy import ndimage as ndi +# from gridbspline.interpolate import BsplineNDInterpolator from .base import ImageSpace, TransformBase from ..funcs import four_to_three @@ -106,7 +107,7 @@ def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True, def map_voxel(self, index, moving=None): return tuple(self._moving[index + self.__s]) - def map_coordinates(self, coordinates, order=3, mode='constant', cval=0.0, + def map_coordinates(self, coordinates, order=3, mode='mirror', cval=0.0, prefilter=True): coordinates = np.array(coordinates) # Extract shapes and dimensions, then flatten @@ -127,8 +128,6 @@ def map_coordinates(self, coordinates, order=3, mode='constant', cval=0.0, cval=cval, prefilter=prefilter) - print(deltas) - deltas = np.moveaxis(deltas[0:3, :].reshape((ndim, ) + output_shape), 0, -1) From 253aed5a70f5863e995836fc49d52bb89891c27a Mon Sep 17 00:00:00 2001 From: oesteban Date: Fri, 10 Aug 2018 18:45:34 -0700 Subject: [PATCH 20/36] starting with bspline transform --- nibabel/transform/base.py | 4 ++ nibabel/transform/nonlinear.py | 117 +++++++++++++++++++++++++++++++-- 2 files changed, 115 insertions(+), 6 deletions(-) diff --git a/nibabel/transform/base.py b/nibabel/transform/base.py index bf4940bc75..fa699a2686 100644 --- a/nibabel/transform/base.py +++ b/nibabel/transform/base.py @@ -97,6 +97,10 @@ def reference(self): 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 diff --git a/nibabel/transform/nonlinear.py b/nibabel/transform/nonlinear.py index 7dca0e47ca..f954381735 100644 --- a/nibabel/transform/nonlinear.py +++ b/nibabel/transform/nonlinear.py @@ -17,23 +17,40 @@ class DeformationFieldTransform(TransformBase): - '''Represents linear transforms on image data''' + '''Represents a dense field of displacements (one vector per voxel)''' __slots__ = ['_field', '_moving', '_moving_space'] __s = (slice(None), ) - def __init__(self, field): + def __init__(self, field, reference=None): ''' Create a dense deformation field transform ''' super(DeformationFieldTransform, self).__init__() - # By definition, a free deformation field has a - # displacement vector per voxel in output (reference) - # space - self.reference = four_to_three(field)[0] 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) @@ -132,3 +149,91 @@ def map_coordinates(self, coordinates, order=3, mode='mirror', cval=0.0, 0, -1) return coordinates + deltas + + +class BSplineFieldTransform(TransformBase): + __slots__ = ['_coeffs', '_knots', '_refknots'] + + def __init__(self, reference, coefficients, order=3): + '''Create a smooth deformation field using B-Spline basis''' + super(BSplineFieldTransform, self).__init__() + 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().reshape(-1, self.ndim) + self._knots = ImageSpace(four_to_three(coefficients)[0]) + + # Calculate the voxel coordinates of the reference image + # in the B-spline basis space. + self._refknots = np.tensordot( + np.linalg.inv(self.knots.affine), + np.vstack((self.reference.ndcoords, np.ones((1, self.reference.nvox)))), + axes=1)[..., :3].reshape(self.reference.shape + (self._knots.ndim, )) + + + def map_voxel(self, index, moving=None): + + indexes = [] + offset = 0.0 if self._order & 1 else 0.5 + for dim in range(self.ndim): + first = int(np.floor(xi[dim] + offset) - self._order // 2) + indexes.append(list(range(first, first + self._order + 1))) + + ndindex = np.moveaxis( + np.array(np.meshgrid(*indexes, indexing='ij')), 0, -1).reshape( + -1, self.ndim) + + vbspl = np.vectorize(cubic) + weights = np.prod(vbspl(ndindex - xi), axis=-1) + ndindex = [tuple(v) for v in ndindex] + + zero = np.zeros(self.ndim) + shape = np.array(self.shape) + coeffs = [] + for ijk in ndindex: + offbounds = (zero > ijk) | (shape <= ijk) + if np.any(offbounds): + # Deal with offbounds samples + if self._off_bounds == 'constant': + coeffs.append([self._fill_value] * self.ncomp) + continue + ijk = np.array(ijk, dtype=int) + ijk[ijk < 0] *= -1 + ijk[ijk >= shape] = (2 * shape[ijk >= shape] - ijk[ijk >= shape] - 1).astype(int) + ijk = tuple(ijk.tolist()) + + coeffs.append(self._coeffs[ijk]) + return weights.dot(np.array(coeffs, dtype=float)) + + + # def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True, + # output_dtype=None): + # ''' + + # 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') + + # ''' + + # 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]) From 47c2a5739b63012f94c8012d79efdb49861246d7 Mon Sep 17 00:00:00 2001 From: oesteban Date: Mon, 13 Aug 2018 08:43:26 -0700 Subject: [PATCH 21/36] finishing b-spline interpolation --- nibabel/transform/nonlinear.py | 33 ++++++++++++++------------------- 1 file changed, 14 insertions(+), 19 deletions(-) diff --git a/nibabel/transform/nonlinear.py b/nibabel/transform/nonlinear.py index f954381735..2e09d208b8 100644 --- a/nibabel/transform/nonlinear.py +++ b/nibabel/transform/nonlinear.py @@ -10,7 +10,7 @@ from __future__ import division, print_function, absolute_import import numpy as np from scipy import ndimage as ndi -# from gridbspline.interpolate import BsplineNDInterpolator +from gridbspline.maths import cubic from .base import ImageSpace, TransformBase from ..funcs import four_to_three @@ -174,39 +174,34 @@ def __init__(self, reference, coefficients, order=3): np.vstack((self.reference.ndcoords, np.ones((1, self.reference.nvox)))), axes=1)[..., :3].reshape(self.reference.shape + (self._knots.ndim, )) - def map_voxel(self, index, moving=None): - - indexes = [] + '''Find the corresponding coordinates for a voxel in reference space''' + 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(xi[dim] + offset) - self._order // 2) - indexes.append(list(range(first, first + self._order + 1))) + first = int(np.floor(index[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(*indexes, indexing='ij')), 0, -1).reshape( + np.array(np.meshgrid(*neighbors, indexing='ij')), 0, -1).reshape( -1, self.ndim) + # Calculate the tensor B-spline weights of each neighbor vbspl = np.vectorize(cubic) - weights = np.prod(vbspl(ndindex - xi), axis=-1) + weights = np.prod(vbspl(ndindex - index), 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.shape) coeffs = [] for ijk in ndindex: offbounds = (zero > ijk) | (shape <= ijk) - if np.any(offbounds): - # Deal with offbounds samples - if self._off_bounds == 'constant': - coeffs.append([self._fill_value] * self.ncomp) - continue - ijk = np.array(ijk, dtype=int) - ijk[ijk < 0] *= -1 - ijk[ijk >= shape] = (2 * shape[ijk >= shape] - ijk[ijk >= shape] - 1).astype(int) - ijk = tuple(ijk.tolist()) - - coeffs.append(self._coeffs[ijk]) + coeffs.append( + self._coeffs[ijk] if not np.any(offbounds) + else [0.0] * self.ndim) return weights.dot(np.array(coeffs, dtype=float)) From 4478c062a336fd58b35e7d517bcfc8eef8d65291 Mon Sep 17 00:00:00 2001 From: oesteban Date: Mon, 11 Mar 2019 17:10:22 -0700 Subject: [PATCH 22/36] Add base implementation of transforms and change base implementation of linear and nonlinear transforms --- nibabel/transform/base.py | 140 +++--------------------------- nibabel/transform/linear.py | 151 +++++++++++++++++++++++++++++++++ nibabel/transform/nonlinear.py | 86 +++++++++++-------- 3 files changed, 211 insertions(+), 166 deletions(-) create mode 100644 nibabel/transform/linear.py diff --git a/nibabel/transform/base.py b/nibabel/transform/base.py index fa699a2686..389872ce1c 100644 --- a/nibabel/transform/base.py +++ b/nibabel/transform/base.py @@ -14,7 +14,8 @@ class ImageSpace(object): '''Class to represent spaces of gridded data (images)''' - __slots__ = ['_affine', '_shape', '_ndim', '_ndindex', '_coords', '_nvox'] + __slots__ = ['_affine', '_shape', '_ndim', '_ndindex', '_coords', '_nvox', + '_inverse'] def __init__(self, image): self._affine = image.affine @@ -23,6 +24,7 @@ def __init__(self, image): 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) @@ -30,6 +32,10 @@ def __init__(self, image): def affine(self): return self._affine + @property + def inverse(self): + return self._inverse + @property def shape(self): return self._shape @@ -166,133 +172,7 @@ def map_voxel(self, index, moving=None): input reference voxel''' raise NotImplementedError - -class Affine(TransformBase): - '''Represents linear transforms on image data''' - __slots__ = ['_matrix'] - - def __init__(self, matrix): - '''Initialize a 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]]) - + def to_filename(self, filename): + '''Store the transform in BIDS-Transforms HDF5 file format (.x5). ''' - self._matrix = np.array(matrix) - assert self._matrix.ndim == 2, 'affine matrix should be 2D' - assert self._matrix.shape[0] == self._matrix.shape[1], 'affine matrix is not square' - super(Affine, self).__init__() - - @property - def matrix(self): - 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') - reference = moving - - # Compose an index to index affine matrix - matrix = reference.affine.dot(self._matrix.dot(np.linalg.inv(moving.affine))) - mdata = moving.get_data() - moved = ndi.affine_transform( - mdata, matrix=matrix[:mdata.ndim, :mdata.ndim], - offset=matrix[:mdata.ndim, mdata.ndim], - output_shape=reference.shape, order=order, mode=mode, - cval=cval, prefilter=prefilter) - - moved_image = moving.__class__(moved, reference.affine, moving.header) - moved_image.header.set_data_dtype(output_dtype) - - return moved_image - - def map_point(self, coords, forward=True): - coords = np.array(coords) - if coords.shape[0] == self._matrix.shape[0] - 1: - coords = np.append(coords, [1]) - affine = self._matrix if forward else np.linalg.inv(self._matrix) - return affine.dot(coords)[:-1] - - def map_voxel(self, index, moving=None): - try: - reference = self.reference - except ValueError: - print('Warning: no reference space defined, using moving as reference') - 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.shape[0] - 1: - index = np.append(index, [1]) - - matrix = reference.affine.dot(self._matrix.dot(np.linalg.inv(moving.affine))) - return tuple(matrix.dot(index)[:-1]) + raise NotImplementedError diff --git a/nibabel/transform/linear.py b/nibabel/transform/linear.py new file mode 100644 index 0000000000..83bfc9d8dd --- /dev/null +++ b/nibabel/transform/linear.py @@ -0,0 +1,151 @@ +# 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 ''' +from __future__ import division, print_function, absolute_import +import numpy as np +from scipy import ndimage as ndi + +from .base import ImageSpace, TransformBase +from ..funcs import four_to_three + + +class Affine(TransformBase): + '''Represents linear transforms on image data''' + __slots__ = ['_matrix'] + + def __init__(self, matrix): + '''Initialize a 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]]) + + ''' + self._matrix = np.array(matrix) + assert self._matrix.ndim == 2, 'affine matrix should be 2D' + assert self._matrix.shape[0] == self._matrix.shape[1], 'affine matrix is not square' + super(Affine, self).__init__() + + @property + def matrix(self): + 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') + reference = moving + + # Compose an index to index affine matrix + matrix = reference.affine.dot(self._matrix.dot(np.linalg.inv(moving.affine))) + mdata = moving.get_data() + moved = ndi.affine_transform( + mdata, matrix=matrix[:mdata.ndim, :mdata.ndim], + offset=matrix[:mdata.ndim, mdata.ndim], + output_shape=reference.shape, order=order, mode=mode, + cval=cval, prefilter=prefilter) + + moved_image = moving.__class__(moved, reference.affine, moving.header) + moved_image.header.set_data_dtype(output_dtype) + + return moved_image + + def map_point(self, coords, forward=True): + coords = np.array(coords) + if coords.shape[0] == self._matrix.shape[0] - 1: + coords = np.append(coords, [1]) + affine = self._matrix if forward else np.linalg.inv(self._matrix) + return affine.dot(coords)[:-1] + + def map_voxel(self, index, moving=None): + try: + reference = self.reference + except ValueError: + print('Warning: no reference space defined, using moving as reference') + 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.shape[0] - 1: + index = np.append(index, [1]) + + matrix = reference.affine.dot(self._matrix.dot(np.linalg.inv(moving.affine))) + return tuple(matrix.dot(index)[:-1]) + + def to_filename(self, filename): + '''Store the transform in BIDS-Transforms HDF5 file format (.x5). + ''' + raise NotImplementedError diff --git a/nibabel/transform/nonlinear.py b/nibabel/transform/nonlinear.py index 2e09d208b8..3b4a8fb6c8 100644 --- a/nibabel/transform/nonlinear.py +++ b/nibabel/transform/nonlinear.py @@ -15,6 +15,8 @@ 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)''' @@ -152,11 +154,13 @@ def map_coordinates(self, coordinates, order=3, mode='mirror', cval=0.0, class BSplineFieldTransform(TransformBase): - __slots__ = ['_coeffs', '_knots', '_refknots'] + __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: @@ -164,23 +168,29 @@ def __init__(self, reference, coefficients, order=3): 'Number of components of the coefficients does ' 'not match the number of dimensions') - self._coeffs = coefficients.get_data().reshape(-1, self.ndim) + self._coeffs = coefficients.get_data() self._knots = ImageSpace(four_to_three(coefficients)[0]) - - # Calculate the voxel coordinates of the reference image - # in the B-spline basis space. - self._refknots = np.tensordot( - np.linalg.inv(self.knots.affine), - np.vstack((self.reference.ndcoords, np.ones((1, self.reference.nvox)))), - axes=1)[..., :3].reshape(self.reference.shape + (self._knots.ndim, )) - - def map_voxel(self, index, moving=None): - '''Find the corresponding coordinates for a voxel in reference space''' + 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(index[dim] + offset) - self._order // 2) + 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 @@ -189,45 +199,49 @@ def map_voxel(self, index, moving=None): -1, self.ndim) # Calculate the tensor B-spline weights of each neighbor - vbspl = np.vectorize(cubic) - weights = np.prod(vbspl(ndindex - index), axis=-1) + 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.shape) + 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) - return weights.dot(np.array(coeffs, dtype=float)) + coords[:3] += weights.dot(np.array(coeffs, dtype=float)) + return self.reference.inverse.dot(np.hstack((coords, 1)))[:3] - # def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True, - # output_dtype=None): - # ''' + def map_voxel(self, index, moving=None): + '''Find the corresponding coordinates for a voxel in reference space''' + return tuple(self._moving[index + self.__s]) - # Examples - # -------- + def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True, + output_dtype=None): + ''' - # >>> 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') + 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') - # return super(BSplineFieldTransform, self).resample( - # moving, order=order, mode=mode, cval=cval, prefilter=prefilter) + ''' + self._cache_moving() + return super(BSplineFieldTransform, self).resample( + moving, order=order, mode=mode, cval=cval, prefilter=prefilter) def _pprint(inlist): From 99fa15d67335df96eb2a1e4db85687e9a299a008 Mon Sep 17 00:00:00 2001 From: oesteban Date: Wed, 13 Mar 2019 10:21:18 -0700 Subject: [PATCH 23/36] export to hdf5 --- nibabel/transform/__init__.py | 5 ++++- nibabel/transform/linear.py | 19 +++++++++++++++++-- nibabel/transform/nonlinear.py | 6 +++--- 3 files changed, 24 insertions(+), 6 deletions(-) diff --git a/nibabel/transform/__init__.py b/nibabel/transform/__init__.py index 96e8c98204..323a353cc7 100644 --- a/nibabel/transform/__init__.py +++ b/nibabel/transform/__init__.py @@ -16,5 +16,8 @@ transform """ from __future__ import absolute_import -from .base import ImageSpace, Affine +from .linear import Affine from .nonlinear import DeformationFieldTransform + + +__all__ = ['Affine', 'DeformationFieldTransform'] diff --git a/nibabel/transform/linear.py b/nibabel/transform/linear.py index 83bfc9d8dd..d9f55fa486 100644 --- a/nibabel/transform/linear.py +++ b/nibabel/transform/linear.py @@ -8,6 +8,7 @@ ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## ''' Linear transforms ''' from __future__ import division, print_function, absolute_import +import h5py import numpy as np from scipy import ndimage as ndi @@ -43,7 +44,7 @@ def __init__(self, matrix): ''' self._matrix = np.array(matrix) - assert self._matrix.ndim == 2, 'affine matrix should be 2D' + assert self._matrix.ndim in (2, 3), 'affine matrix should be 2D or 3D' assert self._matrix.shape[0] == self._matrix.shape[1], 'affine matrix is not square' super(Affine, self).__init__() @@ -148,4 +149,18 @@ def map_voxel(self, index, moving=None): def to_filename(self, filename): '''Store the transform in BIDS-Transforms HDF5 file format (.x5). ''' - raise NotImplementedError + with h5py.File(filename, 'w') as out_file: + root = out_file.create_group('/0') + root.create_dataset('/0/Transform', data=self._matrix[:3, :]) + root.create_dataset('/0/Inverse', data=np.linalg.inv(self._matrix)[:3, :]) + root['Type'] = 'affine' + + if self._reference: + refgrp = root.create_group('/0/Reference') + refgrp.create_dataset('/0/Reference/vox2ras', data=self.reference.affine) + refgrp.create_dataset('/0/Reference/ras2vox', data=self.reference.inverse) + refgrp.create_dataset('/0/Reference/shape', data=self.reference.shape) + refgrp.create_dataset('/0/Reference/ndim', data=self.reference.ndim) + refgrp['Type'] = 'image' + + return filename diff --git a/nibabel/transform/nonlinear.py b/nibabel/transform/nonlinear.py index 3b4a8fb6c8..a538628d99 100644 --- a/nibabel/transform/nonlinear.py +++ b/nibabel/transform/nonlinear.py @@ -10,12 +10,12 @@ from __future__ import division, print_function, absolute_import import numpy as np from scipy import ndimage as ndi -from gridbspline.maths import cubic +# from gridbspline.maths import cubic from .base import ImageSpace, TransformBase from ..funcs import four_to_three -vbspl = np.vectorize(cubic) +# vbspl = np.vectorize(cubic) class DeformationFieldTransform(TransformBase): @@ -199,7 +199,7 @@ def _interp_transform(self, coords): -1, self.ndim) # Calculate the tensor B-spline weights of each neighbor - weights = np.prod(vbspl(ndindex - knots_ijk), axis=-1) + # weights = np.prod(vbspl(ndindex - knots_ijk), axis=-1) ndindex = [tuple(v) for v in ndindex] # Retrieve coefficients and deal with boundary conditions From 5eed01db05b088685e51e4eb2eb1e2cdb071a38e Mon Sep 17 00:00:00 2001 From: oesteban Date: Thu, 14 Mar 2019 12:30:46 -0700 Subject: [PATCH 24/36] corrections to adhere the current x5 format draft --- nibabel/transform/base.py | 14 ++++++++++++- nibabel/transform/linear.py | 39 ++++++++++++++++--------------------- 2 files changed, 30 insertions(+), 23 deletions(-) diff --git a/nibabel/transform/base.py b/nibabel/transform/base.py index 389872ce1c..9283092da4 100644 --- a/nibabel/transform/base.py +++ b/nibabel/transform/base.py @@ -9,6 +9,8 @@ ''' Common interface for transforms ''' from __future__ import division, print_function, absolute_import import numpy as np +import h5py + from scipy import ndimage as ndi @@ -172,7 +174,17 @@ def map_voxel(self, index, moving=None): input reference voxel''' raise NotImplementedError + def _to_hdf5(self, x5_root): + '''Serialize this object into the x5 file format''' + raise NotImplementedError + def to_filename(self, filename): '''Store the transform in BIDS-Transforms HDF5 file format (.x5). ''' - raise NotImplementedError + with h5py.File(filename, 'w') as out_file: + out_file['Format'] = 'X5' + out_file['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 index d9f55fa486..7829defe54 100644 --- a/nibabel/transform/linear.py +++ b/nibabel/transform/linear.py @@ -8,19 +8,17 @@ ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## ''' Linear transforms ''' from __future__ import division, print_function, absolute_import -import h5py import numpy as np from scipy import ndimage as ndi -from .base import ImageSpace, TransformBase -from ..funcs import four_to_three +from .base import TransformBase class Affine(TransformBase): '''Represents linear transforms on image data''' __slots__ = ['_matrix'] - def __init__(self, matrix): + def __init__(self, matrix=None): '''Initialize a transform Parameters @@ -43,6 +41,9 @@ def __init__(self, matrix): [0, 0, 0, 1]]) ''' + if matrix is None: + matrix = np.eye(4) + self._matrix = np.array(matrix) assert self._matrix.ndim in (2, 3), 'affine matrix should be 2D or 3D' assert self._matrix.shape[0] == self._matrix.shape[1], 'affine matrix is not square' @@ -146,21 +147,15 @@ def map_voxel(self, index, moving=None): matrix = reference.affine.dot(self._matrix.dot(np.linalg.inv(moving.affine))) return tuple(matrix.dot(index)[:-1]) - def to_filename(self, filename): - '''Store the transform in BIDS-Transforms HDF5 file format (.x5). - ''' - with h5py.File(filename, 'w') as out_file: - root = out_file.create_group('/0') - root.create_dataset('/0/Transform', data=self._matrix[:3, :]) - root.create_dataset('/0/Inverse', data=np.linalg.inv(self._matrix)[:3, :]) - root['Type'] = 'affine' - - if self._reference: - refgrp = root.create_group('/0/Reference') - refgrp.create_dataset('/0/Reference/vox2ras', data=self.reference.affine) - refgrp.create_dataset('/0/Reference/ras2vox', data=self.reference.inverse) - refgrp.create_dataset('/0/Reference/shape', data=self.reference.shape) - refgrp.create_dataset('/0/Reference/ndim', data=self.reference.ndim) - refgrp['Type'] = 'image' - - return filename + def _to_hdf5(self, x5_root): + print('to hdf5') + x5_root.create_dataset('Transform', data=self._matrix[:3, :]) + x5_root.create_dataset('Inverse', data=np.linalg.inv(self._matrix)[:3, :]) + x5_root['Type'] = 'affine' + + if self._reference: + refgrp = x5_root.create_group('Reference') + refgrp.create_dataset('affine', data=self.reference.affine) + refgrp.create_dataset('shape', data=self.reference.shape) + refgrp.create_dataset('ndim', data=self.reference.ndim) + refgrp['Type'] = 'image' From 632e068877e3f51fd9fe72eeaa017d785433dcc3 Mon Sep 17 00:00:00 2001 From: oesteban Date: Thu, 14 Mar 2019 17:16:16 -0700 Subject: [PATCH 25/36] fix coordinates translation in affines, import itk affines --- nibabel/transform/linear.py | 21 ++++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/nibabel/transform/linear.py b/nibabel/transform/linear.py index 7829defe54..7dad928784 100644 --- a/nibabel/transform/linear.py +++ b/nibabel/transform/linear.py @@ -107,7 +107,7 @@ def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True, reference = moving # Compose an index to index affine matrix - matrix = reference.affine.dot(self._matrix.dot(np.linalg.inv(moving.affine))) + matrix = np.linalg.inv(moving.affine).dot(self._matrix.dot(reference.affine)) mdata = moving.get_data() moved = ndi.affine_transform( mdata, matrix=matrix[:mdata.ndim, :mdata.ndim], @@ -159,3 +159,22 @@ def _to_hdf5(self, x5_root): refgrp.create_dataset('shape', data=self.reference.shape) refgrp.create_dataset('ndim', data=self.reference.ndim) refgrp['Type'] = 'image' + + +def load(filename): + ''' Load a linear transform ''' + from ..affines import from_matvec + with open(filename) as itkfile: + itkxfm = itkfile.read().splitlines() + + parameters = np.fromstring(itkxfm[3].split(':')[-1].strip(), dtype=float, sep=' ') + offset = np.fromstring(itkxfm[4].split(':')[-1].strip(), dtype=float, sep=' ') + if len(parameters) == 12: + matrix = np.vstack((parameters.reshape((4, 3)).T, [0, 0, 0, 1])) + matrix[:2, -1] *= -1.0 # Traslation - LPS to RAS + c_neg = from_matvec(np.eye(3), offset * -1.0) + c_pos = from_matvec(np.eye(3), offset) + matrix = np.diag([-1, -1, 1, 1]).dot( + c_pos.dot(matrix.dot(c_neg.dot(np.diag([-1, -1, 1, 1]))))) + + return Affine(matrix) From f2c062e49da8da6d157e652dcc7501e1315ae58d Mon Sep 17 00:00:00 2001 From: oesteban Date: Fri, 15 Mar 2019 18:04:38 -0700 Subject: [PATCH 26/36] wip: reads & writes ITK's MatrixOffsetTransformBase transforms, with and without FixedParameters --- nibabel/transform/base.py | 2 +- nibabel/transform/linear.py | 28 ++++++++++++++++++++++++---- 2 files changed, 25 insertions(+), 5 deletions(-) diff --git a/nibabel/transform/base.py b/nibabel/transform/base.py index 9283092da4..f952935230 100644 --- a/nibabel/transform/base.py +++ b/nibabel/transform/base.py @@ -178,7 +178,7 @@ def _to_hdf5(self, x5_root): '''Serialize this object into the x5 file format''' raise NotImplementedError - def to_filename(self, filename): + def to_filename(self, filename, format='X5'): '''Store the transform in BIDS-Transforms HDF5 file format (.x5). ''' with h5py.File(filename, 'w') as out_file: diff --git a/nibabel/transform/linear.py b/nibabel/transform/linear.py index 7dad928784..aa9b428c82 100644 --- a/nibabel/transform/linear.py +++ b/nibabel/transform/linear.py @@ -14,6 +14,9 @@ from .base import TransformBase +LPS = np.diag([-1, -1, 1, 1]) + + class Affine(TransformBase): '''Represents linear transforms on image data''' __slots__ = ['_matrix'] @@ -160,6 +163,25 @@ def _to_hdf5(self, x5_root): refgrp.create_dataset('ndim', data=self.reference.ndim) refgrp['Type'] = 'image' + def to_filename(self, filename, fmt='X5'): + '''Store the transform in BIDS-Transforms HDF5 file format (.x5). + ''' + + if fmt.lower() in ['itk', 'ants', 'elastix', 'nifty']: + parameters = LPS.dot(self.matrix.dot(LPS)) + parameters = parameters[:3, :3].reshape(-1).tolist() + parameters[:3, 3].tolist() + with open(filename, 'w') as f: + itkfmt = """\ +#Insight Transform File V1.0 +#Transform 0 +Transform: MatrixOffsetTransformBase_double_3_3 +Parameters: {} +FixedParameters: 0 0 0\n""".format + f.write(itkfmt(' '.join([str(p) for p in parameters]))) + return filename + + return super(Affine, self).to_filename(filename, fmt=fmt) + def load(filename): ''' Load a linear transform ''' @@ -170,11 +192,9 @@ def load(filename): parameters = np.fromstring(itkxfm[3].split(':')[-1].strip(), dtype=float, sep=' ') offset = np.fromstring(itkxfm[4].split(':')[-1].strip(), dtype=float, sep=' ') if len(parameters) == 12: - matrix = np.vstack((parameters.reshape((4, 3)).T, [0, 0, 0, 1])) - matrix[:2, -1] *= -1.0 # Traslation - LPS to RAS + 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) - matrix = np.diag([-1, -1, 1, 1]).dot( - c_pos.dot(matrix.dot(c_neg.dot(np.diag([-1, -1, 1, 1]))))) + matrix = LPS.dot(c_pos.dot(matrix.dot(c_neg.dot(LPS)))) return Affine(matrix) From 8eb670dd4b6d2af1db817becdb11da3be1cd0dba Mon Sep 17 00:00:00 2001 From: oesteban Date: Fri, 15 Mar 2019 18:06:57 -0700 Subject: [PATCH 27/36] fix: print warnings to stderr, sty: minimal fixes --- nibabel/transform/linear.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/nibabel/transform/linear.py b/nibabel/transform/linear.py index aa9b428c82..5f66a3f916 100644 --- a/nibabel/transform/linear.py +++ b/nibabel/transform/linear.py @@ -8,6 +8,7 @@ ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## ''' Linear transforms ''' from __future__ import division, print_function, absolute_import +import sys import numpy as np from scipy import ndimage as ndi @@ -106,7 +107,8 @@ def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True, try: reference = self.reference except ValueError: - print('Warning: no reference space defined, using moving as reference') + print('Warning: no reference space defined, using moving as reference', + file=sys.stderr) reference = moving # Compose an index to index affine matrix @@ -134,7 +136,8 @@ def map_voxel(self, index, moving=None): try: reference = self.reference except ValueError: - print('Warning: no reference space defined, using moving as reference') + print('Warning: no reference space defined, using moving as reference', + file=sys.stderr) reference = moving else: if moving is None: @@ -170,13 +173,13 @@ def to_filename(self, filename, fmt='X5'): if fmt.lower() in ['itk', 'ants', 'elastix', 'nifty']: parameters = LPS.dot(self.matrix.dot(LPS)) parameters = parameters[:3, :3].reshape(-1).tolist() + parameters[:3, 3].tolist() - with open(filename, 'w') as f: - itkfmt = """\ + itkfmt = """\ #Insight Transform File V1.0 #Transform 0 Transform: MatrixOffsetTransformBase_double_3_3 Parameters: {} FixedParameters: 0 0 0\n""".format + with open(filename, 'w') as f: f.write(itkfmt(' '.join([str(p) for p in parameters]))) return filename From 799fb6e6ba8f7865968ff6344947d53b441996c7 Mon Sep 17 00:00:00 2001 From: oesteban Date: Tue, 19 Mar 2019 09:52:47 -0700 Subject: [PATCH 28/36] enh(affines): write out AFNI 12-parameters files --- nibabel/transform/linear.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/nibabel/transform/linear.py b/nibabel/transform/linear.py index 5f66a3f916..45a8e0dcb7 100644 --- a/nibabel/transform/linear.py +++ b/nibabel/transform/linear.py @@ -180,7 +180,15 @@ def to_filename(self, filename, fmt='X5'): Parameters: {} FixedParameters: 0 0 0\n""".format with open(filename, 'w') as f: - f.write(itkfmt(' '.join([str(p) for p in parameters]))) + f.write(itkfmt(' '.join(['%g' for p in parameters]))) + return filename + + if fmt.lower() == 'afni': + parameters = LPS.dot(self.matrix.dot(LPS)) + parameters = parameters[:3, :].reshape(-1).tolist() + np.savetxt(filename, np.atleast_2d(parameters), + delimiter='\t', header="""\ +3dvolreg matrices (DICOM-to-DICOM, row-by-row):""") return filename return super(Affine, self).to_filename(filename, fmt=fmt) From ce36f06c321e52fb925dc58ce496de378615b98b Mon Sep 17 00:00:00 2001 From: oesteban Date: Thu, 21 Mar 2019 14:43:30 -0700 Subject: [PATCH 29/36] enh(transforms): finish up a x5-to-fsl writer of affines --- nibabel/transform/linear.py | 58 ++++++++++++++++++++++++++++++------- 1 file changed, 47 insertions(+), 11 deletions(-) diff --git a/nibabel/transform/linear.py b/nibabel/transform/linear.py index 45a8e0dcb7..5c99867d34 100644 --- a/nibabel/transform/linear.py +++ b/nibabel/transform/linear.py @@ -12,6 +12,8 @@ import numpy as np from scipy import ndimage as ndi +from ..loadsave import load as loadimg +from ..affines import from_matvec, voxel_sizes from .base import TransformBase @@ -22,7 +24,7 @@ class Affine(TransformBase): '''Represents linear transforms on image data''' __slots__ = ['_matrix'] - def __init__(self, matrix=None): + def __init__(self, matrix=None, reference=None): '''Initialize a transform Parameters @@ -53,6 +55,9 @@ def __init__(self, matrix=None): assert self._matrix.shape[0] == self._matrix.shape[1], 'affine matrix is not square' super(Affine, self).__init__() + if reference: + self.reference = reference + @property def matrix(self): return self._matrix @@ -154,19 +159,14 @@ def map_voxel(self, index, moving=None): return tuple(matrix.dot(index)[:-1]) def _to_hdf5(self, x5_root): - print('to hdf5') x5_root.create_dataset('Transform', data=self._matrix[:3, :]) x5_root.create_dataset('Inverse', data=np.linalg.inv(self._matrix)[:3, :]) x5_root['Type'] = 'affine' if self._reference: - refgrp = x5_root.create_group('Reference') - refgrp.create_dataset('affine', data=self.reference.affine) - refgrp.create_dataset('shape', data=self.reference.shape) - refgrp.create_dataset('ndim', data=self.reference.ndim) - refgrp['Type'] = 'image' + self.reference._to_hdf5(x5_root.create_group('Reference')) - def to_filename(self, filename, fmt='X5'): + def to_filename(self, filename, fmt='X5', moving=None): '''Store the transform in BIDS-Transforms HDF5 file format (.x5). ''' @@ -191,12 +191,32 @@ def to_filename(self, filename, fmt='X5'): 3dvolreg matrices (DICOM-to-DICOM, row-by-row):""") 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(post.dot(self.matrix.dot(pre))) + np.savetxt(filename, mat, delimiter=' ', fmt='%g') + return filename return super(Affine, self).to_filename(filename, fmt=fmt) -def load(filename): +def load(filename, reference=None): ''' Load a linear transform ''' - from ..affines import from_matvec with open(filename) as itkfile: itkxfm = itkfile.read().splitlines() @@ -208,4 +228,20 @@ def load(filename): c_pos = from_matvec(np.eye(3), offset) matrix = LPS.dot(c_pos.dot(matrix.dot(c_neg.dot(LPS)))) - return Affine(matrix) + if reference and isinstance(reference, str): + reference = loadimg(reference) + + return Affine(matrix, reference) + + +def _fsl_aff_adapt(space): + """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) From 8e1bf44756a1568c9ab3da87c2cef84d0b22d35d Mon Sep 17 00:00:00 2001 From: oesteban Date: Thu, 21 Mar 2019 14:44:08 -0700 Subject: [PATCH 30/36] enh(transforms): improve generation of x5 structures of reference spaces --- nibabel/transform/base.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/nibabel/transform/base.py b/nibabel/transform/base.py index f952935230..cdc3557dde 100644 --- a/nibabel/transform/base.py +++ b/nibabel/transform/base.py @@ -77,6 +77,12 @@ def map_voxels(self, coordinates): return np.tensordot(np.linalg.inv(self._affine), coordinates, axes=1)[:3, ...] + def _to_hdf5(self, group): + group['Type'] = 'image' + group.create_dataset('affine', data=self.affine) + group.create_dataset('shape', data=self.shape) + group.create_dataset('ndim', data=self.ndim) + def __eq__(self, other): try: return np.allclose(self.affine, other.affine) and self.shape == other.shape @@ -178,7 +184,7 @@ def _to_hdf5(self, x5_root): '''Serialize this object into the x5 file format''' raise NotImplementedError - def to_filename(self, filename, format='X5'): + 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: From e2df7ccdb0fbc46f8b3629ddec5847fc64772bc5 Mon Sep 17 00:00:00 2001 From: oesteban Date: Fri, 22 Mar 2019 11:12:51 -0700 Subject: [PATCH 31/36] fix(transforms): string formating forgotten when writing ITK transforms --- nibabel/transform/linear.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nibabel/transform/linear.py b/nibabel/transform/linear.py index 5c99867d34..383864a1f0 100644 --- a/nibabel/transform/linear.py +++ b/nibabel/transform/linear.py @@ -180,7 +180,7 @@ def to_filename(self, filename, fmt='X5', moving=None): Parameters: {} FixedParameters: 0 0 0\n""".format with open(filename, 'w') as f: - f.write(itkfmt(' '.join(['%g' for p in parameters]))) + f.write(itkfmt(' '.join(['%g' % p for p in parameters]))) return filename if fmt.lower() == 'afni': From c5d10ed8e53e71d695c5fc5e90a3416f244fe9cb Mon Sep 17 00:00:00 2001 From: oesteban Date: Fri, 22 Mar 2019 11:42:59 -0700 Subject: [PATCH 32/36] wip(transforms): writing up some linear transform readers --- nibabel/transform/linear.py | 27 ++++++++++++++++----------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/nibabel/transform/linear.py b/nibabel/transform/linear.py index 383864a1f0..4f0b3f511e 100644 --- a/nibabel/transform/linear.py +++ b/nibabel/transform/linear.py @@ -215,22 +215,27 @@ def to_filename(self, filename, fmt='X5', moving=None): return super(Affine, self).to_filename(filename, fmt=fmt) -def load(filename, reference=None): +def load(filename, fmt='X5', reference=None): ''' Load a linear transform ''' - with open(filename) as itkfile: - itkxfm = itkfile.read().splitlines() - parameters = np.fromstring(itkxfm[3].split(':')[-1].strip(), dtype=float, sep=' ') - offset = np.fromstring(itkxfm[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) - matrix = LPS.dot(c_pos.dot(matrix.dot(c_neg.dot(LPS)))) + if fmt.lower() in ['itk', 'ants', 'elastix', 'nifty']: + with open(filename) as itkfile: + itkxfm = itkfile.read().splitlines() + + parameters = np.fromstring(itkxfm[3].split(':')[-1].strip(), dtype=float, sep=' ') + offset = np.fromstring(itkxfm[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) + matrix = LPS.dot(c_pos.dot(matrix.dot(c_neg.dot(LPS)))) + + # if fmt.lower() == 'afni': + # parameters = LPS.dot(self.matrix.dot(LPS)) + # parameters = parameters[:3, :].reshape(-1).tolist() if reference and isinstance(reference, str): reference = loadimg(reference) - return Affine(matrix, reference) From a18ce1d735b336dba909ba209949e46b94b5d220 Mon Sep 17 00:00:00 2001 From: oesteban Date: Fri, 22 Mar 2019 16:35:38 -0700 Subject: [PATCH 33/36] enh(transform): HDF5 - set values as attributes; generalize affines to N transforms (through time axis) --- nibabel/transform/base.py | 8 +- nibabel/transform/linear.py | 143 ++++++++++++++++++++++++------------ 2 files changed, 102 insertions(+), 49 deletions(-) diff --git a/nibabel/transform/base.py b/nibabel/transform/base.py index cdc3557dde..832e4a1164 100644 --- a/nibabel/transform/base.py +++ b/nibabel/transform/base.py @@ -78,10 +78,10 @@ def map_voxels(self, coordinates): coordinates, axes=1)[:3, ...] def _to_hdf5(self, group): - group['Type'] = 'image' + group.attrs['Type'] = 'image' + group.attrs['ndim'] = self.ndim group.create_dataset('affine', data=self.affine) group.create_dataset('shape', data=self.shape) - group.create_dataset('ndim', data=self.ndim) def __eq__(self, other): try: @@ -188,8 +188,8 @@ 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['Format'] = 'X5' - out_file['Version'] = np.uint16(1) + out_file.attrs['Format'] = 'X5' + out_file.attrs['Version'] = np.uint16(1) root = out_file.create_group('/0') self._to_hdf5(root) diff --git a/nibabel/transform/linear.py b/nibabel/transform/linear.py index 4f0b3f511e..d5b17dfef8 100644 --- a/nibabel/transform/linear.py +++ b/nibabel/transform/linear.py @@ -48,14 +48,19 @@ def __init__(self, matrix=None, reference=None): ''' if matrix is None: - matrix = np.eye(4) + matrix = [np.eye(4)] + + if np.ndim(matrix) == 2: + matrix = [matrix] self._matrix = np.array(matrix) - assert self._matrix.ndim in (2, 3), 'affine matrix should be 2D or 3D' - assert self._matrix.shape[0] == self._matrix.shape[1], 'affine matrix is not square' + assert self._matrix.ndim == 3, 'affine matrix should be 3D' + assert self._matrix.shape[-2] == self._matrix.shape[-1], 'affine matrix is not square' super(Affine, self).__init__() if reference: + if isinstance(reference, str): + reference = loadimg(reference) self.reference = reference @property @@ -116,28 +121,55 @@ def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True, 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 - matrix = np.linalg.inv(moving.affine).dot(self._matrix.dot(reference.affine)) - mdata = moving.get_data() - moved = ndi.affine_transform( - mdata, matrix=matrix[:mdata.ndim, :mdata.ndim], - offset=matrix[:mdata.ndim, mdata.ndim], - output_shape=reference.shape, order=order, mode=mode, - cval=cval, prefilter=prefilter) - - moved_image = moving.__class__(moved, reference.affine, moving.header) + 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, forward=True): + def map_point(self, coords, index=0, forward=True): coords = np.array(coords) - if coords.shape[0] == self._matrix.shape[0] - 1: + if coords.shape[0] == self._matrix[index].shape[0] - 1: coords = np.append(coords, [1]) - affine = self._matrix if forward else np.linalg.inv(self._matrix) + affine = self._matrix[index] if forward else np.linalg.inv(self._matrix[index]) return affine.dot(coords)[:-1] - def map_voxel(self, index, moving=None): + def map_voxel(self, index, nindex=0, moving=None): try: reference = self.reference except ValueError: @@ -152,16 +184,16 @@ def map_voxel(self, index, moving=None): raise ValueError('Reference and moving spaces are both undefined') index = np.array(index) - if index.shape[0] == self._matrix.shape[0] - 1: + if index.shape[0] == self._matrix[nindex].shape[0] - 1: index = np.append(index, [1]) - matrix = reference.affine.dot(self._matrix.dot(np.linalg.inv(moving.affine))) + 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): - x5_root.create_dataset('Transform', data=self._matrix[:3, :]) - x5_root.create_dataset('Inverse', data=np.linalg.inv(self._matrix)[:3, :]) - x5_root['Type'] = 'affine' + 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')) @@ -171,24 +203,26 @@ def to_filename(self, filename, fmt='X5', moving=None): ''' if fmt.lower() in ['itk', 'ants', 'elastix', 'nifty']: - parameters = LPS.dot(self.matrix.dot(LPS)) - parameters = parameters[:3, :3].reshape(-1).tolist() + parameters[:3, 3].tolist() - itkfmt = """\ -#Insight Transform File V1.0 -#Transform 0 + 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: {} +Parameters: {1} FixedParameters: 0 0 0\n""".format - with open(filename, 'w') as f: - f.write(itkfmt(' '.join(['%g' % p for p in parameters]))) + f.write(itkfmt(i, ' '.join(['%g' % p for p in parameters]))) return filename if fmt.lower() == 'afni': - parameters = LPS.dot(self.matrix.dot(LPS)) - parameters = parameters[:3, :].reshape(-1).tolist() - np.savetxt(filename, np.atleast_2d(parameters), - delimiter='\t', header="""\ -3dvolreg matrices (DICOM-to-DICOM, row-by-row):""") + parameters = np.swapaxes(LPS.dot(self.matrix.dot(LPS)), 0, 1) + parameters = parameters[:, :3, :].reshape((self.matrix.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': @@ -209,8 +243,15 @@ def to_filename(self, filename, fmt='X5', moving=None): moving.affine))) # Compose FSL transform - mat = np.linalg.inv(post.dot(self.matrix.dot(pre))) - np.savetxt(filename, mat, delimiter=' ', fmt='%g') + 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) @@ -222,17 +263,29 @@ def load(filename, fmt='X5', reference=None): with open(filename) as itkfile: itkxfm = itkfile.read().splitlines() - parameters = np.fromstring(itkxfm[3].split(':')[-1].strip(), dtype=float, sep=' ') - offset = np.fromstring(itkxfm[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) - matrix = LPS.dot(c_pos.dot(matrix.dot(c_neg.dot(LPS)))) - - # if fmt.lower() == 'afni': + 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) From 9902f505fc36984da38bb1572e3ec8594fd7f46a Mon Sep 17 00:00:00 2001 From: oesteban Date: Thu, 19 Sep 2019 11:41:45 -0700 Subject: [PATCH 34/36] enh: apply some early comments from @effigies. --- .gitignore | 1 - nibabel/transform/__init__.py | 4 ++-- nibabel/transform/base.py | 34 +++++++++++++----------------- nibabel/transform/linear.py | 38 +++++++++++++++++----------------- nibabel/transform/nonlinear.py | 30 ++++++++++++++------------- 5 files changed, 51 insertions(+), 56 deletions(-) diff --git a/.gitignore b/.gitignore index bbd0a11b43..df018f0ead 100644 --- a/.gitignore +++ b/.gitignore @@ -84,4 +84,3 @@ Thumbs.db doc/source/reference venv/ .buildbot.patch -nibabel/maths/bspline.c diff --git a/nibabel/transform/__init__.py b/nibabel/transform/__init__.py index 323a353cc7..be52d50c42 100644 --- a/nibabel/transform/__init__.py +++ b/nibabel/transform/__init__.py @@ -6,7 +6,8 @@ # copyright and license terms. # ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## -"""Geometric transforms +""" +Geometric transforms. .. currentmodule:: nibabel.transform @@ -15,7 +16,6 @@ transform """ -from __future__ import absolute_import from .linear import Affine from .nonlinear import DeformationFieldTransform diff --git a/nibabel/transform/base.py b/nibabel/transform/base.py index 832e4a1164..ada8c6a799 100644 --- a/nibabel/transform/base.py +++ b/nibabel/transform/base.py @@ -6,8 +6,7 @@ # copyright and license terms. # ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## -''' Common interface for transforms ''' -from __future__ import division, print_function, absolute_import +"""Common interface for transforms.""" import numpy as np import h5py @@ -15,7 +14,8 @@ class ImageSpace(object): - '''Class to represent spaces of gridded data (images)''' + """Class to represent spaces of gridded data (images).""" + __slots__ = ['_affine', '_shape', '_ndim', '_ndindex', '_coords', '_nvox', '_inverse'] @@ -92,9 +92,8 @@ def __eq__(self, other): class TransformBase(object): - ''' - Abstract image class to represent transforms - ''' + """Abstract image class to represent transforms.""" + __slots__ = ['_reference'] def __init__(self): @@ -117,11 +116,11 @@ def ndim(self): def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True, output_dtype=None): - '''Resample the moving image in reference space + """ + Resample the moving image in reference space. Parameters ---------- - moving : `spatialimage` The image object containing the data to be resampled in reference space @@ -144,16 +143,14 @@ def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True, 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.header.get_data_dtype() + output_dtype = moving_data.dtype - moving_data = moving.get_data() moved = ndi.geometric_transform( moving_data, mapping=self.map_voxel, @@ -171,22 +168,19 @@ def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True, return moved_image def map_point(self, coords): - '''Find the coordinates in moving space corresponding to the - input reference coordinates''' + """Apply y = f(x), where x is the argument `coords`.""" raise NotImplementedError def map_voxel(self, index, moving=None): - '''Find the voxel indices in the moving image corresponding to the - input reference voxel''' + """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''' + """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). - ''' + """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) diff --git a/nibabel/transform/linear.py b/nibabel/transform/linear.py index d5b17dfef8..f33f6342a1 100644 --- a/nibabel/transform/linear.py +++ b/nibabel/transform/linear.py @@ -6,8 +6,7 @@ # copyright and license terms. # ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## -''' Linear transforms ''' -from __future__ import division, print_function, absolute_import +"""Linear transforms.""" import sys import numpy as np from scipy import ndimage as ndi @@ -21,15 +20,16 @@ class Affine(TransformBase): - '''Represents linear transforms on image data''' + """Represents linear transforms on image data.""" + __slots__ = ['_matrix'] def __init__(self, matrix=None, reference=None): - '''Initialize a transform + """ + Initialize a linear transform. Parameters ---------- - matrix : ndarray The inverse coordinate transformation matrix **in physical coordinates**, mapping coordinates from *reference* space @@ -38,7 +38,6 @@ def __init__(self, matrix=None, reference=None): 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], @@ -46,7 +45,8 @@ def __init__(self, matrix=None, reference=None): [0, 0, 1, 0], [0, 0, 0, 1]]) - ''' + """ + super(Affine, self).__init__() if matrix is None: matrix = [np.eye(4)] @@ -56,7 +56,6 @@ def __init__(self, matrix=None, reference=None): 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' - super(Affine, self).__init__() if reference: if isinstance(reference, str): @@ -69,11 +68,11 @@ def matrix(self): def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True, output_dtype=None): - '''Resample the moving image in reference space + """ + Resample the moving image in reference space. Parameters ---------- - moving : `spatialimage` The image object containing the data to be resampled in reference space @@ -96,21 +95,19 @@ def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True, 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() @@ -163,6 +160,7 @@ def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True, 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]) @@ -170,6 +168,7 @@ def map_point(self, coords, index=0, forward=True): 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: @@ -191,6 +190,7 @@ def map_voxel(self, index, nindex=0, moving=None): 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)) @@ -199,9 +199,7 @@ def _to_hdf5(self, x5_root): 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). - ''' - + """Store the transform in BIDS-Transforms HDF5 file format (.x5).""" if fmt.lower() in ['itk', 'ants', 'elastix', 'nifty']: with open(filename, 'w') as f: f.write('#Insight Transform File V1.0\n') @@ -257,8 +255,7 @@ def to_filename(self, filename, fmt='X5', moving=None): def load(filename, fmt='X5', reference=None): - ''' Load a linear transform ''' - + """Load a linear transform.""" if fmt.lower() in ['itk', 'ants', 'elastix', 'nifty']: with open(filename) as itkfile: itkxfm = itkfile.read().splitlines() @@ -293,7 +290,10 @@ def load(filename, fmt='X5', reference=None): def _fsl_aff_adapt(space): - """Calculates a matrix to convert from the original RAS image + """ + 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 diff --git a/nibabel/transform/nonlinear.py b/nibabel/transform/nonlinear.py index a538628d99..e607bd144f 100644 --- a/nibabel/transform/nonlinear.py +++ b/nibabel/transform/nonlinear.py @@ -6,8 +6,7 @@ # copyright and license terms. # ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## -''' Common interface for transforms ''' -from __future__ import division, print_function, absolute_import +"""Nonlinear transforms.""" import numpy as np from scipy import ndimage as ndi # from gridbspline.maths import cubic @@ -19,14 +18,13 @@ class DeformationFieldTransform(TransformBase): - '''Represents a dense field of displacements (one vector per voxel)''' + """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 - ''' + """Create a dense deformation field transform.""" super(DeformationFieldTransform, self).__init__() self._field = field.get_data() @@ -103,11 +101,11 @@ def _cache_moving(self, moving): 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') @@ -118,16 +116,18 @@ def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True, >>> 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] @@ -154,11 +154,13 @@ def map_coordinates(self, coordinates, order=3, mode='mirror', cval=0.0, 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''' + """Create a smooth deformation field using B-Spline basis.""" super(BSplineFieldTransform, self).__init__() self._order = order self.reference = reference @@ -216,16 +218,16 @@ def _interp_transform(self, coords): return self.reference.inverse.dot(np.hstack((coords, 1)))[:3] def map_voxel(self, index, moving=None): - '''Find the corresponding coordinates for a voxel in reference space''' + """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') @@ -238,7 +240,7 @@ def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True, >>> 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) From fe74efb867e761cf4417cf1c9ca34ec3427e1c3e Mon Sep 17 00:00:00 2001 From: oesteban Date: Thu, 26 Sep 2019 15:42:40 -0700 Subject: [PATCH 35/36] fix: add a first battery of tests Added some initial tests to check the writing of transforms to ITK and FSL formats. --- nibabel/tests/data/affine-LAS-itk.tfm | 1 + nibabel/tests/data/affine-LAS.fsl | 1 + nibabel/tests/data/affine-LPS-itk.tfm | 5 ++ nibabel/tests/data/affine-LPS.fsl | 4 ++ nibabel/tests/data/affine-RAS-itk.tfm | 5 ++ nibabel/tests/data/affine-RAS.fsl | 4 ++ nibabel/tests/data/affine-oblique-itk.tfm | 5 ++ nibabel/tests/data/affine-oblique.fsl | 4 ++ nibabel/tests/data/someones_anatomy.nii.gz | 1 + nibabel/tests/test_transform.py | 71 ++++++++++++++++++++++ nibabel/transform/base.py | 11 +++- nibabel/transform/linear.py | 37 ++++++++++- 12 files changed, 145 insertions(+), 4 deletions(-) create mode 120000 nibabel/tests/data/affine-LAS-itk.tfm create mode 120000 nibabel/tests/data/affine-LAS.fsl create mode 100644 nibabel/tests/data/affine-LPS-itk.tfm create mode 100644 nibabel/tests/data/affine-LPS.fsl create mode 100644 nibabel/tests/data/affine-RAS-itk.tfm create mode 100644 nibabel/tests/data/affine-RAS.fsl create mode 100644 nibabel/tests/data/affine-oblique-itk.tfm create mode 100644 nibabel/tests/data/affine-oblique.fsl create mode 120000 nibabel/tests/data/someones_anatomy.nii.gz create mode 100644 nibabel/tests/test_transform.py 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.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.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.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_transform.py b/nibabel/tests/test_transform.py new file mode 100644 index 0000000000..36756d9e13 --- /dev/null +++ b/nibabel/tests/test_transform.py @@ -0,0 +1,71 @@ +# 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)) + + with InTemporaryDirectory(): + xfm.to_filename('M.tfm', fmt='itk') + xfm.to_filename('M.fsl', fmt='fsl') + + nb_itk = nbl.load('M.tfm', fmt='itk') + nb_fsl = np.loadtxt('M.fsl') + + assert_equal(itk, nb_itk) + assert_almost_equal(fsl, nb_fsl) + +# Create version not aligned to canonical diff --git a/nibabel/transform/base.py b/nibabel/transform/base.py index ada8c6a799..ec9d4e8b6c 100644 --- a/nibabel/transform/base.py +++ b/nibabel/transform/base.py @@ -12,6 +12,8 @@ from scipy import ndimage as ndi +EQUALITY_TOL = 1e-5 + class ImageSpace(object): """Class to represent spaces of gridded data (images).""" @@ -85,7 +87,8 @@ def _to_hdf5(self, group): def __eq__(self, other): try: - return np.allclose(self.affine, other.affine) and self.shape == other.shape + return np.allclose( + self.affine, other.affine, rtol=EQUALITY_TOL) and self.shape == other.shape except AttributeError: pass return False @@ -99,6 +102,12 @@ class TransformBase(object): 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''' diff --git a/nibabel/transform/linear.py b/nibabel/transform/linear.py index f33f6342a1..8a0fc15b05 100644 --- a/nibabel/transform/linear.py +++ b/nibabel/transform/linear.py @@ -10,13 +10,16 @@ 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 +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): @@ -64,6 +67,7 @@ def __init__(self, matrix=None, reference=None): @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, @@ -217,8 +221,35 @@ def to_filename(self, filename, fmt='X5', moving=None): return filename if fmt.lower() == 'afni': - parameters = np.swapaxes(LPS.dot(self.matrix.dot(LPS)), 0, 1) - parameters = parameters[:, :3, :].reshape((self.matrix.shape[0], -1)) + 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 any(obliquity(self.reference.affine) * 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=True) + pre = M.dot(np.linalg.inv(A)) + + if not moving: + moving = self.reference + + if moving and any(obliquity(moving.affine) * 180 / pi > OBLIQUITY_THRESHOLD_DEG): + print('Moving affine axes are oblique.') + M = moving.affine + A = shape_zoom_affine(moving.shape, + voxel_sizes(M), x_flip=True) + post = M.dot(np.linalg.inv(A)) + + # swapaxes is necessary, as axis 0 encodes series of transforms + T = np.swapaxes(post.dot(self.matrix.copy().dot(pre)), 0, 1) + parameters = np.swapaxes(post.dot(T.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 From 596994cbfc9784e39ffbe9fa8e27387d46c568ee Mon Sep 17 00:00:00 2001 From: oesteban Date: Fri, 27 Sep 2019 11:37:14 -0700 Subject: [PATCH 36/36] enh: add tests for affines stored in AFNI format (non-oblique images) --- nibabel/tests/data/affine-LAS.afni | 1 + nibabel/tests/data/affine-LPS.afni | 1 + nibabel/tests/data/affine-RAS.afni | 2 ++ nibabel/tests/test_transform.py | 4 ++++ nibabel/transform/linear.py | 21 ++++++++++----------- nibabel/volumeutils.py | 8 +++++++- 6 files changed, 25 insertions(+), 12 deletions(-) create mode 120000 nibabel/tests/data/affine-LAS.afni create mode 120000 nibabel/tests/data/affine-LPS.afni create mode 100644 nibabel/tests/data/affine-RAS.afni 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-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-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/test_transform.py b/nibabel/tests/test_transform.py index 36756d9e13..f8160e2ee1 100644 --- a/nibabel/tests/test_transform.py +++ b/nibabel/tests/test_transform.py @@ -57,15 +57,19 @@ def test_affines_save(image_orientation): 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/linear.py b/nibabel/transform/linear.py index 8a0fc15b05..84017b203a 100644 --- a/nibabel/transform/linear.py +++ b/nibabel/transform/linear.py @@ -204,7 +204,7 @@ def _to_hdf5(self, x5_root): 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', 'nifty']: + if fmt.lower() in ['itk', 'ants', 'elastix']: with open(filename, 'w') as f: f.write('#Insight Transform File V1.0\n') @@ -229,26 +229,25 @@ def to_filename(self, filename, fmt='X5', moving=None): T = self.matrix.copy() pre = LPS post = LPS - if any(obliquity(self.reference.affine) * 180 / pi > OBLIQUITY_THRESHOLD_DEG): + 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=True) - pre = M.dot(np.linalg.inv(A)) + 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 any(obliquity(moving.affine) * 180 / pi > OBLIQUITY_THRESHOLD_DEG): + if moving and obliquity(moving.affine).min() * 180 / pi > OBLIQUITY_THRESHOLD_DEG: print('Moving affine axes are oblique.') - M = moving.affine - A = shape_zoom_affine(moving.shape, - voxel_sizes(M), x_flip=True) - post = M.dot(np.linalg.inv(A)) + 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 - T = np.swapaxes(post.dot(self.matrix.copy().dot(pre)), 0, 1) - parameters = np.swapaxes(post.dot(T.dot(pre)), 0, 1) + 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') 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)