Merge pull request #82 from mwaskom/regspace_mat

Save func2anat registration matrix to reg directory for epi registrations

doc/release/v0.0.10.txt

@@ -0,0 +1,16 @@
+
+v0.0.10 (Unreleased)
+--------------------
+
+Command line interface
+~~~~~~~~~~~~~~~~~~~~~~
+
+- The full command-line argument namespace is now saved out with the other
+  experiment parameters for better reproducibility.
+
+Registration workflow
+~~~~~~~~~~~~~~~~~~~~~
+
+- The correct registration matrix to go from the epi space to the anatomy is
+  now written out during the registration so that downstream results from
+  experiments that use the ``-regexp`` flag are correct.

lyman/frontend.py

@@ -33,7 +33,7 @@ def gather_project_info():
     return project_dict
 
 
-def gather_experiment_info(exp_name=None, altmodel=None):
+def gather_experiment_info(exp_name=None, altmodel=None, args=None):
     """Import an experiment module and add some formatted information."""
     lyman_dir = os.environ["LYMAN_DIR"]
 
@@ -89,6 +89,10 @@ def keep(k):
     # Build contrasts list if neccesary
     exp_dict["contrast_names"] = [c[0] for c in exp_dict["contrasts"]]
 
+    # Add command line arguments for reproducibility
+    if args is not None:
+        exp_dict["command_line"] = vars(args)
+
     return exp_dict
 
 

lyman/tools/fileutils.py

@@ -22,6 +22,7 @@ class SaveParameters(BaseInterface):
 
     input_spec = SaveParametersInput
     output_spec = SaveParametersOutput
+    _always_run = True
 
     def _run_interface(self, runtime):
 

lyman/workflows/registration.py

@@ -15,6 +15,7 @@
 """
 import os
 import os.path as op
+import shutil
 import numpy as np
 
 from nipype import Workflow, Node, IdentityInterface
@@ -40,12 +41,16 @@ def create_reg_workflow(name="reg", space="mni",
         fields = ["copes", "varcopes", "sumsquares"]
     elif regtype == "timeseries":
         fields = ["timeseries"]
+
     if cross_exp:
         fields.extend(["first_rigid"])
+
     fields.extend(["means", "masks", "rigids"])
 
     if space == "mni":
         fields.extend(["affine", "warpfield"])
+    else:
+        fields.extend(["tkreg_rigid"])
 
     inputnode = Node(IdentityInterface(fields), "inputnode")
 
@@ -112,6 +117,7 @@ class EPIRegInput(BaseInterfaceInputSpec):
 
     rigids = InputMultiPath(File(exists=True))
     first_rigid = File(exists=True)
+    tkreg_rigid = File(exists=True)
 
 
 class MNIModelRegInput(MNIRegInput,
@@ -368,6 +374,12 @@ def _run_interface(self, runtime):
                 runtime = self.apply_fsl_rigid(runtime, in_file,
                                                out_file, full_rigid)
 
+            # Copy the matrix to go from this space to the anatomy
+            if not i:
+                tkreg_fname = op.basename(self.inputs.tkreg_rigid)
+                out_tkreg = op.join(out_dir, tkreg_fname)
+                shutil.copyfile(self.inputs.tkreg_rigid, out_tkreg)
+
         self.out_files = out_files
         return runtime
 
@@ -471,5 +483,11 @@ def _run_interface(self, runtime):
             runtime = self.apply_fsl_rigid(runtime, run_mean,
                                            out_mean, full_rigid)
 
+            # Copy the matrix to go from this space to the anatomy
+            if not i:
+                tkreg_fname = op.basename(self.inputs.tkreg_rigid)
+                out_tkreg = op.join(out_dir, tkreg_fname)
+                shutil.copyfile(self.inputs.tkreg_rigid, out_tkreg)
+
         self.out_files = out_files
         return runtime

lyman/workflows/surfols.py

@@ -132,7 +132,7 @@ class RemoveEmpty(BaseInterface):
     def _run_interface(self, runtime):
 
         good_images = [f for f in self.inputs.in_files
-                       if not np.allclose(nib.load(f).get_data(), 0).any()]
+                       if not np.allclose(nib.load(f).get_data(), 0)]
         self.good_images = good_images
         return runtime
 

scripts/run_fmri.py

@@ -32,7 +32,7 @@ def main(arglist):
 
     # Get and process specific information
     project = lyman.gather_project_info()
-    exp = lyman.gather_experiment_info(args.experiment, args.altmodel)
+    exp = lyman.gather_experiment_info(args.experiment, args.altmodel, args)
 
     # Set up the SUBJECTS_DIR for Freesurfer
     os.environ["SUBJECTS_DIR"] = project["data_dir"]
@@ -224,6 +224,14 @@ def main(arglist):
                                              "normalization/warpfield.nii.gz")
         reg_templates["affine"] = op.join(data_dir, "{subject_id}",
                                           "normalization/affine." + aff_ext)
+    else:
+        if args.regexp is None:
+            tkreg_base = analysis_dir
+        else:
+            tkreg_base = op.join(project["analysis_dir"], args.regexp)
+        reg_templates["tkreg_rigid"] = op.join(tkreg_base,
+                                               "{subject_id}", "preproc",
+                                               "run_1", "func2anat_tkreg.dat")
 
     # Rigid (6dof) functional-to-anatomical matrices
     rigid_stem = "{subject_id}/preproc/run_*/func2anat_"
@@ -286,7 +294,8 @@ def main(arglist):
     wf_name = space + "_ffx"
     ffx, ffx_input, ffx_output = wf.create_ffx_workflow(wf_name,
                                                         space,
-                                                        exp["contrast_names"])
+                                                        exp["contrast_names"],
+                                                        exp_info=exp)
 
     ext = "_warp.nii.gz" if space == "mni" else "_xfm.nii.gz"
     ffx_base = op.join("{subject_id}/reg", space, "{smoothing}/run_*")
@@ -307,8 +316,9 @@ def main(arglist):
         reg = op.join(os.environ["FREESURFER_HOME"],
                       "average/mni152.register.dat")
     else:
-        bg = "{subject_id}/preproc/run_1/mean_func.nii.gz"
-        reg = "{subject_id}/preproc/run_1/func2anat_tkreg.dat"
+        reg_dir = "{subject_id}/reg/epi/{smoothing}/run_1"
+        bg = op.join(reg_dir, "mean_func_xfm.nii.gz")
+        reg = op.join(reg_dir, "func2anat_tkreg.dat")
     ffx_templates["anatomy"] = bg
     ffx_templates["reg_file"] = reg
 

scripts/run_group.py

@@ -28,7 +28,7 @@ def main(arglist):
 
     # Get and process specific information
     project = lyman.gather_project_info()
-    exp = lyman.gather_experiment_info(args.experiment, args.altmodel)
+    exp = lyman.gather_experiment_info(args.experiment, args.altmodel, args)
 
     if args.experiment is None:
         args.experiment = project["default_exp"]
@@ -83,8 +83,9 @@ def main(arglist):
             dofs=op.join(mfx_base, "tdof_t1.nii.gz")))
     else:
         templates.update(dict(
-            reg_file=op.join(anal_dir_base, "{subject_id}/preproc/run_1",
-                             "func2anat_tkreg.dat")))
+            reg_file=op.join(anal_dir_base,
+                             "{subject_id}/reg/epi/", ffxsmooth,
+                             "run_1/func2anat_tkreg.dat")))
 
     # Workflow source node
     mfx_source = MapNode(SelectFiles(templates,