[FIX] Make hsvs workflow work

.circleci/FreeSurferIngressTest.sh

@@ -1,44 +1,21 @@
 #!/bin/bash
 
-cat << DOC
-
-Test the FreeSurfer Ingress workflow
-====================================
-
-Test the case where QSIPrep has processed the anatomical data and an external
-(likely run as part of fmriprep) freesurfer run has happened as well
-
-This tests the following features:
- - Blip-up + Blip-down DWI series for TOPUP/Eddy
- - Eddy is run on a CPU
- - Denoising is skipped
- - A follow-up reconstruction using the dsi_studio_gqi workflow
-
-Inputs:
--------
-
- - freesurfer SUBJECTS_DIR (data/freesurfer)
- - qsiprep multi shell results with anatomical outputs (data/qsiprep_with_anat)
-
-DOC
-
 set +e
 source ./get_data.sh
 TESTDIR=${PWD}
-#get_config_data ${TESTDIR}
+get_config_data ${TESTDIR}
 #get_bids_data ${TESTDIR} freesurfer
 #get_bids_data ${TESTDIR} abcd_output
 
 CFG=${TESTDIR}/data/nipype.cfg
-EDDY_CFG=${TESTDIR}/data/eddy_config.json
 export FS_LICENSE=${TESTDIR}/data/license.txt
 
 # Test dipy_mapmri
 TESTNAME=fs_ingress_test
-# setup_dir ${TESTDIR}/${TESTNAME}
+setup_dir ${TESTDIR}/${TESTNAME}
 TEMPDIR=${TESTDIR}/${TESTNAME}/work
 OUTPUT_DIR=${TESTDIR}/${TESTNAME}/derivatives
-BIDS_INPUT_DIR=${TESTDIR}/data/qsiprep_with_anat
+BIDS_INPUT_DIR=${TESTDIR}/data/araikes/qsiprep
 SUBJECTS_DIR=${TESTDIR}/data/freesurfer
 QSIPREP_CMD=$(run_qsiprep_cmd ${BIDS_INPUT_DIR} ${OUTPUT_DIR})
 
@@ -49,7 +26,7 @@ ${QSIPREP_CMD} \
 	 --recon-spec ${PWD}/test_5tt_hsv.json \
 	 --freesurfer-input ${SUBJECTS_DIR} \
 	 --recon-only \
-	 -vv
+	 -vv --debug pdb
 
 
 

.circleci/get_data.sh

@@ -484,9 +484,9 @@ setup_dir(){
     DIR=$1
     mkdir -p ${DIR}/derivatives
     mkdir -p ${DIR}/work
-    # setfacl -d -m group:$(id -gn):rwx ${DIR}/derivatives && \
-    #     setfacl -m group:$(id -gn):rwx ${DIR}/derivatives
-    # setfacl -d -m group:$(id -gn):rwx ${DIR}/work && \
-    #     setfacl -m group:$(id -gn):rwx ${DIR}/work
+    setfacl -d -m group:$(id -gn):rwx ${DIR}/derivatives && \
+        setfacl -m group:$(id -gn):rwx ${DIR}/derivatives
+    setfacl -d -m group:$(id -gn):rwx ${DIR}/work && \
+        setfacl -m group:$(id -gn):rwx ${DIR}/work
 
 }

.circleci/test_5tt_hsv.json

@@ -0,0 +1,102 @@
+{
+  "name": "mrtrix_multishell_msmt_hsvs",
+  "space": "T1w",
+  "atlases": ["schaefer100"],
+  "anatomical": [
+    "mrtrix_5tt_hsvs"
+  ],
+  "nodes": [
+    {
+      "name": "msmt_csd",
+      "software": "MRTrix3",
+      "action": "csd",
+      "qsirecon_suffix": "MRtrix3_act-HSVS",
+      "input": "qsiprep",
+      "parameters": {
+        "mtnormalize": true,
+        "response": {
+          "algorithm": "dhollander"
+        },
+        "fod": {
+          "algorithm": "msmt_csd",
+          "max_sh": [
+            8,
+            8,
+            8
+          ]
+        }
+      }
+    },
+    {
+      "name": "track_ifod2",
+      "software": "MRTrix3",
+      "action": "tractography",
+      "qsirecon_suffix": "MRtrix3_act-HSVS",
+      "input": "msmt_csd",
+      "parameters": {
+        "use_5tt": true,
+        "method_5tt": "hsvs",
+        "use_sift2": true,
+        "tckgen": {
+          "algorithm": "iFOD2",
+          "select": 10000000,
+          "max_length": 250,
+          "min_length": 30,
+          "power": 0.33,
+          "crop_at_gmwmi": true,
+          "backtrack": true,
+          "quiet": true
+        },
+        "sift2": {}
+      }
+    },
+    {
+      "name": "mrtrix_conn",
+      "software": "MRTrix3",
+      "action": "connectivity",
+      "qsirecon_suffix": "MRtrix3_act-HSVS",
+      "input": "track_ifod2",
+      "parameters": {
+        "tck2connectome": [
+          {
+            "zero_diagonal": false,
+            "search_radius": 2,
+            "scale_invnodevol": true,
+            "symmetric": true,
+            "use_sift_weights": true,
+            "stat_edge": "sum",
+            "measure": "sift_invnodevol_radius2_count"
+          },
+          {
+            "zero_diagonal": false,
+            "search_radius": 2,
+            "scale_invnodevol": false,
+            "symmetric": true,
+            "length_scale": "length",
+            "use_sift_weights": false,
+            "stat_edge": "mean",
+            "measure": "radius2_meanlength"
+          },
+          {
+            "zero_diagonal": false,
+            "search_radius": 2,
+            "scale_invnodevol": false,
+            "symmetric": true,
+            "use_sift_weights": false,
+            "stat_edge": "sum",
+            "measure": "radius2_count"
+          },
+          {
+            "zero_diagonal": false,
+            "search_radius": 2,
+            "scale_invnodevol": false,
+            "symmetric": true,
+            "use_sift_weights": true,
+            "stat_edge": "sum",
+            "measure": "sift_radius2_count"
+          }
+        ]
+      }
+    }
+  ]
+}

qsiprep/config.py

@@ -437,6 +437,8 @@ class execution(_Config):
     """Root of QSIRecon BIDS Derivatives dataset."""
     recon_input = None
     """Directory containing QSIPrep derivatives to run through recon workflows."""
+    freesurfer_input = None
+    """Directory containing FreeSurfer directories to use for recon workflows."""
     recon_only = False
     """Run only recon workflows."""
     reportlets_dir = None

qsiprep/interfaces/bids.py

@@ -382,14 +382,22 @@ def get_recon_output_name(
     datatype_dir = output_bids_entities.get("datatype", op.basename(op.dirname(source_file)))
     out_path += f"/{datatype_dir}"
     _, source_fname, _ = split_filename(source_file)
+    # Remove the suffix
     source_fname, _ = source_fname.rsplit("_", 1)
     _, _, extension = split_filename(derivative_file)
 
     # It may be that the space has changed. Check if it has
     if "space" in output_bids_entities:
-        source_fname = re.sub(
-            "_space-[a-zA-Z0-9]+", "_space-" + output_bids_entities["space"], source_fname
-        )
+        if "_space-" in source_fname:
+            source_fname = re.sub(
+                "_space-[a-zA-Z0-9]+", "_space-" + output_bids_entities["space"], source_fname
+            )
+        elif "_desc-preproc" in source_fname:
+            source_fname = source_fname.replace(
+                "_desc-preproc", f"_space-{output_bids_entities['space']}_desc-preproc"
+            )
+        else:
+            source_fname += f"_space-{output_bids_entities['space']}"
     base_fname = op.join(out_path, source_fname)
 
     # Add the new bids entities for the output file
@@ -444,6 +452,8 @@ def _run_interface(self, runtime):
         output_bids = {}
         if self.inputs.atlas:
             output_bids["atlas"] = self.inputs.atlas
+        if self.inputs.space:
+            output_bids["space"] = self.inputs.space
         if self.inputs.bundles:
             output_bids["bundles"] = self.inputs.bundles
         if self.inputs.bundle:

qsiprep/interfaces/freesurfer.py

@@ -20,6 +20,7 @@
 
 """
 import os.path as op
+from pathlib import Path
 
 import nibabel as nb
 import numpy as np
@@ -454,10 +455,10 @@ def find_fs_path(freesurfer_dir, subject_id):
         return None
     nosub = op.join(freesurfer_dir, subject_id)
     if op.exists(nosub):
-        return nosub
+        return Path(nosub)
     withsub = op.join(freesurfer_dir, "sub-" + subject_id)
     if op.exists(withsub):
-        return withsub
+        return Path(withsub)
     return None
 
 

qsiprep/interfaces/interchange.py

@@ -9,7 +9,12 @@
 from qsiprep.interfaces.ingress import QsiReconDWIIngress
 
 # Anatomical (t1w/t2w) slots
-FS_FILES_TO_REGISTER = ["brain", "aseg"]
+FS_FILES_TO_REGISTER = [
+    # The skull-stripped brain produced by FreeSurfer
+    "brain",
+    # Automated volumetric segmentation file from FreeSurfer
+    "aseg",
+]
 CREATEABLE_ANATOMICAL_OUTPUTS = [
     "fs_5tt_hsvs",
     "qsiprep_5tt_hsvs",

qsiprep/workflows/recon/anatomical.py

@@ -41,6 +41,7 @@
 HSV_REQUIREMENTS = [
     "mri/aparc+aseg.mgz",
     "mri/brainmask.mgz",
+    "mri/norm.mgz",
     "mri/transforms/talairach.xfm",
     "surf/lh.white",
     "surf/lh.pial",
@@ -79,13 +80,13 @@ def init_highres_recon_anatomical_wf(
     outputnode = pe.Node(
         niu.IdentityInterface(fields=anatomical_workflow_outputs), name="outputnode"
     )
-
+    workflow.__desc__ = ""
     # "Gather" the input data. ``status`` is a dict that reflects which anatomical data
     # are present. The anat_ingress_node is a nipype node that ensures that qsiprep-style
     # anatomical data is available. In the case where ``pipeline_source`` is not "qsiprep",
     # the data is converted in this node to be qsiprep-like.
     pipeline_source = config.workflow.recon_input_pipeline
-    freesurfer_dir = config.execution.fs_subjects_dir
+    freesurfer_dir = config.execution.freesurfer_input
     if pipeline_source == "qsiprep":
         anat_ingress_node, status = gather_qsiprep_anatomical_data(subject_id)
     elif pipeline_source == "ukb":
@@ -126,12 +127,15 @@ def init_highres_recon_anatomical_wf(
             "Freesurfer directory %s exists for %s", subject_freesurfer_path, subject_id
         )
         fs_source = pe.Node(
-            nio.FreeSurferSource(subjects_dir=freesurfer_dir, subject_id="sub-" + subject_id),
+            nio.FreeSurferSource(
+                subjects_dir=str(subject_freesurfer_path.parent),
+                subject_id=subject_freesurfer_path.name,
+            ),
             name="fs_source",
         )
-        workflow.connect([
-            (fs_source, outputnode, [
-                (field, field) for field in FS_FILES_TO_REGISTER])])  # fmt:skip
+        # workflow.connect([
+        #     (fs_source, outputnode, [
+        #         (field, field) for field in FS_FILES_TO_REGISTER])])  # fmt:skip
 
     # Do we need to calculate anything else on the
     if "mrtrix_5tt_hsvs" in extras_to_make:
@@ -143,28 +147,39 @@ def init_highres_recon_anatomical_wf(
         config.loggers.workflow.info("FreeSurfer data will be used to create a HSVS 5tt image.")
         status["has_freesurfer_5tt_hsvs"] = True
         create_5tt_hsvs = pe.Node(
-            GenerateMasked5tt(algorithm="hsvs", in_file=str(subject_freesurfer_path)),
+            GenerateMasked5tt(
+                algorithm="hsvs",
+                in_file=str(subject_freesurfer_path),
+                nthreads=config.nipype.omp_nthreads,
+            ),
             name="create_5tt_hsvs",
+            n_procs=config.nipype.omp_nthreads,
         )
-        workflow.connect([
-            (create_5tt_hsvs, outputnode, [('out_file', 'fs_5tt_hsvs')])])  # fmt:skip
         ds_qsiprep_5tt_hsvs = pe.Node(
-            ReconDerivativesDataSink(atlas="hsvs", suffix="dseg", qsirecon_suffix="anat"),
+            ReconDerivativesDataSink(
+                atlas="hsvs", space="T1w", suffix="dseg", qsirecon_suffix="anat"
+            ),
             name="ds_qsiprep_5tt_hsvs",
             run_without_submitting=True,
         )
         ds_fs_5tt_hsvs = pe.Node(
             ReconDerivativesDataSink(
-                desc="hsvs", space="fsnative", suffix="dseg", qsirecon_suffix="anat"
+                desc="hsvs", space="fsnative", suffix="dseg", qsirecon_suffix="anat", compress=True
             ),
             name="ds_fs_5tt_hsvs",
             run_without_submitting=True,
         )
+        workflow.connect([
+            (anat_ingress_node, ds_fs_5tt_hsvs, [("t1_preproc", "source_file")]),
+            (anat_ingress_node, ds_qsiprep_5tt_hsvs, [("t1_preproc", "source_file")]),
+            (create_5tt_hsvs, outputnode, [('out_file', 'fs_5tt_hsvs')]),
+            (create_5tt_hsvs, ds_fs_5tt_hsvs, [("out_file", "in_file")]),
+        ])  # fmt:skip
 
         # Transform the 5tt image so it's registered to the QSIPrep AC-PC T1w
         if status["has_qsiprep_t1w"]:
             config.loggers.workflow.info(
-                "HSVS 5tt imaged will be registered to the " "QSIPrep T1w image."
+                "HSVS 5tt imaged will be registered to the QSIPrep T1w image."
             )
             status["has_qsiprep_5tt_hsvs"] = True
             register_fs_to_qsiprep_wf = init_register_fs_to_qsiprep_wf(
@@ -184,8 +199,6 @@ def init_highres_recon_anatomical_wf(
                         "fs_to_qsiprep_transform_itk")] + [
                     ("outputnode." + field, field) for field in FS_FILES_TO_REGISTER]),
                 (create_5tt_hsvs, apply_header_to_5tt, [("out_file", "in_image")]),
-                (create_5tt_hsvs, ds_fs_5tt_hsvs, [("out_file", "in_file")]),
-
                 (register_fs_to_qsiprep_wf, apply_header_to_5tt, [
                     ("outputnode.fs_to_qsiprep_transform_mrtrix", "transform_file")]),
                 (apply_header_to_5tt, outputnode, [

qsiprep/workflows/recon/base.py

@@ -256,8 +256,7 @@ def init_single_subject_recon_wf(subject_id):
                 else config.execution.output_dir
             )
             workflow.get_node(_node).inputs.base_directory = base_dir
-            # workflow.get_node(_node).inputs.source_file = \
-            #     "anat/sub-{}_desc-preproc_T1w.nii.gz".format(subject_id)
+
     return workflow