Merge pull request #2572 from effigies/enh/reorient

ENH: Reorient interface

Author: oesteban

nipype/interfaces/image.py

@@ -0,0 +1,145 @@
+# -*- coding: utf-8 -*-
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+
+import numpy as np
+import nibabel as nb
+
+from ..utils.filemanip import fname_presuffix
+from .base import (SimpleInterface, TraitedSpec, BaseInterfaceInputSpec,
+                   traits, File)
+
+_axes = ('RL', 'AP', 'SI')
+_orientations = tuple(
+    ''.join((x[i], y[j], z[k]))
+    for x in _axes for y in _axes for z in _axes
+    if x != y != z != x
+    for i in (0, 1) for j in (0, 1) for k in (0, 1))
+
+
+class ReorientInputSpec(BaseInterfaceInputSpec):
+    in_file = File(exists=True, mandatory=True, desc='Input image')
+    orientation = traits.Enum(_orientations, usedefault=True,
+                              desc='Target axis orientation')
+
+
+class ReorientOutputSpec(TraitedSpec):
+    out_file = File(exists=True, desc='Reoriented image')
+    transform = File(exists=True,
+                     desc='Affine transform from input orientation to output')
+
+
+class Reorient(SimpleInterface):
+    """Conform an image to a given orientation
+
+    Flips and reorder the image data array so that the axes match the
+    directions indicated in ``orientation``.
+    The default ``RAS`` orientation corresponds to the first axis being ordered
+    from left to right, the second axis from posterior to anterior, and the
+    third axis from inferior to superior.
+
+    For oblique images, the original orientation is considered to be the
+    closest plumb orientation.
+
+    No resampling is performed, and thus the output image is not de-obliqued
+    or registered to any other image or template.
+
+    The effective transform is calculated from the original affine matrix to
+    the reoriented affine matrix.
+
+    Examples
+    --------
+
+    If an image is not reoriented, the original file is not modified
+
+    >>> import numpy as np
+    >>> from nipype.interfaces.image import Reorient
+    >>> reorient = Reorient(orientation='LPS')
+    >>> reorient.inputs.in_file = 'segmentation0.nii.gz'
+    >>> res = reorient.run()
+    >>> res.outputs.out_file
+    'segmentation0.nii.gz'
+
+    >>> print(np.loadtxt(res.outputs.transform))
+    [[1. 0. 0. 0.]
+     [0. 1. 0. 0.]
+     [0. 0. 1. 0.]
+     [0. 0. 0. 1.]]
+
+    >>> reorient.inputs.orientation = 'RAS'
+    >>> res = reorient.run()
+    >>> res.outputs.out_file  # doctest: +ELLIPSIS
+    '.../segmentation0_ras.nii.gz'
+
+    >>> print(np.loadtxt(res.outputs.transform))
+    [[-1.  0.  0. 60.]
+     [ 0. -1.  0. 72.]
+     [ 0.  0.  1.  0.]
+     [ 0.  0.  0.  1.]]
+
+    """
+    input_spec = ReorientInputSpec
+    output_spec = ReorientOutputSpec
+
+    def _run_interface(self, runtime):
+        from nibabel.orientations import (
+            axcodes2ornt, ornt_transform, inv_ornt_aff)
+
+        fname = self.inputs.in_file
+        orig_img = nb.load(fname)
+
+        # Find transform from current (approximate) orientation to
+        # target, in nibabel orientation matrix and affine forms
+        orig_ornt = nb.io_orientation(orig_img.affine)
+        targ_ornt = axcodes2ornt(self.inputs.orientation)
+        transform = ornt_transform(orig_ornt, targ_ornt)
+        affine_xfm = inv_ornt_aff(transform, orig_img.shape)
+
+        # Check can be eliminated when minimum nibabel version >= 2.2
+        if hasattr(orig_img, 'as_reoriented'):
+            reoriented = orig_img.as_reoriented(transform)
+        else:
+            reoriented = _as_reoriented_backport(orig_img, transform)
+
+        # Image may be reoriented
+        if reoriented is not orig_img:
+            suffix = '_' + self.inputs.orientation.lower()
+            out_name = fname_presuffix(fname, suffix=suffix,
+                                       newpath=runtime.cwd)
+            reoriented.to_filename(out_name)
+        else:
+            out_name = fname
+
+        mat_name = fname_presuffix(fname, suffix='.mat',
+                                   newpath=runtime.cwd, use_ext=False)
+        np.savetxt(mat_name, affine_xfm, fmt='%.08f')
+
+        self._results['out_file'] = out_name
+        self._results['transform'] = mat_name
+
+        return runtime
+
+
+def _as_reoriented_backport(img, ornt):
+    """Backport of img.as_reoriented as of nibabel 2.2.0"""
+    from nibabel.orientations import inv_ornt_aff
+    if np.array_equal(ornt, [[0, 1], [1, 1], [2, 1]]):
+        return img
+
+    t_arr = nb.apply_orientation(img.get_data(), ornt)
+    new_aff = img.affine.dot(inv_ornt_aff(ornt, img.shape))
+    reoriented = img.__class__(t_arr, new_aff, img.header)
+
+    if isinstance(reoriented, nb.Nifti1Pair):
+        # Also apply the transform to the dim_info fields
+        new_dim = list(reoriented.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]
+
+        reoriented.header.set_dim_info(*new_dim)
+
+    return reoriented

nipype/interfaces/tests/test_auto_Reorient.py

@@ -0,0 +1,30 @@
+# AUTO-GENERATED by tools/checkspecs.py - DO NOT EDIT
+from __future__ import unicode_literals
+from ..image import Reorient
+
+
+def test_Reorient_inputs():
+    input_map = dict(
+        ignore_exception=dict(
+            deprecated='1.0.0',
+            nohash=True,
+            usedefault=True,
+        ),
+        in_file=dict(mandatory=True, ),
+        orientation=dict(usedefault=True, ),
+    )
+    inputs = Reorient.input_spec()
+
+    for key, metadata in list(input_map.items()):
+        for metakey, value in list(metadata.items()):
+            assert getattr(inputs.traits()[key], metakey) == value
+def test_Reorient_outputs():
+    output_map = dict(
+        out_file=dict(),
+        transform=dict(),
+    )
+    outputs = Reorient.output_spec()
+
+    for key, metadata in list(output_map.items()):
+        for metakey, value in list(metadata.items()):
+            assert getattr(outputs.traits()[key], metakey) == value

nipype/interfaces/tests/test_image.py

@@ -0,0 +1,64 @@
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+import numpy as np
+import nibabel as nb
+import pytest
+
+from nibabel.orientations import axcodes2ornt, ornt_transform
+
+from ..image import _as_reoriented_backport, _orientations
+from ... import LooseVersion
+
+nibabel22 = LooseVersion(nb.__version__) >= LooseVersion('2.2.0')
+
+
+@pytest.mark.skipif(not nibabel22,
+                    reason="Old nibabel - can't directly compare")
+def test_reorientation_backport():
+    pixdims = ((1, 1, 1), (2, 2, 3))
+    data = np.random.normal(size=(17, 18, 19, 2))
+
+    for pixdim in pixdims:
+        # Generate a randomly rotated affine
+        angles = np.random.uniform(-np.pi, np.pi, 3) * [1, 0.5, 1]
+        rot = nb.eulerangles.euler2mat(*angles)
+        scale = np.diag(pixdim)
+        translation = np.array((17, 18, 19)) / 2
+        affine = nb.affines.from_matvec(rot.dot(scale), translation)
+
+        # Create image
+        img = nb.Nifti1Image(data, affine)
+        dim_info = {'freq': 0, 'phase': 1, 'slice': 2}
+        img.header.set_dim_info(**dim_info)
+
+        # Find a random, non-identity transform
+        targ_ornt = orig_ornt = nb.io_orientation(affine)
+        while np.array_equal(targ_ornt, orig_ornt):
+            new_code = np.random.choice(_orientations)
+            targ_ornt = axcodes2ornt(new_code)
+
+        identity = ornt_transform(orig_ornt, orig_ornt)
+        transform = ornt_transform(orig_ornt, targ_ornt)
+
+        # Identity transform returns exact image
+        assert img.as_reoriented(identity) is img
+        assert _as_reoriented_backport(img, identity) is img
+
+        reoriented_a = img.as_reoriented(transform)
+        reoriented_b = _as_reoriented_backport(img, transform)
+
+        flips_only = img.shape == reoriented_a.shape
+
+        # Reorientation changes affine and data array
+        assert not np.allclose(img.affine, reoriented_a.affine)
+        assert not (flips_only and
+                    np.allclose(img.get_data(), reoriented_a.get_data()))
+        # Dimension info changes iff axes are reordered
+        assert flips_only == np.array_equal(img.header.get_dim_info(),
+                                            reoriented_a.header.get_dim_info())
+
+        # Both approaches produce equivalent images
+        assert np.allclose(reoriented_a.affine, reoriented_b.affine)
+        assert np.array_equal(reoriented_a.get_data(), reoriented_b.get_data())
+        assert np.array_equal(reoriented_a.header.get_dim_info(),
+                              reoriented_b.header.get_dim_info())