ENH,WIP: New transforms module

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

nibabel/affines.py

@@ -2,7 +2,6 @@
 # vi: set ft=python sts=4 ts=4 sw=4 et:
 """ Utility routines for working with points and affine transforms
 """
-
 import numpy as np
 
 from functools import reduce
@@ -296,3 +295,31 @@ def voxel_sizes(affine):
     """
     top_left = affine[:-1, :-1]
     return np.sqrt(np.sum(top_left ** 2, axis=0))
+
+
+def obliquity(affine):
+    r"""
+    Estimate the *obliquity* an affine's axes represent.
+
+    The term *obliquity* is defined here as the rotation of those axes with
+    respect to the cardinal axes.
+    This implementation is inspired by `AFNI's implementation
+    <https://github.com/afni/afni/blob/b6a9f7a21c1f3231ff09efbd861f8975ad48e525/src/thd_coords.c#L660-L698>`_.
+    For further details about *obliquity*, check `AFNI's documentation
+    <https://sscc.nimh.nih.gov/sscc/dglen/Obliquity>_.
+
+    Parameters
+    ----------
+    affine : 2D array-like
+        Affine transformation array.  Usually shape (4, 4), but can be any 2D
+        array.
+
+    Returns
+    -------
+    angles : 1D array-like
+        The *obliquity* of each axis with respect to the cardinal axes, in radians.
+
+    """
+    vs = voxel_sizes(affine)
+    best_cosines = np.abs((affine[:-1, :-1] / vs).max(axis=1))
+    return np.arccos(best_cosines)

nibabel/tests/data/affine-LAS-itk.tfm

@@ -0,0 +1 @@
+affine-RAS-itk.tfm

nibabel/tests/data/affine-LAS.afni

@@ -0,0 +1 @@
+affine-RAS.afni

nibabel/tests/data/affine-LAS.fsl

@@ -0,0 +1 @@
+affine-RAS.fsl

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

nibabel/tests/data/affine-LPS.afni

@@ -0,0 +1 @@
+affine-RAS.afni

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

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

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

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

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

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

nibabel/tests/data/someones_anatomy.nii.gz

@@ -0,0 +1 @@
+../../../doc/source/downloads/someones_anatomy.nii.gz

nibabel/tests/test_affines.py

@@ -7,7 +7,7 @@
 
 from ..eulerangles import euler2mat
 from ..affines import (AffineError, apply_affine, append_diag, to_matvec,
-                       from_matvec, dot_reduce, voxel_sizes)
+                       from_matvec, dot_reduce, voxel_sizes, obliquity)
 
 
 from nose.tools import assert_equal, assert_raises
@@ -178,3 +178,15 @@ def test_voxel_sizes():
             rot_affine[:3, :3] = rotation
             full_aff = rot_affine.dot(aff)
             assert_almost_equal(voxel_sizes(full_aff), vox_sizes)
+
+
+def test_obliquity():
+    """Check the calculation of inclination of an affine axes."""
+    from math import pi
+    aligned = np.diag([2.0, 2.0, 2.3, 1.0])
+    aligned[:-1, -1] = [-10, -10, -7]
+    R = from_matvec(euler2mat(x=0.09, y=0.001, z=0.001), [0.0, 0.0, 0.0])
+    oblique = R.dot(aligned)
+    assert_almost_equal(obliquity(aligned), [0.0, 0.0, 0.0])
+    assert_almost_equal(obliquity(oblique) * 180 / pi,
+                        [0.0810285, 5.1569949, 5.1569376])

nibabel/tests/test_transform.py

@@ -0,0 +1,75 @@
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+"""Tests of the transform module."""
+import os
+import numpy as np
+from numpy.testing import assert_array_equal, assert_almost_equal, \
+    assert_array_almost_equal
+import pytest
+
+from ..loadsave import load as loadimg
+from ..nifti1 import Nifti1Image
+from ..eulerangles import euler2mat
+from ..affines import from_matvec
+from ..volumeutils import shape_zoom_affine
+from ..transform import linear as nbl
+from ..testing import (assert_equal, assert_not_equal, assert_true,
+                       assert_false, assert_raises, data_path,
+                       suppress_warnings, assert_dt_equal)
+from ..tmpdirs import InTemporaryDirectory
+
+
+SOMEONES_ANATOMY = os.path.join(data_path, 'someones_anatomy.nii.gz')
+# SOMEONES_ANATOMY = os.path.join(data_path, 'someones_anatomy.nii.gz')
+
+
+@pytest.mark.parametrize('image_orientation', ['RAS', 'LAS', 'LPS', 'oblique'])
+def test_affines_save(image_orientation):
+    """Check implementation of exporting affines to formats."""
+    # Generate test transform
+    img = loadimg(SOMEONES_ANATOMY)
+    imgaff = img.affine
+
+    if image_orientation == 'LAS':
+        newaff = imgaff.copy()
+        newaff[0, 0] *= -1.0
+        newaff[0, 3] = imgaff.dot(np.hstack((np.array(img.shape[:3]) - 1, 1.0)))[0]
+        img = Nifti1Image(np.flip(img.get_fdata(), 0), newaff, img.header)
+    elif image_orientation == 'LPS':
+        newaff = imgaff.copy()
+        newaff[0, 0] *= -1.0
+        newaff[1, 1] *= -1.0
+        newaff[:2, 3] = imgaff.dot(np.hstack((np.array(img.shape[:3]) - 1, 1.0)))[:2]
+        img = Nifti1Image(np.flip(np.flip(img.get_fdata(), 0), 1), newaff, img.header)
+    elif image_orientation == 'oblique':
+        A = shape_zoom_affine(img.shape, img.header.get_zooms(), x_flip=False)
+        R = from_matvec(euler2mat(x=0.09, y=0.001, z=0.001))
+        newaff = R.dot(A)
+        img = Nifti1Image(img.get_fdata(), newaff, img.header)
+        img.header.set_qform(newaff, 1)
+        img.header.set_sform(newaff, 1)
+
+    T = from_matvec(euler2mat(x=0.9, y=0.001, z=0.001), [4.0, 2.0, -1.0])
+
+    xfm = nbl.Affine(T)
+    xfm.reference = img
+
+    itk = nbl.load(os.path.join(data_path, 'affine-%s-itk.tfm' % image_orientation),
+                   fmt='itk')
+    fsl = np.loadtxt(os.path.join(data_path, 'affine-%s.fsl' % image_orientation))
+    afni = np.loadtxt(os.path.join(data_path, 'affine-%s.afni' % image_orientation))
+
+    with InTemporaryDirectory():
+        xfm.to_filename('M.tfm', fmt='itk')
+        xfm.to_filename('M.fsl', fmt='fsl')
+        xfm.to_filename('M.afni', fmt='afni')
+
+        nb_itk = nbl.load('M.tfm', fmt='itk')
+        nb_fsl = np.loadtxt('M.fsl')
+        nb_afni = np.loadtxt('M.afni')
+
+    assert_equal(itk, nb_itk)
+    assert_almost_equal(fsl, nb_fsl)
+    assert_almost_equal(afni, nb_afni)
+
+# Create version not aligned to canonical

nibabel/transform/__init__.py

@@ -0,0 +1,23 @@
+# emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
+#
+#   See COPYING file distributed along with the NiBabel package for the
+#   copyright and license terms.
+#
+### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
+"""
+Geometric transforms.
+
+.. currentmodule:: nibabel.transform
+
+.. autosummary::
+   :toctree: ../generated
+
+   transform
+"""
+from .linear import Affine
+from .nonlinear import DeformationFieldTransform
+
+
+__all__ = ['Affine', 'DeformationFieldTransform']