Merge pull request #420 from matthew-brett/master

MRG: add Gregory Lee's new PARREC data

Test parrec data with multiple rotations.

Modify translation order to minimize error between our rotations and
those from the NIfTI images in the conversion.

Author: matthew-brett

.gitmodules

@@ -10,3 +10,6 @@
 [submodule "nibabel-data/nitest-freesurfer"]
 	path = nibabel-data/nitest-freesurfer
 	url = https://bitbucket.org/nipy/nitest-freesurfer.git
+[submodule "nibabel-data/parrec_oblique"]
+	path = nibabel-data/parrec_oblique
+	url = https://github.com/grlee77/parrec_oblique.git

nibabel-data/parrec_oblique

@@ -126,6 +126,8 @@
                       [0, 0, 0, 1]])
 )
 
+DEG2RAD = np.pi / 180.
+
 # General information dict definitions
 # assign props to PAR header entries
 # values are: (shortname[, dtype[, shape]])
@@ -892,11 +894,13 @@ def get_affine(self, origin='scanner'):
                 "Unknown slice orientation ({0}).".format(slice_orientation))
         # hdr has deg, we need radians
         # Order is [ap, fh, rl]
-        ang_rad = self.general_info['angulation'] * np.pi / 180.0
-        # euler2mat accepts z, y, x angles and does rotation around z, y, x
-        # axes in that order. It's possible that PAR assumes rotation in a
-        # different order, we still need some relevant data to test this
-        rot = from_matvec(euler2mat(*ang_rad[::-1]), [0, 0, 0])
+        ap_rot, fh_rot, rl_rot = self.general_info['angulation'] * DEG2RAD
+        Mx = euler2mat(x=ap_rot)
+        My = euler2mat(y=fh_rot)
+        Mz = euler2mat(z=rl_rot)
+        # By trial and error, this unexpected order of rotations seem to give
+        # the closest to the observed (converted NIfTI) affine.
+        rot = from_matvec(dot_reduce(Mz, Mx, My))
         # compose the PSL affine
         psl_aff = dot_reduce(rot, permute_to_psl, zoomer, to_center)
         if origin == 'scanner':

nibabel/parrec.py

@@ -9,6 +9,7 @@
 
 from .. import load as top_load
 from ..parrec import load
+from ..affines import voxel_sizes
 
 from .nibabel_data import get_nibabel_data, needs_nibabel_data
 
@@ -19,9 +20,10 @@
 from numpy.testing import assert_almost_equal
 
 BALLS = pjoin(get_nibabel_data(), 'nitest-balls1')
+OBLIQUE = pjoin(get_nibabel_data(), 'parrec_oblique')
 
 # Amount by which affine translation differs from NIFTI conversion
-AFF_OFF = [-0.93644031, -0.95572686, 0.03288748]
+AFF_OFF = [-0.93575081, -0.95657335, 0.03264122]
 
 
 @needs_nibabel_data('nitest-balls1')
@@ -45,8 +47,8 @@ def test_loading():
             aff_off = pimg.affine[:3, 3] - nimg.affine[:3, 3]
             assert_almost_equal(aff_off, AFF_OFF, 4)
             # The difference is max in the order of 0.5 voxel
-            vox_sizes = np.sqrt((nimg.affine[:3, :3] ** 2).sum(axis=0))
-            assert_true(np.all(np.abs(aff_off / vox_sizes) <= 0.5))
+            vox_sizes = voxel_sizes(nimg.affine)
+            assert_true(np.all(np.abs(aff_off / vox_sizes) <= 0.501))
             # The data is very close, unless it's the fieldmap
             if par_root != 'fieldmap':
                 assert_true(np.allclose(pimg.dataobj, nimg.dataobj))
@@ -63,3 +65,25 @@ def test_fieldmap():
     load(fieldmap_par)
     top_load(fieldmap_nii)
     raise SkipTest('Fieldmap remains puzzling')
+
+
+@needs_nibabel_data('parrec_oblique')
+def test_oblique_loading():
+    # Test loading of oblique parrec files
+    for par in glob(pjoin(OBLIQUE, 'PARREC', '*.PAR')):
+        par_root, ext = splitext(basename(par))
+        # Check we can load the image
+        pimg = load(par)
+        assert_equal(pimg.shape, (560, 560, 1))
+        # Compare against NIFTI if present
+        nifti_fname = pjoin(OBLIQUE, 'NIFTI', par_root + '.nii')
+        nimg = top_load(nifti_fname)
+        assert_almost_equal(nimg.affine[:3, :3], pimg.affine[:3, :3], 3)
+        # The translation part is always off
+        # The ammount differs by rotation
+        aff_off = pimg.affine[:3, 3] - nimg.affine[:3, 3]
+        # The difference is max in the order of 0.5 voxel
+        vox_sizes = voxel_sizes(nimg.affine)
+        assert_true(np.all(np.abs(aff_off / vox_sizes) <= 0.5))
+        # The data is very close
+        assert_true(np.allclose(pimg.dataobj, nimg.dataobj))

nibabel/tests/test_parrec_data.py

@@ -263,7 +263,7 @@ def test_parrec2nii_with_data():
                 assert_almost_equal(aff_off, AFF_OFF, 3)
                 # The difference is max in the order of 0.5 voxel
                 vox_sizes = vox_size(philips_img.affine)
-                assert_true(np.all(np.abs(aff_off / vox_sizes) <= 0.5))
+                assert_true(np.all(np.abs(aff_off / vox_sizes) <= 0.501))
                 # The data is very close, unless it's the fieldmap
                 if par_root != 'fieldmap':
                     conved_data_lps = flip_axis(conved_img.dataobj, 1)