Merge pull request #544 from effigies/improve_as_closest_canonical

RF: Move reorientation into image classes

Author: effigies

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
 ---------

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

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
 ===============

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):

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

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)

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))