diff --git a/Changelog b/Changelog index e676d25b3d..fdd3ef48b4 100644 --- a/Changelog +++ b/Changelog @@ -49,6 +49,8 @@ Enhancements ArrayProxy reshaping (pr/521) (CM) * Allow unknown NIfTI intent codes, add FSL codes (pr/528) (Paul McCarthy) * Improve error handling for ``img.__getitem__`` (pr/533) (Ariel Rokem) +* Delegate reorientation to SpatialImage classes (pr/544) (Mark Hymers, CM, + reviewed by MB) Bug fixes --------- diff --git a/Makefile b/Makefile index 326212a3c9..2190f815fc 100644 --- a/Makefile +++ b/Makefile @@ -288,4 +288,4 @@ rm-orig: # Remove .orig temporary diff files generated by git find . -name "*.orig" -print | grep -v "fsaverage" | xargs rm -.PHONY: orig-src pylint +.PHONY: orig-src pylint all build diff --git a/doc/source/index.rst b/doc/source/index.rst index 1c6bc503d6..9ce39b4254 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -77,6 +77,7 @@ contributed code and discussion (in rough order of appearance): * Paul McCarthy * Fernando Pérez García * Venky Reddy +* Mark Hymers License reprise =============== diff --git a/nibabel/funcs.py b/nibabel/funcs.py index 6225a40573..240b20f802 100644 --- a/nibabel/funcs.py +++ b/nibabel/funcs.py @@ -10,8 +10,7 @@ ''' Processor functions for images ''' import numpy as np -from .orientations import (io_orientation, inv_ornt_aff, - apply_orientation, OrientationError) +from .orientations import io_orientation, OrientationError from .loadsave import load @@ -206,25 +205,14 @@ def as_closest_canonical(img, enforce_diag=False): already has the correct data ordering, we just return `img` unmodified. ''' - aff = img.affine - ornt = io_orientation(aff) - if np.all(ornt == [[0, 1], - [1, 1], - [2, 1]]): # canonical already - # however, the affine may not be diagonal - if enforce_diag and not _aff_is_diag(aff): - raise OrientationError('Transformed affine is not diagonal') - return img - shape = img.shape - t_aff = inv_ornt_aff(ornt, shape) - out_aff = np.dot(aff, t_aff) - # check if we are going to end up with something diagonal - if enforce_diag and not _aff_is_diag(aff): + # Get the image class to transform the data for us + img = img.as_reoriented(io_orientation(img.affine)) + + # however, the affine may not be diagonal + if enforce_diag and not _aff_is_diag(img.affine): raise OrientationError('Transformed affine is not diagonal') - # we need to transform the data - arr = img.get_data() - t_arr = apply_orientation(arr, ornt) - return img.__class__(t_arr, out_aff, img.header) + + return img def _aff_is_diag(aff): diff --git a/nibabel/nifti1.py b/nibabel/nifti1.py index 88702594f5..51f5892f1b 100644 --- a/nibabel/nifti1.py +++ b/nibabel/nifti1.py @@ -1274,7 +1274,8 @@ def set_dim_info(self, freq=None, phase=None, slice=None): This is stored in one byte in the header ''' for inp in (freq, phase, slice): - if inp not in (None, 0, 1, 2): + # Don't use == on None to avoid a FutureWarning in python3 + if inp is not None and inp not in (0, 1, 2): raise HeaderDataError('Inputs must be in [None, 0, 1, 2]') info = 0 if freq is not None: @@ -1968,6 +1969,40 @@ def set_sform(self, affine, code=None, **kwargs): else: self._affine[:] = self._header.get_best_affine() + def as_reoriented(self, ornt): + """Apply an orientation change and return a new image + + If ornt is identity transform, return the original image, unchanged + + Parameters + ---------- + ornt : (n,2) orientation array + orientation transform. ``ornt[N,1]` is flip of axis N of the + array implied by `shape`, where 1 means no flip and -1 means + flip. For example, if ``N==0`` and ``ornt[0,1] == -1``, and + there's an array ``arr`` of shape `shape`, the flip would + correspond to the effect of ``np.flipud(arr)``. ``ornt[:,0]`` is + the transpose that needs to be done to the implied array, as in + ``arr.transpose(ornt[:,0])`` + """ + img = super(Nifti1Pair, self).as_reoriented(ornt) + + if img is self: + return img + + # Also apply the transform to the dim_info fields + new_dim = list(img.header.get_dim_info()) + for idx, value in enumerate(new_dim): + # For each value, leave as None if it was that way, + # otherwise check where we have mapped it to + if value is None: + continue + new_dim[idx] = np.where(ornt[:, 0] == idx)[0] + + img.header.set_dim_info(*new_dim) + + return img + class Nifti1Image(Nifti1Pair): """ Class for single file NIfTI1 format image diff --git a/nibabel/spatialimages.py b/nibabel/spatialimages.py index 86265697b4..2554600649 100644 --- a/nibabel/spatialimages.py +++ b/nibabel/spatialimages.py @@ -141,6 +141,7 @@ from .viewers import OrthoSlicer3D from .volumeutils import shape_zoom_affine from .deprecated import deprecate_with_version +from .orientations import apply_orientation, inv_ornt_aff class HeaderDataError(Exception): @@ -483,3 +484,34 @@ def orthoview(self): """ return OrthoSlicer3D(self.dataobj, self.affine, title=self.get_filename()) + + + def as_reoriented(self, ornt): + """Apply an orientation change and return a new image + + If ornt is identity transform, return the original image, unchanged + + Parameters + ---------- + ornt : (n,2) orientation array + orientation transform. ``ornt[N,1]` is flip of axis N of the + array implied by `shape`, where 1 means no flip and -1 means + flip. For example, if ``N==0`` and ``ornt[0,1] == -1``, and + there's an array ``arr`` of shape `shape`, the flip would + correspond to the effect of ``np.flipud(arr)``. ``ornt[:,0]`` is + the transpose that needs to be done to the implied array, as in + ``arr.transpose(ornt[:,0])`` + + Notes + ----- + Subclasses should override this if they have additional requirements + when re-orienting an image. + """ + + if np.array_equal(ornt, [[0, 1], [1, 1], [2, 1]]): + return self + + t_arr = apply_orientation(self.get_data(), ornt) + new_aff = self.affine.dot(inv_ornt_aff(ornt, self.shape)) + + return self.__class__(t_arr, new_aff, self.header) diff --git a/nibabel/tests/test_funcs.py b/nibabel/tests/test_funcs.py index 1841dbdd9a..6032c08672 100644 --- a/nibabel/tests/test_funcs.py +++ b/nibabel/tests/test_funcs.py @@ -12,6 +12,7 @@ import numpy as np from ..funcs import concat_images, as_closest_canonical, OrientationError +from ..analyze import AnalyzeImage from ..nifti1 import Nifti1Image from ..loadsave import save @@ -128,17 +129,40 @@ def test_concat(): def test_closest_canonical(): - arr = np.arange(24).reshape((2, 3, 4, 1)) - # no funky stuff, returns same thing + # Use 32-bit data so that the AnalyzeImage class doesn't complain + arr = np.arange(24).reshape((2, 3, 4, 1)).astype(np.int32) + + # Test with an AnalyzeImage first + img = AnalyzeImage(arr, np.eye(4)) + xyz_img = as_closest_canonical(img) + assert_true(img is xyz_img) + + # And a case where the Analyze image has to be flipped + img = AnalyzeImage(arr, np.diag([-1, 1, 1, 1])) + xyz_img = as_closest_canonical(img) + assert_false(img is xyz_img) + out_arr = xyz_img.get_data() + assert_array_equal(out_arr, np.flipud(arr)) + + # Now onto the NIFTI cases (where dim_info also has to be updated) + + # No funky stuff, returns same thing img = Nifti1Image(arr, np.eye(4)) + # set freq/phase/slice dim so that we can check that we + # re-order them properly + img.header.set_dim_info(0, 1, 2) xyz_img = as_closest_canonical(img) assert_true(img is xyz_img) + # a axis flip img = Nifti1Image(arr, np.diag([-1, 1, 1, 1])) + img.header.set_dim_info(0, 1, 2) xyz_img = as_closest_canonical(img) assert_false(img is xyz_img) + assert_true(img.header.get_dim_info() == xyz_img.header.get_dim_info()) out_arr = xyz_img.get_data() assert_array_equal(out_arr, np.flipud(arr)) + # no error for enforce_diag in this case xyz_img = as_closest_canonical(img, True) # but there is if the affine is not diagonal @@ -150,3 +174,22 @@ def test_closest_canonical(): assert_true(img is xyz_img) # it's still not diagnonal assert_raises(OrientationError, as_closest_canonical, img, True) + + # an axis swap + aff = np.diag([1, 0, 0, 1]) + aff[1, 2] = 1; aff[2, 1] = 1 + img = Nifti1Image(arr, aff) + img.header.set_dim_info(0, 1, 2) + + xyz_img = as_closest_canonical(img) + assert_false(img is xyz_img) + # Check both the original and new objects + assert_true(img.header.get_dim_info() == (0, 1, 2)) + assert_true(xyz_img.header.get_dim_info() == (0, 2, 1)) + out_arr = xyz_img.get_data() + assert_array_equal(out_arr, np.transpose(arr, (0, 2, 1, 3))) + + # same axis swap but with None dim info (except for slice dim) + img.header.set_dim_info(None, None, 2) + xyz_img = as_closest_canonical(img) + assert_true(xyz_img.header.get_dim_info() == (None, None, 1))