Replace MCFLIRT with 3dVolReg

Author: chrisgorgo

fmriprep/utils/misc.py

@@ -7,6 +7,58 @@
 """
 
 
+def split_and_deoblique_func(in_file):
+    import os
+    from nilearn.image import iter_img
+    import nibabel as nb
+    import numpy as np
+    out_files = []
+    for i, img in enumerate(iter_img(in_file)):
+        out_file = os.path.abspath('vol%04d.nii.gz'%i)
+        affine = img.affine
+        affine[:3,:3] = np.diag(np.diag(affine[:3, :3]))
+        nb.Nifti1Image(np.asanyarray(img.dataobj), affine, img.header).to_filename(out_file)
+        out_files.append(out_file)
+    return out_files
+
+
+def afni2itk_func(in_file):
+    import os
+    from scipy.io import loadmat, savemat
+    from numpy import loadtxt, around, hstack, vstack, zeros, float64
+
+    def read_afni_affine(input_file, debug=False):
+        orig_afni_mat = loadtxt(input_file)
+        if debug:
+            print(orig_afni_mat)
+
+        output = []
+        for i in range(orig_afni_mat.shape[0]):
+            output.append(vstack((orig_afni_mat[i,:].reshape(3, 4, order='C'), [0, 0, 0, 1])))
+        return output
+
+    def get_ants_dict(affine, debug=False):
+        out_dict = {}
+        out_dict['AffineTransform_double_3_3'] = hstack(
+            (affine[:3, :3].reshape(1, -1), affine[:3, 3].reshape(1, -1))).reshape(-1, 1).astype(float64)
+        out_dict['fixed'] = zeros((3, 1))
+        if debug:
+            print(out_dict)
+
+        return out_dict
+
+    out_file = os.path.abspath('mc4d.txt')
+    with open(out_file, 'w') as fp:
+        fp.write("#Insight Transform File V1.0\n")
+        for i, affine in enumerate(read_afni_affine(in_file)):
+            fp.write("#Transform %d\n"%i)
+            fp.write("Transform: AffineTransform_double_3_3\n")
+            trans_dict = get_ants_dict(affine)
+            fp.write("Parameters: " + ' '.join(["%g"%i for i in list(trans_dict['AffineTransform_double_3_3'])]) + "\n")
+            fp.write("FixedParameters: " + ' '.join(["%g"%i for i in list(trans_dict['fixed'])]) + "\n")
+    return out_file
+
+
 def fix_multi_T1w_source_name(in_files):
     """
     Make up a generic source name when there are multiple T1s

fmriprep/workflows/bold/base.py

@@ -19,6 +19,7 @@
 from nipype.pipeline import engine as pe
 from nipype.interfaces import utility as niu
 
+from fmriprep.utils.misc import split_and_deoblique_func
 from ...interfaces import (
     DerivativesDataSink,
     GiftiNameSource,
@@ -371,7 +372,8 @@ def init_func_preproc_wf(bold_file, ignore, freesurfer,
         omp_nthreads=omp_nthreads, enhance_t2=True)
 
     # Top-level BOLD splitter
-    bold_split = pe.Node(FSLSplit(dimension='t'), name='bold_split',
+    bold_split = pe.Node(niu.Function(function=split_and_deoblique_func, input_names=['in_file'],
+                                      output_names=['out_files']), name='bold_split',
                          mem_gb=mem_gb['filesize'] * 3)
 
     # HMC on the BOLD

fmriprep/workflows/bold/hmc.py

@@ -10,8 +10,9 @@
 """
 
 from nipype.pipeline import engine as pe
-from nipype.interfaces import utility as niu, fsl
+from nipype.interfaces import utility as niu, fsl, afni
 from niworkflows.interfaces import NormalizeMotionParams
+from fmriprep.utils.misc import afni2itk_func
 from ...engine import Workflow
 from ...interfaces import MCFLIRT2ITK
 
@@ -73,24 +74,24 @@ def init_bold_hmc_wf(mem_gb, omp_nthreads, name='bold_hmc_wf'):
         name='outputnode')
 
     # Head motion correction (hmc)
-    mcflirt = pe.Node(fsl.MCFLIRT(save_mats=True, save_plots=True),
-                      name='mcflirt', mem_gb=mem_gb * 3)
+    mc = pe.Node(afni.Volreg(args='-prefix NULL -twopass',
+                             zpad=4, outputtype='NIFTI_GZ'), name="mc", mem_gb=mem_gb * 3)
 
-    fsl2itk = pe.Node(MCFLIRT2ITK(), name='fsl2itk',
-                      mem_gb=0.05, n_procs=omp_nthreads)
+    afni2itk = pe.Node(niu.Function(function=afni2itk_func,
+                                    input_names=["in_file"],
+                                    output_names=["out_files"]), name='afni2itk',
+                       mem_gb=0.05)
 
-    normalize_motion = pe.Node(NormalizeMotionParams(format='FSL'),
+    normalize_motion = pe.Node(NormalizeMotionParams(format='AFNI'),
                                name="normalize_motion",
                                mem_gb=DEFAULT_MEMORY_MIN_GB)
 
     workflow.connect([
-        (inputnode, mcflirt, [('raw_ref_image', 'ref_file'),
-                              ('bold_file', 'in_file')]),
-        (inputnode, fsl2itk, [('raw_ref_image', 'in_source'),
-                              ('raw_ref_image', 'in_reference')]),
-        (mcflirt, fsl2itk, [('mat_file', 'in_files')]),
-        (mcflirt, normalize_motion, [('par_file', 'in_file')]),
-        (fsl2itk, outputnode, [('out_file', 'xforms')]),
+        (inputnode, mc, [('raw_ref_image', 'basefile'),
+                         ('bold_file', 'in_file')]),
+        (mc, afni2itk, [('oned_matrix_save', 'in_file')]),
+        (mc, normalize_motion, [('oned_file', 'in_file')]),
+        (afni2itk, outputnode, [('out_files', 'xforms')]),
         (normalize_motion, outputnode, [('out_file', 'movpar_file')]),
     ])
 

fmriprep/workflows/bold/util.py

@@ -102,16 +102,18 @@ def init_bold_reference_wf(omp_nthreads, bold_file=None, name='bold_reference_wf
 
     gen_ref = pe.Node(EstimateReferenceImage(), name="gen_ref",
                       mem_gb=1)  # OE: 128x128x128x50 * 64 / 8 ~ 900MB.
+    deoblique = pe.Node(afni.Refit(deoblique=True), name='deoblique')
     enhance_and_skullstrip_bold_wf = init_enhance_and_skullstrip_bold_wf(omp_nthreads=omp_nthreads)
 
     workflow.connect([
         (inputnode, validate, [('bold_file', 'in_file')]),
         (validate, gen_ref, [('out_file', 'in_file')]),
-        (gen_ref, enhance_and_skullstrip_bold_wf, [('ref_image', 'inputnode.in_file')]),
+        (gen_ref, deoblique, [('ref_image', 'in_file')]),
+        (deoblique, enhance_and_skullstrip_bold_wf, [('out_file', 'inputnode.in_file')]),
         (validate, outputnode, [('out_file', 'bold_file'),
                                 ('out_report', 'validation_report')]),
-        (gen_ref, outputnode, [('ref_image', 'raw_ref_image'),
-                               ('n_volumes_to_discard', 'skip_vols')]),
+        (deoblique, outputnode, [('out_file', 'raw_ref_image')]),
+        (gen_ref, outputnode, [('n_volumes_to_discard', 'skip_vols')]),
         (enhance_and_skullstrip_bold_wf, outputnode, [
             ('outputnode.bias_corrected_file', 'ref_image'),
             ('outputnode.mask_file', 'bold_mask'),