IT WORKS

PennLINC · Dec 16, 2024 · 6ed525f · 6ed525f
1 parent e2076c0
commit 6ed525f
Show file tree

Hide file tree

Showing 7 changed files with 96 additions and 44 deletions.
diff --git a/qsirecon/data/io_spec.json b/qsirecon/data/io_spec.json
@@ -14,6 +14,10 @@
             "name": "seg",
             "pattern": "(?:^|_)seg-([a-zA-Z0-9]+)"
         },
+        {
+            "name": "dsistudiotemplate",
+            "pattern": "([a-zA-Z0-9_]+)"
+        },
         {
             "name": "model",
             "pattern": "(?:^|_)model-([a-zA-Z0-9]+)"
@@ -40,6 +44,7 @@
         "sub-{subject}[/ses-{session}]/{datatype<dwi>|dwi}/sub-{subject}[_ses-{session}][_acq-{acquisition}][_ce-{ceagent}][_dir-{direction}][_rec-{reconstruction}][_run-{run}][_space-{space}][_cohort-{cohort}][_seg-{seg}][_model-{model}][_bundles-{bundles}][_param-{param}][_bundle-{bundle}][_label-{label}][_desc-{desc}]_{suffix<streamlines>}.{extension<tck|tck.gz|json>|tck.gz}",
         "sub-{subject}[/ses-{session}]/{datatype<dwi>|dwi}/sub-{subject}[_ses-{session}][_acq-{acquisition}][_ce-{ceagent}][_dir-{direction}][_rec-{reconstruction}][_run-{run}][_space-{space}][_cohort-{cohort}][_seg-{seg}][_model-{model}][_bundles-{bundles}][_param-{param}][_bundle-{bundle}][_label-{label}][_desc-{desc}]_{suffix<scalarstats|tdistats|bundlestats|streamlineweights>}.{extension<tsv|csv|json>|tsv}",
         "sub-{subject}[/ses-{session}]/{datatype<dwi>|dwi}/sub-{subject}[_ses-{session}][_acq-{acquisition}][_ce-{ceagent}][_dir-{direction}][_rec-{reconstruction}][_run-{run}][_space-{space}][_cohort-{cohort}][_seg-{seg}][_model-{model}][_bundles-{bundles}][_param-{param}][_bundle-{bundle}][_label-{label}][_desc-{desc}]_{suffix<mapping>}.{extension<map|map.gz|json>|map.gz}",
+        "sub-{subject}[/ses-{session}]/{datatype<dwi>|dwi}/sub-{subject}[_ses-{session}][_acq-{acquisition}][_ce-{ceagent}][_dir-{direction}][_rec-{reconstruction}][_run-{run}][_space-{space}][_cohort-{cohort}][_seg-{seg}][_model-{model}][_bundles-{bundles}][_param-{param}][_bundle-{bundle}][_label-{label}][_desc-{desc}]_{suffix<dwimap>}.fib.gz[.{dsistudiotemplate}].{extension<map|map.gz|json>|map.gz}",
         "sub-{subject}[/ses-{session}]/{datatype<dwi|anat>|dwi}/sub-{subject}[_ses-{session}][_acq-{acquisition}][_ce-{ceagent}][_dir-{direction}][_rec-{reconstruction}][_run-{run}][_space-{space}][_cohort-{cohort}][_seg-{seg}][_model-{model}][_bundles-{bundles}][_param-{param}][_bundle-{bundle}][_label-{label}][_desc-{desc}]_{suffix<mask>}.{extension<nii|nii.gz|json>|nii.gz}",
         "sub-{subject}[/ses-{session}]/{datatype<dwi|anat>|dwi}/sub-{subject}[_ses-{session}][_acq-{acquisition}][_ce-{ceagent}][_dir-{direction}][_rec-{reconstruction}][_run-{run}][_space-{space}][_cohort-{cohort}][_seg-{seg}][_model-{model}][_bundles-{bundles}][_param-{param}][_bundle-{bundle}][_label-{label}][_desc-{desc}]_{suffix<dseg|probseg>}.{extension<mif|mif.gz|nii.gz|nii|txt|json>|nii.gz}",
         "sub-{subject}[/ses-{session}]/{datatype<dwi>|dwi}/sub-{subject}[_ses-{session}][_acq-{acquisition}][_ce-{ceagent}][_dir-{direction}][_rec-{reconstruction}][_run-{run}][_space-{space}][_cohort-{cohort}][_seg-{seg}][_model-{model}][_bundles-{bundles}][_param-{param}][_bundle-{bundle}][_label-{label}][_desc-{desc}]_{suffix<mu>}.{extension<txt|json>|txt}",

diff --git a/qsirecon/data/pipelines/ss3t_fod_autotrack.yaml b/qsirecon/data/pipelines/ss3t_fod_autotrack.yaml
@@ -10,6 +10,20 @@ nodes:
     qsirecon_suffix: SS3TAutoTrack
     software: DSI Studio
 
+    # Get 3D images of DSI Studio's scalar maps
+-   action: export
+    input: dsistudio_gqi
+    name: gqi_scalars
+    qsirecon_suffix: SS3TAutoTrack
+    software: DSI Studio
+
+    # Perform the registration using the GQI-based QA+ISO
+-   action: autotrack_registration
+    input: dsistudio_gqi
+    name: autotrack_gqi_registration
+    # qsirecon_suffix: Don't include here - the map.gz is saved in autotrack
+    software: DSI Studio
+
     # Select just one shell + b=0 for autotrack
 -   action: select_gradients
     input: qsirecon
@@ -33,31 +47,15 @@ nodes:
             algorithm: dhollander
     qsirecon_suffix: SS3TAutoTrack
 
-    # Get 3D images of DSI Studio's scalar maps
--   action: export
-    input: dsistudio_gqi
-    name: gqi_scalars
-    qsirecon_suffix: SS3TAutoTrack
-    software: DSI Studio
-
-# Uncomment to map scalar maps to template space
-# -   action: template_map
-#     input: qsirecon
-#     name: template_map
-#     parameters:
-#         interpolation: NearestNeighbor
-#     scalars_from:
-#     - gqi_scalars
-#     - dipy_dki
-#     software: qsirecon
-
 -   action: fod_fib_merge
     name: create_fod_fib
-    # to include the fib file and the map file
-    input: dsistudio_gqi
+    # to include the fib file AND the map file
+    input: autotrack_gqi_registration
     csd_input: ss3t_csd
     # outputs include the FOD fib file and the map file is passed through
     qsirecon_suffix: SS3TAutoTrack
+    parameters:
+        model: ss3t
 
 -   action: autotrack
     input: create_fod_fib
@@ -67,5 +65,16 @@ nodes:
         track_id: Association,Projection,Commissure,Cerebellum,CranialNerve
         track_voxel_ratio: 2.0
         yield_rate: 1.0e-06
+        model: ss3t
     qsirecon_suffix: SS3TAutoTrack
     software: DSI Studio
+
+# Uncomment to map scalar maps to template space
+# -   action: template_map
+#     input: qsirecon
+#     name: template_map
+#     parameters:
+#         interpolation: NearestNeighbor
+#     scalars_from:
+#     - gqi_scalars
+#     software: qsirecon
diff --git a/qsirecon/interfaces/converters.py b/qsirecon/interfaces/converters.py
@@ -5,6 +5,7 @@
 import os
 import os.path as op
 import re
+import shutil
 import subprocess
 from pathlib import Path
 
@@ -245,7 +246,28 @@ def _run_interface(self, runtime):
             path_fod_fib=self.inputs.csd_fib_file,
             merged_fib=merged_fib_file,
         )
-        self._results["fibgz"] = merged_fib_file
+
+        # gzip the merged file
+        merged_fibgz_file = merged_fib_file + ".gz"
+
+        p = subprocess.Popen(
+            ["gzip", "-v", merged_fib_file], stdout=subprocess.PIPE, stderr=subprocess.PIPE
+        )
+        response = p.communicate()
+        if not p.returncode == 0:
+            raise Exception(f"Gzip exitted with code {p.returncode}: {response}")
+        if not Path(merged_fibgz_file).exists():
+            raise Exception(f"Failed to gzip {merged_fib_file}")
+        self._results["fibgz"] = merged_fibgz_file
+
+        # Handle the map file if it was provided
+        if isdefined(self.inputs.fibgz_map):
+            LOGGER.info(f"Creating new map file to match {merged_fib_file}.")
+            # DSI Studio stores the template of the mapping file like icbm_adult.map.gz
+            dsistudiotemplate = self.inputs.fibgz_map.split(".")[-3]
+            new_mapping_file = merged_fibgz_file + f".{dsistudiotemplate}.map.gz"
+            shutil.copyfile(self.inputs.fibgz_map, new_mapping_file)
+            self._results["fibgz_map"] = new_mapping_file
 
         return runtime
 

diff --git a/qsirecon/interfaces/dsi_studio.py b/qsirecon/interfaces/dsi_studio.py
@@ -31,6 +31,14 @@
 
 LOGGER = logging.getLogger("nipype.interface")
 DSI_STUDIO_VERSION = "94b9c79"
+DSI_STUDIO_TEMPLATES = [
+    "c57bl6_mouse",
+    "icbm152_adult",
+    "indi_rhesus",
+    "pitt_marmoset",
+    "whs_sd_rat",
+    "dhcp_neonate",
+]
 
 
 class DSIStudioCommandLineInputSpec(CommandLineInputSpec):
@@ -712,6 +720,7 @@ class _AutoTrackOutputSpec(TraitedSpec):
     native_trk_files = OutputMultiObject(File(exists=True))
     stat_files = OutputMultiObject(File(exists=True))
     map_file = File(exists=True)
+    dsistudiotemplate = traits.Str(desc="DSI Studio's name for the template used for registration")
 
 
 class AutoTrack(CommandLine):
@@ -738,7 +747,16 @@ def _list_outputs(self):
             raise Exception("Too many map files generated")
         if not map_files:
             raise Exception("No map files found in " + str(cwd.absolute()))
-        outputs["map_file"] = str(map_files[0].absolute())
+        map_path = map_files[0]
+        outputs["map_file"] = str(map_path.absolute())
+
+        # Which of the template spaces was used?
+        template_space = traits.Undefined
+        for _template_name in DSI_STUDIO_TEMPLATES:
+            if _template_name in map_path.name:
+                template_space = _template_name
+        outputs["dsistudiotemplate"] = template_space
+
         return outputs
 
 

diff --git a/qsirecon/workflows/recon/build_workflow.py b/qsirecon/workflows/recon/build_workflow.py
@@ -13,6 +13,7 @@
 )
 from .dsi_studio import (
     init_dsi_studio_autotrack_wf,
+    init_dsi_studio_autotrack_registration_wf,
     init_dsi_studio_connectivity_wf,
     init_dsi_studio_export_wf,
     init_dsi_studio_recon_wf,
@@ -271,6 +272,8 @@ def workflow_from_spec(inputs_dict, node_spec):
             return init_dsi_studio_tractography_wf(**kwargs)
         if node_spec["action"] == "connectivity":
             return init_dsi_studio_connectivity_wf(**kwargs)
+        if node_spec["action"] == "autotrack_registration":
+            return init_dsi_studio_autotrack_registration_wf(**kwargs)
         if node_spec["action"] == "autotrack":
             return init_dsi_studio_autotrack_wf(**kwargs)
 

diff --git a/qsirecon/workflows/recon/converters.py b/qsirecon/workflows/recon/converters.py
@@ -152,14 +152,13 @@ def init_fod_fib_wf(inputs_dict, name="fod_fib", qsirecon_suffix="", params={}):
                 dismiss_entities=("desc",),
                 suffix="dwimap",
                 extension=".fib.gz",
-                model="ss3t",
+                model=params.get("model", "csd"),
                 compress=True,
             ),
             name="ds_fibgz",
             run_without_submitting=True,
         )
-        workflow.connect(outputnode, 'fibgz',
-                         ds_fibgz, 'in_file')  # fmt:skip
+        workflow.connect(outputnode, 'fibgz', ds_fibgz, 'in_file')  # fmt:skip
 
     return clean_datasinks(workflow, qsirecon_suffix)
 

diff --git a/qsirecon/workflows/recon/dsi_studio.py b/qsirecon/workflows/recon/dsi_studio.py
@@ -304,26 +304,20 @@ def init_dsi_studio_autotrack_registration_wf(
     outputnode.inputs.recon_scalars = []
 
     omp_nthreads = config.nipype.omp_nthreads
-    bundle_names = _get_dsi_studio_bundles(params.get("track_id", ""))
-    bundle_desc = (
-        "AutoTrack attempted to reconstruct the following bundles:\n  * "
-        + "\n  * ".join(bundle_names)
-        + "\n\n"
-    )
-    LOGGER.info(bundle_desc)
-
     workflow = Workflow(name=name)
 
-    # Run autotrack on only one bundle
-    actual_trk = pe.Node(
-        AutoTrack(num_threads=omp_nthreads, **params), name="actual_trk", n_procs=omp_nthreads
-    )  # An extra thread is needed
+    # Run autotrack on only one bundle. The important part is getting the map.gz
+    registration_atk = pe.Node(
+        AutoTrack(num_threads=omp_nthreads, track_id="Association_ArcuateFasciculusL"),
+        name="registration_atk",
+        n_procs=omp_nthreads,
+    )
 
     workflow.connect([
-        (inputnode, actual_trk, [
+        (inputnode, registration_atk, [
             ('fibgz', 'fib_file')]),
         (inputnode, outputnode, [('fibgz', 'fibgz')]),
-        (actual_trk, outputnode, [('map_file', 'fibgz_map')]),
+        (registration_atk, outputnode, [('map_file', 'fibgz_map')]),
     ])  # fmt:skip
 
     return clean_datasinks(workflow, qsirecon_suffix)
@@ -396,7 +390,7 @@ def init_dsi_studio_autotrack_wf(
         "DSI Studio (version %s) and bundle shape statistics were calculated [@autotrack]. "
         % DSI_STUDIO_VERSION
     )
-    model_name = params.get("model_name", "gqi")
+    model_name = params.pop("model", "gqi")
     omp_nthreads = config.nipype.omp_nthreads
     bundle_names = _get_dsi_studio_bundles(params.get("track_id", ""))
     bundle_desc = (
@@ -460,9 +454,9 @@ def init_dsi_studio_autotrack_wf(
     ds_mapping = pe.Node(
         DerivativesDataSink(
             dismiss_entities=("desc",),
-            suffix="mapping",
+            suffix="dwimap",
             model=model_name,
-            extension=".fib.gz.map.gz",
+            extension="map.gz",
             compress=True,
         ),
         name="ds_mapping",
@@ -472,10 +466,12 @@ def init_dsi_studio_autotrack_wf(
     workflow.connect([
         (inputnode, actual_trk, [
             ('fibgz', 'fib_file'),
-            ('fibgz_map', 'fibgz_map')]),
+            ('fibgz_map', 'map_file')]),
         (inputnode, aggregate_atk_results, [('dwi_file', 'source_file')]),
         (inputnode, convert_to_tck, [('dwi_file', 'reference_nifti')]),
-        (actual_trk, ds_mapping, [('map_file', 'in_file')]),
+        (actual_trk, ds_mapping, [
+            ('map_file', 'in_file'),
+            ('dsistudiotemplate', 'dsistudiotemplate')]),
         (actual_trk, aggregate_atk_results, [
             ("native_trk_files", "trk_files"),
             ("stat_files", "stat_files"),