FIX: Remove trace by default

Author: larsoner

bin/parrec2nii

@@ -103,6 +103,11 @@ def get_opt_parser():
                    store the corresponding settings in the NIfTI header, unless
                    non-uniform scaling is used, in which case the data is
                    stored in the file in scaled form. Default: 'dv'""")))
+    p.add_option(
+        Option('--keep-trace', action="store_true", dest='keep_trace',
+               default=False,
+               help=one_line("""Do not discard the diagnostic Philips DTI
+                             trace volume, if it exists in the data.""")))
     p.add_option(
         Option("--overwrite", action="store_true", dest="overwrite",
                default=False,
@@ -171,12 +176,26 @@ def proc_file(infile, opts):
     ornt = io_orientation(np.diag([-1, 1, 1, 1]).dot(affine))
     if np.all(ornt == [[0, 1],
                        [1, 1],
-                       [2, 1]]): # already in LAS+
+                       [2, 1]]):  # already in LAS+
         t_aff = np.eye(4)
-    else: # Not in LAS+
+    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)
+
+    bvals, bvecs = pr_hdr.get_bvals_bvecs()
+    if not opts.keep_trace:  # discard Philips DTI trace if present
+        if bvals is not None:
+            bad_mask = np.logical_and(bvals != 0, (bvecs == 0).all(axis=1))
+            if bad_mask.sum() > 0:
+                pl = 's' if bad_mask.sum() != 1 else ''
+                verbose('Removing %s DTI trace volume%s'
+                        % (bad_mask.sum(), pl))
+                good_mask = ~bad_mask
+                in_data = in_data[..., good_mask]
+                bvals = bvals[good_mask]
+                bvecs = bvecs[good_mask]
+
     # Make corresponding NIfTI image
     nimg = nifti1.Nifti1Image(in_data, affine, pr_hdr)
     nhdr = nimg.header
@@ -208,7 +227,6 @@ def proc_file(infile, opts):
 
     # 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:

nibabel/tests/test_scripts.py

@@ -221,12 +221,17 @@ def test_parrec2nii_with_data():
         assert_true(exists('DTI.nii'))
         assert_false(exists('DTI.bvals'))
         assert_false(exists('DTI.bvecs'))
+        # ensure trace was removed
+        img = load('DTI.nii')
+        assert_equal(img.shape[-1], len(DTI_PAR_BVALS) - 1)
+        del img
         # Does not overwrite unless option given
         code, stdout, stderr = run_command(['parrec2nii', dti_par],
                                            check_code=False)
         assert_equal(code, 1)
         # Writes bvals, bvecs files if asked
-        run_command(['parrec2nii', '--overwrite', '--bvs', dti_par])
+        run_command(['parrec2nii', '--overwrite', '--keep-trace',
+                     '--bvs', dti_par])
         assert_almost_equal(np.loadtxt('DTI.bvals'), DTI_PAR_BVALS)
         # Bvecs in header, transposed from PSL to LPS
         bvecs_LPS = DTI_PAR_BVECS[:, [2, 0, 1]]