enh: add tests for affines stored in AFNI format (non-oblique images)

Author: oesteban

nibabel/tests/data/affine-LAS.afni

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

nibabel/tests/data/affine-LPS.afni

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

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/test_transform.py

@@ -57,15 +57,19 @@ def test_affines_save(image_orientation):
     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/linear.py

@@ -204,7 +204,7 @@ def _to_hdf5(self, x5_root):
 
     def to_filename(self, filename, fmt='X5', moving=None):
         """Store the transform in BIDS-Transforms HDF5 file format (.x5)."""
-        if fmt.lower() in ['itk', 'ants', 'elastix', 'nifty']:
+        if fmt.lower() in ['itk', 'ants', 'elastix']:
             with open(filename, 'w') as f:
                 f.write('#Insight Transform File V1.0\n')
 
@@ -229,26 +229,25 @@ def to_filename(self, filename, fmt='X5', moving=None):
             T = self.matrix.copy()
             pre = LPS
             post = LPS
-            if any(obliquity(self.reference.affine) * 180 / pi > OBLIQUITY_THRESHOLD_DEG):
+            if obliquity(self.reference.affine).min() * 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))
+                                      voxel_sizes(M), x_flip=False, y_flip=False)
+                pre = M.dot(np.linalg.inv(A)).dot(LPS)
 
                 if not moving:
                     moving = self.reference
 
-            if moving and any(obliquity(moving.affine) * 180 / pi > OBLIQUITY_THRESHOLD_DEG):
+            if moving and obliquity(moving.affine).min() * 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))
+                M2 = moving.affine
+                A2 = shape_zoom_affine(moving.shape,
+                                       voxel_sizes(M2), x_flip=True, y_flip=True)
+                post = A2.dot(np.linalg.inv(M2))
 
             # 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 = np.swapaxes(post.dot(self.matrix.copy().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')

nibabel/volumeutils.py

@@ -1456,7 +1456,7 @@ def finite_range(arr, check_nan=False):
     return np.nanmin(mins), np.nanmax(maxes)
 
 
-def shape_zoom_affine(shape, zooms, x_flip=True):
+def shape_zoom_affine(shape, zooms, x_flip=True, y_flip=False):
     ''' Get affine implied by given shape and zooms
 
     We get the translations from the center of the image (implied by
@@ -1471,6 +1471,9 @@ def shape_zoom_affine(shape, zooms, x_flip=True):
     x_flip : {True, False}
        whether to flip the X row of the affine.  Corresponds to
        radiological storage on disk.
+    y_flip : {False, True}
+       whether to flip the Y row of the affine.  Corresponds to
+       DICOM storage on disk when x_flip is also True.
 
     Returns
     -------
@@ -1510,6 +1513,9 @@ def shape_zoom_affine(shape, zooms, x_flip=True):
         zooms = full_zooms
     if x_flip:
         zooms[0] *= -1
+
+    if y_flip:
+        zooms[1] *= -1
     # Get translations from center of image
     origin = (shape - 1) / 2.0
     aff = np.eye(4)