ENH: Replace MCFLIRT with 3dVolReg

fmriprep_rodents/interfaces/mc.py

@@ -0,0 +1,53 @@
+from pathlib import Path
+
+from nipype.utils.filemanip import fname_presuffix
+from nipype.interfaces.base import (
+    BaseInterfaceInputSpec,
+    TraitedSpec,
+    SimpleInterface,
+    File,
+)
+import numpy as np
+
+
+class Volreg2ITKInputSpec(BaseInterfaceInputSpec):
+    in_file = File(exists=True, mandatory=True,
+                   desc='mat file generated by AFNI\'s 3dVolreg')
+
+
+class Volreg2ITKOutputSpec(TraitedSpec):
+    out_file = File(desc='the output ITKTransform file')
+
+
+class Volreg2ITK(SimpleInterface):
+
+    """
+    Convert an AFNI's mat file into an ITK Transform file.
+    """
+    input_spec = Volreg2ITKInputSpec
+    output_spec = Volreg2ITKOutputSpec
+
+    def _run_interface(self, runtime):
+        # Load AFNI mat entries and reshape appropriately
+        orig_afni_mat = np.loadtxt(self.inputs.in_file)
+        afni_affines = [mat.reshape(3, 4, order='C') for mat in orig_afni_mat]
+
+        out_file = Path(fname_presuffix(self.inputs.in_file, use_ext=False,
+                        suffix='_mc4d_itk.txt', newpath=runtime.cwd))
+
+        fixed_params = 'FixedParameters: 0 0 0'  # Center of rotation does not change
+        lines = ["#Insight Transform File V1.0"]
+        for i, affine in enumerate(afni_affines):
+            lines.append("#Transform %d" % i)
+            lines.append("Transform: AffineTransform_double_3_3")
+
+            ants_affine_2d = np.hstack((affine[:3, :3].reshape(1, -1),
+                                        affine[:3, 3].reshape(1, -1)))
+            params = ants_affine_2d.reshape(-1).astype('float64')
+            params_list = ["%g" % i for i in params.tolist()]
+            lines.append("Parameters: %s" % ' '.join(params_list))
+            lines.append(fixed_params)
+
+        out_file.write_text('\n'.join(lines))
+        self._results['out_file'] = str(out_file)
+        return runtime

fmriprep_rodents/workflows/bold/hmc.py

@@ -9,7 +9,8 @@
 """
 
 from nipype.pipeline import engine as pe
-from nipype.interfaces import utility as niu, fsl
+from nipype.interfaces import utility as niu, afni
+from ...interfaces.mc import Volreg2ITK
 
 from ...config import DEFAULT_MEMORY_MIN_GB
 
@@ -67,8 +68,8 @@ def init_bold_hmc_wf(mem_gb, omp_nthreads, name='bold_hmc_wf'):
 Head-motion parameters with respect to the BOLD reference
 (transformation matrices, and six corresponding rotation and translation
 parameters) are estimated before any spatiotemporal filtering using
-`mcflirt` [FSL {fsl_ver}, @mcflirt].
-""".format(fsl_ver=fsl.Info().version() or '<ver>')
+`3dVolReg` [AFNI {afni_ver}, @afni, RRID:SCR_005927].
+""".format(afni_ver=afni.Info().version() or '<ver>')
 
     inputnode = pe.Node(
         niu.IdentityInterface(fields=['bold_file', 'raw_ref_image']),
@@ -79,29 +80,51 @@ 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, save_rms=True),
-        name='mcflirt', mem_gb=mem_gb * 3)
-
-    fsl2itk = pe.Node(MCFLIRT2ITK(), name='fsl2itk',
-                      mem_gb=0.05, n_procs=omp_nthreads)
-
-    normalize_motion = pe.Node(NormalizeMotionParams(format='FSL'),
+    # mcflirt = pe.Node(
+    #     fsl.MCFLIRT(save_mats=True, save_plots=True, save_rms=True),
+    #     name='mcflirt', mem_gb=mem_gb * 3)
+
+    # fsl2itk = pe.Node(MCFLIRT2ITK(), name='fsl2itk',
+    #                 mem_gb=0.05, n_procs=omp_nthreads)
+
+    mc = pe.Node(
+        afni.Volreg(
+            zpad=4,
+            outputtype='NIFTI_GZ',
+            args='-prefix NULL -twopass'
+        ),
+        name='mc',
+        mem_gb=mem_gb * 3,
+    )
+
+    mc2itk = pe.Node(Volreg2ITK(), name='mcitk', mem_gb=0.05)
+
+    normalize_motion = pe.Node(NormalizeMotionParams(format='AFNI'),
                                name="normalize_motion",
                                mem_gb=DEFAULT_MEMORY_MIN_GB)
 
-    def _pick_rel(rms_files):
-        return rms_files[-1]
-
+    # def _pick_rel(rms_files):
+    #     return rms_files[-1]
+
+    # 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')]),
+    #     (mcflirt, outputnode, [(('rms_files', _pick_rel), 'rmsd_file')]),
+    #     (fsl2itk, outputnode, [('out_file', 'xforms')]),
+    #     (normalize_motion, outputnode, [('out_file', 'movpar_file')]),
+    # ])
     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')]),
-        (mcflirt, outputnode, [(('rms_files', _pick_rel), 'rmsd_file')]),
-        (fsl2itk, outputnode, [('out_file', 'xforms')]),
+        (inputnode, mc, [
+            ('raw_ref_image', 'basefile'),
+            ('bold_file', 'in_file'),
+        ]),
+        (mc, mc2itk, [('oned_matrix_save', 'in_file')]),
+        (mc, normalize_motion, [('oned_file', 'in_file')]),
+        (mc2itk, outputnode, [('out_file', 'xforms')]),
         (normalize_motion, outputnode, [('out_file', 'movpar_file')]),
     ])