Merge pull request #286 from matthew-brett/refactor-parrec

matthew-brett · matthew-brett · commit 04e6465d43a5 · 2014-12-01T13:27:33.000-06:00
MRG: refactoring parrec2nii while adding tests

Refactoring parrec2nii to make it possible to save nifti images in LAS voxel orientation. 

* add tests for minmax,  scaling options, writing PAR in nifti extension;
* fix writing PAR into nifti extension;
* add reorientation to LAS;
* test reorientation to LAS in voxel and bvecs.
diff --git a/bin/parrec2nii b/bin/parrec2nii
@@ -5,16 +5,19 @@ from __future__ import division, print_function, absolute_import
 
 from optparse import OptionParser, Option
 import numpy as np
+import numpy.linalg as npl
 import sys
 import os
-import gzip
 import nibabel
 import nibabel.parrec as pr
 from nibabel.parrec import one_line
 from nibabel.mriutils import calculate_dwell_time, MRIError
 import nibabel.nifti1 as nifti1
 from nibabel.filename_parser import splitext_addext
-from nibabel.volumeutils import array_to_file, fname_ext_ul_case
+from nibabel.volumeutils import fname_ext_ul_case
+from nibabel.orientations import (io_orientation, inv_ornt_aff,
+                                  apply_orientation)
+from nibabel.affines import apply_affine, from_matvec, to_matvec
 
 
 def get_opt_parser():
@@ -136,82 +139,84 @@ def proc_file(infile, opts):
                       'overwrite it' % outfilename)
 
     # load the PAR header and data
-    scaling = None if opts.scaling == 'off' else opts.scaling
+    scaling = 'dv' if opts.scaling == 'off' else opts.scaling
     infile = fname_ext_ul_case(infile)
     pr_img = pr.load(infile,
                      permit_truncated=opts.permit_truncated,
                      scaling=scaling)
     pr_hdr = pr_img.header
-    raw_data = pr_img.dataobj.get_unscaled()
     affine = pr_hdr.get_affine(origin=opts.origin)
-    nimg = nifti1.Nifti1Image(raw_data, affine, pr_hdr)
+    slope, intercept = pr_hdr.get_data_scaling(scaling)
+    if opts.scaling != 'off':
+        verbose('Using data scaling "%s"' % opts.scaling)
+    # get original scaling, and decide if we scale in-place or not
+    if opts.scaling == 'off':
+        slope = np.array([1.])
+        intercept = np.array([0.])
+        in_data = pr_img.dataobj.get_unscaled()
+        out_dtype = pr_hdr.get_data_dtype()
+    elif not np.any(np.diff(slope)) and not np.any(np.diff(intercept)):
+        # Single scalefactor case
+        slope = slope.ravel()[0]
+        intercept = intercept.ravel()[0]
+        in_data = pr_img.dataobj.get_unscaled()
+        out_dtype = pr_hdr.get_data_dtype()
+    else:
+        # Multi scalefactor case
+        slope = np.array([1.])
+        intercept = np.array([0.])
+        in_data = np.array(pr_img.dataobj)
+        out_dtype = np.float64
+    # Reorient data block to LAS+ if necessary
+    ornt = io_orientation(np.diag([-1, 1, 1, 1]).dot(affine))
+    if np.all(ornt == [[0, 1],
+                       [1, 1],
+                       [2, 1]]): # already in LAS+
+        t_aff = np.eye(4)
+    else: # Not in LAS+
+        t_aff = inv_ornt_aff(ornt, pr_img.shape)
+        affine = np.dot(affine, t_aff)
+        in_data = apply_orientation(in_data, ornt)
+    # Make corresponding NIfTI image
+    nimg = nifti1.Nifti1Image(in_data, affine, pr_hdr)
     nhdr = nimg.header
+    nhdr.set_data_dtype(out_dtype)
+    nhdr.set_slope_inter(slope, intercept)
 
     if 'parse' in opts.minmax:
         # need to get the scaled data
         verbose('Loading (and scaling) the data to determine value range')
-        if opts.scaling == 'off':
-            scaled_data = raw_data
-        else:
-            slope, intercept = pr_hdr.get_data_scaling(method=opts.scaling)
-            scaled_data = slope * raw_data + intercept
     if opts.minmax[0] == 'parse':
-        nhdr.structarr['cal_min'] = scaled_data.min()
+        nhdr['cal_min'] = in_data.min() * slope + intercept
     else:
-        nhdr.structarr['cal_min'] = float(opts.minmax[0])
+        nhdr['cal_min'] = float(opts.minmax[0])
     if opts.minmax[1] == 'parse':
-        nhdr.structarr['cal_max'] = scaled_data.max()
+        nhdr['cal_max'] = in_data.max() * slope + intercept
     else:
-        nhdr.structarr['cal_max'] = float(opts.minmax[1])
+        nhdr['cal_max'] = float(opts.minmax[1])
 
     # container for potential NIfTI1 header extensions
     if opts.store_header:
-        exts = nifti1.Nifti1Extensions()
         # dump the full PAR header content into an extension
         with open(infile, 'r') as fobj:
             hdr_dump = fobj.read()
             dump_ext = nifti1.Nifti1Extension('comment', hdr_dump)
-        exts.append(dump_ext)
-    # finalize the header: set proper data offset, pixdims, ...
-    nimg.update_header()
-
-    # prep a file
-    if opts.compressed:
-        verbose('Using gzip compression')
-        outfile = gzip.open(outfilename, 'wb')
-    else:
-        outfile = open(outfilename, 'wb')
+        nhdr.extensions.append(dump_ext)
 
-    # get original scaling, and decide if we scale in-place or not
-    if opts.scaling == 'off':
-        slope = np.array([1.])
-        intercept = np.array([0.])
-    else:
-        verbose('Using data scaling "%s"' % opts.scaling)
-        slope, intercept = pr_hdr.get_data_scaling(method=opts.scaling)
     verbose('Writing %s' % outfilename)
-    if not np.any(np.diff(slope)) and not np.any(np.diff(intercept)):
-        # Single scalefactor case
-        nhdr.set_slope_inter(slope.ravel()[0], intercept.ravel()[0])
-        data_obj = raw_data
-    else:
-        # Multi scalefactor case
-        nhdr.set_slope_inter(1, 0)
-        nhdr.set_data_dtype(np.float64)
-        data_obj = pr_img.dataobj
-    nhdr.write_to(outfile)
-    offset = nhdr.get_data_offset()
-    array_to_file(data_obj,
-                  outfile,
-                  out_dtype=nhdr.get_data_dtype(), # for endianness
-                  offset=offset)
+    nibabel.save(nimg, outfilename)
+
     # write out bvals/bvecs if requested
     if opts.bvs:
         bvals, bvecs = pr_hdr.get_bvals_bvecs()
         if bvals is None and bvecs is None:
             verbose('No DTI volumes detected, bvals and bvecs not written')
         else:
             verbose('Writing .bvals and .bvecs files')
+            # Transform bvecs with reorientation affine
+            orig2new = npl.inv(t_aff)
+            bv_reorient = from_matvec(to_matvec(orig2new)[0], [0, 0, 0])
+            bvecs = apply_affine(bv_reorient, bvecs)
             with open(basefilename + '.bvals', 'w') as fid:
                 # np.savetxt could do this, but it's just a loop anyway
                 for val in bvals:
diff --git a/nibabel/parrec.py b/nibabel/parrec.py
@@ -929,7 +929,7 @@ def get_data_scaling(self, method="dv"):
             slope = 1.0 / scale_slope
             intercept = rescale_intercept / (rescale_slope * scale_slope)
         else:
-            raise ValueError("Unknown scling method '%s'." % method)
+            raise ValueError("Unknown scaling method '%s'." % method)
         reorder = self.get_sorted_slice_indices()
         slope = slope[reorder]
         intercept = intercept[reorder]
diff --git a/nibabel/tests/test_scripts.py b/nibabel/tests/test_scripts.py
@@ -16,16 +16,20 @@
 
 from ..tmpdirs import InTemporaryDirectory
 from ..loadsave import load
+from ..orientations import flip_axis, aff2axcodes, inv_ornt_aff
 
 from nose.tools import (assert_true, assert_false, assert_not_equal,
                         assert_equal)
 
-from numpy.testing import assert_almost_equal
+from numpy.testing import assert_almost_equal, assert_array_equal
 
 from .scriptrunner import ScriptRunner
 from .nibabel_data import needs_nibabel_data
-from .test_parrec import DTI_PAR_BVECS, DTI_PAR_BVALS
+from ..testing import assert_dt_equal
+from .test_parrec import (DTI_PAR_BVECS, DTI_PAR_BVALS,
+                          EXAMPLE_IMAGES as PARREC_EXAMPLES)
 from .test_parrec_data import BALLS, AFF_OFF
+from .test_helpers import assert_data_similar
 
 
 def _proc_stdout(stdout):
@@ -83,31 +87,88 @@ def vox_size(affine):
     return np.sqrt(np.sum(affine[:3,:3] ** 2, axis=0))
 
 
+def check_conversion(cmd, pr_data, out_fname):
+    run_command(cmd)
+    img = load(out_fname)
+    # Check orientations always LAS
+    assert_equal(aff2axcodes(img.affine), tuple('LAS'))
+    data = img.get_data()
+    assert_true(np.allclose(data, pr_data))
+    assert_true(np.allclose(img.header['cal_min'], data.min()))
+    assert_true(np.allclose(img.header['cal_max'], data.max()))
+    # Check minmax options
+    run_command(cmd + ['--minmax', '1', '2'])
+    img = load(out_fname)
+    assert_true(np.allclose(data, pr_data))
+    assert_true(np.allclose(img.header['cal_min'], 1))
+    assert_true(np.allclose(img.header['cal_max'], 2))
+    run_command(cmd + ['--minmax', 'parse', '2'])
+    img = load(out_fname)
+    assert_true(np.allclose(data, pr_data))
+    assert_true(np.allclose(img.header['cal_min'], data.min()))
+    assert_true(np.allclose(img.header['cal_max'], 2))
+    run_command(cmd + ['--minmax', '1', 'parse'])
+    img = load(out_fname)
+    assert_true(np.allclose(data, pr_data))
+    assert_true(np.allclose(img.header['cal_min'], 1))
+    assert_true(np.allclose(img.header['cal_max'], data.max()))
+
+
 @script_test
 def test_parrec2nii():
     # Test parrec2nii script
     cmd = ['parrec2nii', '--help']
     code, stdout, stderr = run_command(cmd)
     assert_true(stdout.startswith('Usage'))
-    in_fname = pjoin(DATA_PATH, 'phantom_EPI_asc_CLEAR_2_1.PAR')
-    out_froot = 'phantom_EPI_asc_CLEAR_2_1.nii'
     with InTemporaryDirectory():
-        run_command(['parrec2nii', in_fname])
-        img = load(out_froot)
-        assert_equal(img.shape, (64, 64, 9, 3))
-        assert_equal(img.get_data_dtype(), np.dtype(np.uint16))
-        # Check against values from Philips converted nifti image
-        data = img.get_data()
-        assert_true(np.allclose(
-            (data.min(), data.max(), data.mean()),
-            (0.0, 2299.4110643863678, 194.95876256117265)))
-        assert_almost_equal(vox_size(img.get_affine()), (3.75, 3.75, 8))
+        for eg_dict in PARREC_EXAMPLES:
+            fname = eg_dict['fname']
+            run_command(['parrec2nii', fname])
+            out_froot = splitext(basename(fname))[0] + '.nii'
+            img = load(out_froot)
+            assert_equal(img.shape, eg_dict['shape'])
+            assert_dt_equal(img.get_data_dtype(), eg_dict['dtype'])
+            # Check against values from Philips converted nifti image
+            data = img.get_data()
+            assert_data_similar(data, eg_dict)
+            assert_almost_equal(img.header.get_zooms(), eg_dict['zooms'])
+            # Standard save does not save extensions
+            assert_equal(len(img.header.extensions), 0)
+            # Does not overwrite unless option given
+            code, stdout, stderr = run_command(
+                ['parrec2nii', fname], check_code=False)
+            assert_equal(code, 1)
+            # Default scaling is dv
+            pr_img = load(fname)
+            flipped_data = flip_axis(pr_img.get_data(), 1)
+            base_cmd = ['parrec2nii', '--overwrite', fname]
+            check_conversion(base_cmd, flipped_data, out_froot)
+            check_conversion(base_cmd + ['--scaling=dv'],
+                             flipped_data,
+                             out_froot)
+            # fp
+            pr_img = load(fname, scaling='fp')
+            flipped_data = flip_axis(pr_img.get_data(), 1)
+            check_conversion(base_cmd + ['--scaling=fp'],
+                             flipped_data,
+                             out_froot)
+            # no scaling
+            unscaled_flipped = flip_axis(pr_img.dataobj.get_unscaled(), 1)
+            check_conversion(base_cmd + ['--scaling=off'],
+                             unscaled_flipped,
+                             out_froot)
+            # Save extensions
+            run_command(base_cmd + ['--store-header'])
+            img = load(out_froot)
+            assert_equal(len(img.header.extensions), 1)
 
 
 @script_test
 @needs_nibabel_data('nitest-balls1')
 def test_parrec2nii_with_data():
     # Use nibabel-data to test conversion
+    # Premultiplier to relate our affines to Philips conversion
+    LAS2LPS = inv_ornt_aff([[0, 1], [1, -1], [2, 1]], (80, 80, 10))
     with InTemporaryDirectory():
         for par in glob(pjoin(BALLS, 'PARREC', '*.PAR')):
             par_root, ext = splitext(basename(par))
@@ -118,23 +179,31 @@ def test_parrec2nii_with_data():
             # Do conversion
             run_command(['parrec2nii', par])
             conved_img = load(par_root + '.nii')
+            # Confirm parrec2nii conversions are LAS
+            assert_equal(aff2axcodes(conved_img.affine), tuple('LAS'))
+            # Shape same whether LPS or LAS
             assert_equal(conved_img.shape[:3], (80, 80, 10))
-            # Test against converted NIfTI
+            # Test against original converted NIfTI
             nifti_fname = pjoin(BALLS, 'NIFTI', par_root + '.nii.gz')
             if exists(nifti_fname):
-                nimg = load(nifti_fname)
-                assert_almost_equal(nimg.affine[:3, :3],
-                                    conved_img.affine[:3, :3], 3)
+                philips_img = load(nifti_fname)
+                # Confirm Philips converted image always LPS
+                assert_equal(aff2axcodes(philips_img.affine), tuple('LPS'))
+                # Equivalent to Philips LPS affine
+                equiv_affine = conved_img.affine.dot(LAS2LPS)
+                assert_almost_equal(philips_img.affine[:3, :3],
+                                    equiv_affine[:3, :3], 3)
                 # The translation part is always off by the same ammout
-                aff_off = conved_img.affine[:3, 3] - nimg.affine[:3, 3]
-                assert_almost_equal(aff_off, AFF_OFF, 4)
+                aff_off = equiv_affine[:3, 3] - philips_img.affine[:3, 3]
+                assert_almost_equal(aff_off, AFF_OFF, 3)
                 # The difference is max in the order of 0.5 voxel
-                vox_sizes = np.sqrt((nimg.affine[:3, :3] ** 2).sum(axis=0))
+                vox_sizes = vox_size(philips_img.affine)
                 assert_true(np.all(np.abs(aff_off / vox_sizes) <= 0.5))
                 # The data is very close, unless it's the fieldmap
                 if par_root != 'fieldmap':
-                    assert_true(np.allclose(conved_img.dataobj,
-                                            nimg.dataobj))
+                    conved_data_lps = flip_axis(conved_img.dataobj, 1)
+                    assert_true(np.allclose(conved_data_lps,
+                                            philips_img.dataobj))
     with InTemporaryDirectory():
         # Test some options
         dti_par = pjoin(BALLS, 'PARREC', 'DTI.PAR')
@@ -149,8 +218,12 @@ def test_parrec2nii_with_data():
         # Writes bvals, bvecs files if asked
         run_command(['parrec2nii', '--overwrite', '--bvs', dti_par])
         assert_almost_equal(np.loadtxt('DTI.bvals'), DTI_PAR_BVALS)
-        assert_almost_equal(np.loadtxt('DTI.bvecs'),
-                            DTI_PAR_BVECS[:, [2, 0, 1]].T)
+        # Bvecs in header, transposed from PSL to LPS
+        bvecs_LPS = DTI_PAR_BVECS[:, [2, 0, 1]]
+        # Adjust for output flip of Y axis in data and bvecs
+        bvecs_LAS = bvecs_LPS * [1, -1, 1]
+        assert_almost_equal(np.loadtxt('DTI.bvecs'), bvecs_LAS.T)
+        # Dwell time
         assert_false(exists('DTI.dwell_time'))
         # Need field strength if requesting dwell time
         code, _, _, = run_command(
