MRG: add Gregory Lee's new PARREC data

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

@@ -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]])
@@ -876,11 +878,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/tests/test_parrec_data.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_scripts.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)