fix: add a first battery of tests

Added some initial tests to check the writing of transforms to ITK and
FSL formats.

Author: oesteban

Diff for: 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
-    <https://github.com/afni/afni/blob/b6a9f7a21c1f3231ff09efbd861f8975ad48e525/src/thd_coords.c#L660-L698>`_
+    <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
     ----------

Diff for: nibabel/tests/data/affine-LAS-itk.tfm

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

Diff for: nibabel/tests/data/affine-LAS.fsl

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

Diff for: 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 for: 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 for: 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 for: 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 for: 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 for: 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 for: nibabel/tests/data/someones_anatomy.nii.gz

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

Diff for: 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 for: 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 for: 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