proposal for slice timing correction

[chrplr]

Currently there is a 'slice_order' variable which specifies the slice order and is applied to all functional files.

To handle more complex cases (e.g. multiband acquisition), we can pass individual slice timings (in millseconds), as well as a reference time (in millisec). SPM12 slice timing correction handles this (see below).

Maybe (?), in cases where there are different slice timings for different functional files, we could have:

session_1_slice_order = ... session_2_slice_order = ...

in the same way that, currently, we write:

session_1_func = session_2_func =

matlabbatch var to control slice timing correction in spm12: %%% slice matlabbatch{1}.spm.temporal.st.scans = {''}; matlabbatch{1}.spm.temporal.st.nslices = 10; matlabbatch{1}.spm.temporal.st.tr = 2; matlabbatch{1}.spm.temporal.st.ta = 1.8; matlabbatch{1}.spm.temporal.st.so = [100 50 30 10 20 100 50 30 10 20]; % either order or timings in ms matlabbatch{1}.spm.temporal.st.refslice = 200; % either order or a time point in ms matlabbatch{1}.spm.temporal.st.prefix = 'a';

[bthirion]

This sounds good. My comment is that in 95% of the cases, the slice order is fixed across sessions, hence could be fixed one and for all. So, the session_XXX_slice_order specification should not be required, but only taken into account when provided.

Also, i would first enforce SPM12 convention on slice timing specification before adding session-specific slice order.

[chrplr]

Here is a short python to extract timing information from dicom files (thanks to Martin Perez)

#! /usr/bin/env python

from __future__ import print_function
import os
import glob
import dicom
import sys

try:
    dicom_path = sys.argv[1]

    if os.path.isdir(dicom_path):
        dicom_ref = sorted(glob.glob(os.path.join(dicom_path, '*.dcm')))[4]
    else:
        if (os.path.isfile(dicom_path)):
            dicom_ref = dicom_path

    TR = dicom.read_file(dicom_ref).RepetitionTime
    slice_times = dicom.read_file(dicom_ref)[0x19, 0x1029].value
    nb_slices = len(slice_times)

except:
    print("Unexpected error: %s" % sys.exc_info()[0])
    print("\nUsage:\n    %s dicom_path\nWhere dicom_path is a dicom directory" % sys.argv[0])
    sys.exit(-1)

print("TR = %.3f" % (TR/1000.))
#print("nb_slices = %d" % nb_slices)
print("slice_timings = ", end=" ")
for v in slice_times:
    print("%.1f" % v, end=" ")
print("\n", end=" ")
print("refslice = %.3f" % (TR/2000.))

[chrplr]

For info: I have put an example dataset (just one EPI session, both in dicom and nifti format) on the www.nitrc.org website. It is available at the link:

    https://www.nitrc.org/frs/download.php/9029/epi-example-dicom.tar.gz

Note: This tarball also contains an ini file with the slice times. This can be used as a testbed for slice timing correction. Todo: integrate as an example of an *.ini in pypreprocess, with the dataset fetcher.

[mrahim]

In SPM12, Nipype can take slice timepoints in milliseconds instead of slice order. This does not change anything in the API. We can set it in slice_order of the config.ini file.

[isadenghien]

Hi, I've tried to use slice_order with slice in time. I used the slice_order into the config.ini like this:

# Can be ascending, descending, or an explicitly specified sequence.
# slice_order = ascending

slice_order = 0, 60, 120, 180, 240, 300, 362, 422, 482, 542, 602, 662, 722, 782, 842, 902, 962, 1025, 1085, 1145, 1205, 1265, 1325, 1385, 1445, 1505, 1565, 1627, 1687, 1747, 1807, 1867, 1927, 1987, 2047, 2107, 2167, 2227, 2290, 2350

# Where the EPI slices interleaved ?
interleaved = False

But it failed with this kind of message:

  File "/home/id999999/.local/lib/python2.7/site-packages/pypreprocess/slice_timing.py", line 69, in get_slice_indices
    slice_order < n_slices)), slice_order
AssertionError: [  -1   59  119  179  239  299  361  421  481  541  601  661  721  781  841
  901  961 1024 1084 1144 1204 1264 1324 1384 1444 1504 1564 1626 1686 1746
 1806 1866 1926 1986 2046 2106 2166 2226 2289 2349]

I suppose that pypreprocess uses these values like number of slice and not in time. I updated all packages like traits, nipy and nipype. Perhaps something is wrong in my config.ini ? Or the slice_order is not taken into account with values in time ?

[mrahim]

There is an assert to check whether slice_order values are less than the total number of slices. So we need to consider this case if we want time slice instead of slice order.

[chrplr]

I could advance a somewhat extremist proposal, that is, to force the user enter slice timings in msec in all cases.

My experience is that keywords like ascending, intererleaved (coming from SPM usage) can be confusing for some users, and I have indeed had students making mistakes. Also, about hte use of sequences of number to specify the order, I am not sure about what happens if we reuse the same number for MB sequences.

IMHO, it is much cleaner to enforce the use of timings.

Alas, this would only work out of the box for spm12 but not for spm8, so you probably do not want to accept this proposal. :-)

[chrplr]

To further test slice timing algorithms, I have created an image file with sinewave timeseries with different phases along z.

See nifti/sinewave.nii in https://www.nitrc.org/frs/download.php/9035/epi-example-dicom.tar.gz

[bthirion]

You're probably right, but we cannot break the compatibility with SPM8...

[chrplr]

A quick hack of nipype_preproc_spm_utils.py to run spm12 slice timing with 'slice_order' variable set as a list of times of acquisition in millieconds. Seems to run, but needs to be tested on a synthetic data set. Need to cleaning and check that is does not break spm8.

cp983411@is151622:~/code/pypreprocess/pypreprocess$ git diff nipype_preproc_spm_utils.py
diff --git a/pypreprocess/nipype_preproc_spm_utils.py b/pypreprocess/nipype_preproc_spm_utils.py
index a5f6df8..dfc2896 100644
--- a/pypreprocess/nipype_preproc_spm_utils.py
+++ b/pypreprocess/nipype_preproc_spm_utils.py
@@ -143,20 +143,24 @@ def _do_subject_slice_timing(subject_data, TR, TA=None, spm_dir=None,

     # compute nslices
     nslices = load_vols(subject_data.func[0])[0].shape[2]
-    assert 1 <= ref_slice <= nslices, ref_slice
+#    assert 1 <= ref_slice <= nslices, ref_slice

-    # compute slice indices / order
-    if not isinstance(slice_order, basestring):
-        slice_order = np.array(slice_order) - 1
-    slice_order = get_slice_indices(nslices, slice_order=slice_order,
-                                    interleaved=interleaved)
+#    # compute slice indices / order
+#    if not isinstance(slice_order, basestring):
+#        slice_order = np.array(slice_order) - 1
+#    slice_order = get_slice_indices(nslices, slice_order=slice_order,
+#                                    inter  leaved=interleaved)

     # use pure python (pp) code ?
+#    import pdb
+#    pdb.set_trace()
+    print " !!!!!!! DISPLAY !!!!!!!"
+    print ref_slice
+    print slice_order
     if software == "python":
         return _pp_do_subject_slice_timing(subject_data, ref_slice=ref_slice,
                                            slice_order=slice_order,
-                                           caching=caching)
-
+                                           caching=caching)                                 
     # sanitize software choice
     if software != "spm":
         raise NotImplementedError(
@@ -193,8 +197,9 @@ def _do_subject_slice_timing(subject_data, TR, TA=None, spm_dir=None,
     for sess_func in subject_data.func:
         stc_result = stc(in_files=sess_func, time_repetition=TR,
                          time_acquisition=TA, num_slices=nslices,
-                         ref_slice=ref_slice + 1,
-                         slice_order=list(slice_order + 1),  # SPM
+                         ref_slice=ref_slice,
+                         #slice_order=list(slice_order + 1),  # SPM8 ?
+                         slice_order=[int(x) for x in slice_order.split()],  # SPM12
                          ignore_exception=False
                          )
         if stc_result.outputs is None:

[bthirion]

This is almost a PR ;-) The only daunting part is to check/enforce SPM8 compliance

[chrplr]

Life gets a bit simpler if we dont enforce SPM8 compliance. To handle the input as either slice order, or slice times, One idea is that two variables could coexist: slice_order and slice_times