[ENH] Replace FSL mcflirt with AFNI 3dVolReg

fmriprep/interfaces/images.py

@@ -293,6 +293,7 @@ def _run_interface(self, runtime):
 
 class ValidateImageInputSpec(BaseInterfaceInputSpec):
     in_file = File(exists=True, mandatory=True, desc='input image')
+    deoblique = traits.Bool(False, usedefault=True)
 
 
 class ValidateImageOutputSpec(TraitedSpec):
@@ -342,6 +343,15 @@ class ValidateImage(SimpleInterface):
     output_spec = ValidateImageOutputSpec
 
     def _run_interface(self, runtime):
+
+        def deoblique(img):
+            import os
+            import nibabel as nb
+            import numpy as np
+            affine = img.affine
+            affine[:3, :3] = np.diag(np.diag(affine[:3, :3]))
+            return nb.Nifti1Image(np.asanyarray(img.dataobj), affine, img.header)
+
         img = nb.load(self.inputs.in_file)
         out_report = os.path.join(runtime.cwd, 'report.html')
 
@@ -370,7 +380,13 @@ def _run_interface(self, runtime):
 
         # Both match, qform valid (implicit with match), codes okay -> do nothing, empty report
         if matching_affines and qform_code > 0 and sform_code > 0:
-            self._results['out_file'] = self.inputs.in_file
+            if self.inputs.deoblique:
+                out_fname = fname_presuffix(self.inputs.in_file, suffix='_valid', newpath=runtime.cwd)
+                img = deoblique(img)
+                img.to_filename(out_fname)
+                self._results['out_file'] = out_fname
+            else:
+                self._results['out_file'] = self.inputs.in_file
             open(out_report, 'w').close()
             self._results['out_report'] = out_report
             return runtime
@@ -418,6 +434,9 @@ def _run_interface(self, runtime):
 </p>
 """
         snippet = '<h3 class="elem-title">%s</h3>\n%s\n' % (warning_txt, description)
+
+        if self.inputs.deoblique:
+            img = deoblique(img)
         # Store new file and report
         img.to_filename(out_fname)
         with open(out_report, 'w') as fobj:

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,
@@ -396,7 +397,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

@@ -104,7 +104,7 @@ 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.
     # Re-run validation; no effect if no sbref; otherwise apply same validation to sbref as bold
-    validate_ref = pe.Node(ValidateImage(), name='validate_ref', mem_gb=DEFAULT_MEMORY_MIN_GB)
+    validate_ref = pe.Node(ValidateImage(deoblique=True), name='validate_ref', mem_gb=DEFAULT_MEMORY_MIN_GB)
     enhance_and_skullstrip_bold_wf = init_enhance_and_skullstrip_bold_wf(omp_nthreads=omp_nthreads)
 
     workflow.connect([