From 58c748e34c805e25ef4b417293a8f0cc5be107a0 Mon Sep 17 00:00:00 2001 From: Matthew Brett Date: Thu, 10 Mar 2016 17:38:32 -0800 Subject: [PATCH] TST: 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. --- .gitmodules | 3 +++ nibabel-data/parrec_oblique | 1 + nibabel/parrec.py | 14 +++++++++----- nibabel/tests/test_parrec_data.py | 30 +++++++++++++++++++++++++++--- nibabel/tests/test_scripts.py | 2 +- 5 files changed, 41 insertions(+), 9 deletions(-) create mode 160000 nibabel-data/parrec_oblique diff --git a/.gitmodules b/.gitmodules index f11370ad49..9e03a837ce 100644 --- a/.gitmodules +++ b/.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 diff --git a/nibabel-data/parrec_oblique b/nibabel-data/parrec_oblique new file mode 160000 index 0000000000..45d4d44e17 --- /dev/null +++ b/nibabel-data/parrec_oblique @@ -0,0 +1 @@ +Subproject commit 45d4d44e1783a814cc50990c2a0cca2fd38b245c diff --git a/nibabel/parrec.py b/nibabel/parrec.py index 5e85c411f5..1c2ed034c8 100644 --- a/nibabel/parrec.py +++ b/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': diff --git a/nibabel/tests/test_parrec_data.py b/nibabel/tests/test_parrec_data.py index e428f5e974..630e66cab8 100644 --- a/nibabel/tests/test_parrec_data.py +++ b/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)) diff --git a/nibabel/tests/test_scripts.py b/nibabel/tests/test_scripts.py index 33b0bceee5..fd73e9d809 100644 --- a/nibabel/tests/test_scripts.py +++ b/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)