From 67cebd9314eebc141729e535cec3e1e35f84d114 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Wed, 29 Mar 2023 14:22:30 -0400
Subject: [PATCH 01/27] Make --basil and --scorescrub mutually exclusive.

---
 aslprep/cli/parser.py | 157 +++++++++++++++++++++++++++---------------
 1 file changed, 100 insertions(+), 57 deletions(-)

diff --git a/aslprep/cli/parser.py b/aslprep/cli/parser.py
index 65a520499..3635916c8 100644
--- a/aslprep/cli/parser.py
+++ b/aslprep/cli/parser.py
@@ -76,8 +76,10 @@ def _bids_filter(value):
         "bids_dir",
         action="store",
         type=PathExists,
-        help="the root folder of a BIDS valid dataset (sub-XXXXX folders should "
-        "be found at the top level in this folder).",
+        help=(
+            "the root folder of a BIDS valid dataset (sub-XXXXX folders should "
+            "be found at the top level in this folder)."
+        ),
     )
     parser.add_argument(
         "output_dir",
@@ -88,8 +90,10 @@ def _bids_filter(value):
     parser.add_argument(
         "analysis_level",
         choices=["participant"],
-        help='processing stage to be run, only "participant" in the case of '
-        "ASLPREP (see BIDS-Apps specification).",
+        help=(
+            'processing stage to be run, only "participant" in the case of '
+            "ASLPREP (see BIDS-Apps specification)."
+        ),
     )
 
     # optional arguments
@@ -109,8 +113,10 @@ def _bids_filter(value):
         action="store",
         nargs="+",
         type=_drop_sub,
-        help="a space delimited list of participant identifiers or a single "
-        "identifier (the sub- prefix can be removed)",
+        help=(
+            "a space delimited list of participant identifiers or a single "
+            "identifier (the sub- prefix can be removed)"
+        ),
     )
     # Re-enable when option is actually implemented
     # g_bids.add_argument('-s', '--session-id', action='store', default='single_session',
@@ -146,8 +152,10 @@ def _bids_filter(value):
         action="store",
         metavar="PATH",
         type=PathExists,
-        help="Reuse the anatomical derivatives from another ASLPrep run or calculated "
-        "with an alternative processing tool (NOT RECOMMENDED).",
+        help=(
+            "Reuse the anatomical derivatives from another ASLPrep run or calculated "
+            "with an alternative processing tool (NOT RECOMMENDED)."
+        ),
     )
 
     g_perfm = parser.add_argument_group("Options to handle performance")
@@ -179,7 +187,7 @@ def _bids_filter(value):
     g_perfm.add_argument(
         "--low-mem",
         action="store_true",
-        help="attempt to reduce memory usage (will increase disk usage " "in working directory)",
+        help="attempt to reduce memory usage (will increase disk usage in working directory)",
     )
     g_perfm.add_argument(
         "--use-plugin",
@@ -189,7 +197,11 @@ def _bids_filter(value):
         type=IsFile,
         help="nipype plugin configuration file",
     )
-    g_perfm.add_argument("--anat-only", action="store_true", help="run anatomical workflows only")
+    g_perfm.add_argument(
+        "--anat-only",
+        action="store_true",
+        help="run anatomical workflows only",
+    )
     g_perfm.add_argument(
         "--boilerplate_only",
         action="store_true",
@@ -219,8 +231,10 @@ def _bids_filter(value):
         nargs="+",
         default=[],
         choices=["fieldmaps", "slicetiming", "sbref"],
-        help="ignore selected aspects of the input dataset to disable corresponding "
-        "parts of the workflow (a space delimited list)",
+        help=(
+            "ignore selected aspects of the input dataset to disable corresponding "
+            "parts of the workflow (a space delimited list)"
+        ),
     )
     g_conf.add_argument(
         "--longitudinal",
@@ -231,16 +245,18 @@ def _bids_filter(value):
         "--output-spaces",
         nargs="*",
         action=OutputReferencesAction,
-        help="""\
-Standard and non-standard spaces to resample anatomical and functional images to. \
-Standard spaces may be specified by the form \
-``<SPACE>[:cohort-<label>][:res-<resolution>][...]``, where ``<SPACE>`` is \
-a keyword designating a spatial reference, and may be followed by optional, \
-colon-separated parameters. \
-Non-standard spaces imply specific orientations and sampling grids. \
-Important to note, the ``res-*`` modifier does not define the resolution used for \
-the spatial normalization. To generate no ASL outputs, use this option without specifying \
-any spatial references.""",
+        help=(
+            "Standard and non-standard spaces to resample anatomical and functional images to. "
+            "Standard spaces may be specified by the form "
+            "``<SPACE>[:cohort-<label>][:res-<resolution>][...]``, where ``<SPACE>`` is "
+            "a keyword designating a spatial reference, and may be followed by optional, "
+            "colon-separated parameters. "
+            "Non-standard spaces imply specific orientations and sampling grids. "
+            "Important to note, the ``res-*`` modifier does not define the resolution used for "
+            "the spatial normalization. "
+            "To generate no ASL outputs, use this option without specifying any spatial "
+            "references."
+        ),
     )
 
     g_conf.add_argument(
@@ -248,8 +264,10 @@ def _bids_filter(value):
         action="store",
         default="register",
         choices=["register", "header"],
-        help='Either "register" (the default) to initialize volumes at center or "header"'
-        " to use the header information when coregistering ASL to T1w images.",
+        help=(
+            'Either "register" (the default) to initialize volumes at center or "header" '
+            "to use the header information when coregistering ASL to T1w images."
+        ),
     )
     g_conf.add_argument(
         "--asl2t1w-dof",
@@ -257,17 +275,21 @@ def _bids_filter(value):
         default=6,
         choices=[6, 9, 12],
         type=int,
-        help="Degrees of freedom when registering ASL to T1w images. "
-        "6 degrees (rotation and translation) are used by default.",
+        help=(
+            "Degrees of freedom when registering ASL to T1w images. "
+            "6 degrees (rotation and translation) are used by default."
+        ),
     )
-    g_conf.add_argument(
+
+    bbr_group = g_conf.add_mutually_exclusive_group()
+    bbr_group.add_argument(
         "--force-bbr",
         action="store_true",
         dest="use_bbr",
         default=False,
         help="Always use boundary-based registration (no goodness-of-fit checks)",
     )
-    g_conf.add_argument(
+    bbr_group.add_argument(
         "--force-no-bbr",
         action="store_false",
         dest="use_bbr",
@@ -296,7 +318,7 @@ def _bids_filter(value):
         action="store",
         default=0,
         type=int,
-        help="Number of initial volumes to ignore",
+        help="Number of initial volumes to ignore.",
     )
     g_conf.add_argument(
         "--smooth_kernel",
@@ -306,23 +328,28 @@ def _bids_filter(value):
         help="Smoothing kernel for the M0 image(s)",
     )
 
-    g_conf.add_argument(
+    # CBF calculation settings
+    cbf_group = g_conf.add_mutually_exclusive_group()
+    cbf_group.add_argument(
         "--scorescrub",
         action="store_true",
         default=False,
-        help=" Sudipto algoritms for denoising CBF",
+        help=(
+            "Apply the Sudipto algorithms (SCRUB and SCORE) for denoising CBF. "
+            "These algorithms will be applied to the standard CBF results."
+        ),
     )
-
-    g_conf.add_argument(
+    cbf_group.add_argument(
         "--basil",
         action="store_true",
         default=False,
-        help=" FSL's CBF computation with spatial regularization and \
-          partial volume correction",
+        help=(
+            "Use FSL's CBF computation with spatial regularization and partial volume correction, "
+            "rather than ASLPrep's standard CBF estimation method."
+        ),
     )
-    # Confounds options
 
-    #  ANTs options
+    # ANTs options
     g_ants = parser.add_argument_group("Specific options for ANTs registrations")
     g_ants.add_argument(
         "--skull-strip-template",
@@ -333,18 +360,22 @@ def _bids_filter(value):
     g_ants.add_argument(
         "--skull-strip-fixed-seed",
         action="store_true",
-        help="do not use a random seed for skull-stripping - will ensure "
-        "run-to-run replicability when used with --omp-nthreads 1 and "
-        "matching --random-seed <int>",
+        help=(
+            "do not use a random seed for skull-stripping - will ensure "
+            "run-to-run replicability when used with --omp-nthreads 1 and "
+            "matching --random-seed <int>"
+        ),
     )
     g_ants.add_argument(
         "--skull-strip-t1w",
         action="store",
         choices=("auto", "skip", "force"),
         default="force",
-        help="determiner for T1-weighted skull stripping ('force' ensures skull "
-        "stripping, 'skip' ignores skull stripping, and 'auto' applies brain extraction "
-        "based on the outcome of a heuristic to check whether the brain is already masked).",
+        help=(
+            "determiner for T1-weighted skull stripping ('force' ensures skull "
+            "stripping, 'skip' ignores skull stripping, and 'auto' applies brain extraction "
+            "based on the outcome of a heuristic to check whether the brain is already masked)."
+        ),
     )
 
     # Fieldmap options
@@ -374,8 +405,10 @@ def _bids_filter(value):
         "--force-syn",
         action="store_true",
         default=False,
-        help="EXPERIMENTAL/TEMPORARY: Use SyN correction in addition to "
-        "fieldmap correction, if available",
+        help=(
+            "EXPERIMENTAL/TEMPORARY: Use SyN correction in addition to "
+            "fieldmap correction, if available"
+        ),
     )
 
     # FreeSurfer options
@@ -384,8 +417,10 @@ def _bids_filter(value):
         "--fs-license-file",
         metavar="FILE",
         type=IsFile,
-        help="Path to FreeSurfer license key file. Get it (for free) by registering"
-        " at https://surfer.nmr.mgh.harvard.edu/registration.html",
+        help=(
+            "Path to FreeSurfer license key file. Get it (for free) by registering "
+            "at https://surfer.nmr.mgh.harvard.edu/registration.html"
+        ),
     )
 
     g_other = parser.add_argument_group("Other options")
@@ -401,8 +436,10 @@ def _bids_filter(value):
         "--clean-workdir",
         action="store_true",
         default=False,
-        help="Clears working directory of contents. Use of this flag is not"
-        "recommended when running concurrent processes of aslprep.",
+        help=(
+            "Clears working directory of contents. Use of this flag is not "
+            "recommended when running concurrent processes of aslprep."
+        ),
     )
     g_other.add_argument(
         "--resource-monitor",
@@ -414,15 +451,19 @@ def _bids_filter(value):
         "--reports-only",
         action="store_true",
         default=False,
-        help="only generate reports, don't run workflows. This will only rerun report "
-        "aggregation, not reportlet generation for specific nodes.",
+        help=(
+            "only generate reports, don't run workflows. This will only rerun report "
+            "aggregation, not reportlet generation for specific nodes."
+        ),
     )
     g_other.add_argument(
         "--run-uuid",
         action="store",
         default=None,
-        help="Specify UUID of previous run, to include error logs in report. "
-        "No effect without --reports-only.",
+        help=(
+            "Specify UUID of previous run, to include error logs in report. "
+            "No effect without --reports-only."
+        ),
     )
     g_other.add_argument(
         "--write-graph",
@@ -440,10 +481,12 @@ def _bids_filter(value):
         "--notrack",
         action="store_true",
         default=False,
-        help="Opt-out of sending tracking information of this run to "
-        "the aslprep developers. This information helps to "
-        "improve aslprep and provides an indicator of real "
-        "world usage crucial for obtaining funding.",
+        help=(
+            "Opt-out of sending tracking information of this run to "
+            "the aslprep developers. This information helps to "
+            "improve aslprep and provides an indicator of real "
+            "world usage crucial for obtaining funding."
+        ),
     )
     g_other.add_argument(
         "--sloppy",

From d5b03ef92985088394c4cf25e8886d01a649ed4d Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Tue, 4 Apr 2023 12:14:34 -0400
Subject: [PATCH 02/27] Rename CBF computation workflows.

---
 aslprep/workflows/asl/base.py  | 52 ++++++++++++++++++----------------
 aslprep/workflows/asl/cbf.py   | 18 ++++++------
 aslprep/workflows/asl/gecbf.py | 48 +++++++++++++++----------------
 3 files changed, 60 insertions(+), 58 deletions(-)

diff --git a/aslprep/workflows/asl/base.py b/aslprep/workflows/asl/base.py
index 292aaa541..f753986d0 100644
--- a/aslprep/workflows/asl/base.py
+++ b/aslprep/workflows/asl/base.py
@@ -17,7 +17,7 @@
 from aslprep.utils.bids import collect_run_data
 from aslprep.utils.meepi import combine_meepi_source
 from aslprep.utils.misc import _create_mem_gb, _get_series_len, _get_wf_name
-from aslprep.workflows.asl.cbf import init_cbf_compt_wf, init_cbfroiquant_wf
+from aslprep.workflows.asl.cbf import init_compute_cbf_wf, init_cbfroiquant_wf
 from aslprep.workflows.asl.confounds import init_asl_confs_wf, init_carpetplot_wf
 from aslprep.workflows.asl.hmc import init_asl_hmc_wf
 from aslprep.workflows.asl.outputs import init_asl_derivatives_wf
@@ -527,8 +527,8 @@ def init_asl_preproc_wf(asl_file):
         ])
         # fmt:on
 
-    # compute  the CBF here
-    compt_cbf_wf = init_cbf_compt_wf(
+    # compute the CBF here
+    compute_cbf_wf = init_compute_cbf_wf(
         name_source=asl_file,
         aslcontext=run_data["aslcontext"],
         dummy_vols=dummyvols,
@@ -537,24 +537,26 @@ def init_asl_preproc_wf(asl_file):
         basil=basil,
         smooth_kernel=smoothkernel,
         metadata=metadata,
-        name="compt_cbf_wf",
+        name="compute_cbf_wf",
     )
 
     # cbf computation workflow
     # fmt:off
     workflow.connect([
-        (inputnode, compt_cbf_wf, [
+        (inputnode, compute_cbf_wf, [
             ("t1w_tpms", "inputnode.t1w_tpms"),
             ("m0scan", "inputnode.m0scan"),
             ("m0scan_metadata", "inputnode.m0scan_metadata"),
+            ("t1w_mask", "inputnode.t1w_mask"),
         ]),
-        (asl_asl_trans_wf, compt_cbf_wf, [
+        (asl_asl_trans_wf, compute_cbf_wf, [
             ("outputnode.asl", "inputnode.asl_file"),
             ("outputnode.asl_mask", "inputnode.asl_mask"),
         ]),
-        (asl_reg_wf, compt_cbf_wf, [("outputnode.itk_t1_to_asl", "inputnode.t1_asl_xform")]),
-        (asl_reg_wf, compt_cbf_wf, [("outputnode.itk_asl_to_t1", "inputnode.itk_asl_to_t1")]),
-        (inputnode, compt_cbf_wf, [("t1w_mask", "inputnode.t1w_mask")]),
+        (asl_reg_wf, compute_cbf_wf, [
+            ("outputnode.itk_t1_to_asl", "inputnode.t1_asl_xform"),
+            ("outputnode.itk_asl_to_t1", "inputnode.itk_asl_to_t1"),
+        ]),
         (asl_reg_wf, outputnode, [
             ("outputnode.itk_t1_to_asl", "itk_t1_to_asl"),
             ("outputnode.itk_asl_to_t1", "itk_asl_to_t1"),
@@ -671,7 +673,7 @@ def init_asl_preproc_wf(asl_file):
             (asl_t1_trans_wf, aslmask_to_t1w, [("outputnode.asl_mask_t1", "reference_image")]),
             (refine_mask, aslmask_to_t1w, [("out_mask", "input_image")]),
             (aslmask_to_t1w, outputnode, [("output_image", "asl_mask_t1")]),
-            (compt_cbf_wf, asl_t1_trans_wf, [
+            (compute_cbf_wf, asl_t1_trans_wf, [
                 ("outputnode.out_cbf", "inputnode.cbf"),
                 ("outputnode.out_mean", "inputnode.meancbf"),
             ]),
@@ -685,10 +687,10 @@ def init_asl_preproc_wf(asl_file):
         if scorescrub:
             # fmt:off
             workflow.connect([
-                (compt_cbf_wf, asl_t1_trans_wf, [
                     ("outputnode.out_score", "inputnode.score"),
                     ("outputnode.out_avgscore", "inputnode.avgscore"),
                     ("outputnode.out_scrub", "inputnode.scrub"),
+                (compute_cbf_wf, asl_t1_trans_wf, [
                 ]),
                 (asl_t1_trans_wf, asl_derivatives_wf, [
                     ("outputnode.scrub_t1", "inputnode.scrub_t1"),
@@ -701,7 +703,7 @@ def init_asl_preproc_wf(asl_file):
         if basil:
             # fmt:off
             workflow.connect([
-                (compt_cbf_wf, asl_t1_trans_wf, [
+                (compute_cbf_wf, asl_t1_trans_wf, [
                     ("outputnode.out_cbfb", "inputnode.basil"),
                     ("outputnode.out_cbfpv", "inputnode.pv"),
                     ("outputnode.out_cbfpvwm", "inputnode.pvwm"),
@@ -726,7 +728,7 @@ def init_asl_preproc_wf(asl_file):
             (refine_mask, asl_derivatives_wf, [
                 ("out_mask", "inputnode.asl_mask_native"),
             ]),
-            (compt_cbf_wf, asl_derivatives_wf, [
+            (compute_cbf_wf, asl_derivatives_wf, [
                 ("outputnode.out_cbf", "inputnode.cbf"),
                 ("outputnode.out_mean", "inputnode.meancbf"),
             ]),
@@ -736,7 +738,7 @@ def init_asl_preproc_wf(asl_file):
         if scorescrub:
             # fmt:off
             workflow.connect([
-                (compt_cbf_wf, asl_derivatives_wf, [
+                (compute_cbf_wf, asl_derivatives_wf, [
                     ("outputnode.out_score", "inputnode.score"),
                     ("outputnode.out_avgscore", "inputnode.avgscore"),
                     ("outputnode.out_scrub", "inputnode.scrub"),
@@ -747,7 +749,7 @@ def init_asl_preproc_wf(asl_file):
         if basil:
             # fmt:off
             workflow.connect([
-                (compt_cbf_wf, asl_derivatives_wf, [
+                (compute_cbf_wf, asl_derivatives_wf, [
                     ("outputnode.out_cbfb", "inputnode.basil"),
                     ("outputnode.out_cbfpv", "inputnode.pv"),
                     ("outputnode.out_cbfpvwm", "inputnode.pvwm"),
@@ -782,7 +784,7 @@ def init_asl_preproc_wf(asl_file):
             ]),
             (refine_mask, asl_std_trans_wf, [("out_mask", "inputnode.asl_mask")]),
             (asl_sdc_wf, asl_std_trans_wf, [("outputnode.out_warp", "inputnode.fieldwarp")]),
-            (compt_cbf_wf, asl_std_trans_wf, [
+            (compute_cbf_wf, asl_std_trans_wf, [
                 ("outputnode.out_cbf", "inputnode.cbf"),
                 ("outputnode.out_mean", "inputnode.meancbf"),
             ]),
@@ -792,7 +794,7 @@ def init_asl_preproc_wf(asl_file):
         if scorescrub:
             # fmt:off
             workflow.connect([
-                (compt_cbf_wf, asl_std_trans_wf, [
+                (compute_cbf_wf, asl_std_trans_wf, [
                     ("outputnode.out_score", "inputnode.score"),
                     ("outputnode.out_avgscore", "inputnode.avgscore"),
                     ("outputnode.out_scrub", "inputnode.scrub"),
@@ -803,7 +805,7 @@ def init_asl_preproc_wf(asl_file):
         if basil:
             # fmt:off
             workflow.connect([
-                (compt_cbf_wf, asl_std_trans_wf, [
+                (compute_cbf_wf, asl_std_trans_wf, [
                     ("outputnode.out_cbfb", "inputnode.basil"),
                     ("outputnode.out_cbfpv", "inputnode.pv"),
                     ("outputnode.out_cbfpvwm", "inputnode.pvwm"),
@@ -878,7 +880,7 @@ def init_asl_preproc_wf(asl_file):
         (inputnode, compt_qccbf_wf, [("t1w_tpms", "inputnode.t1w_tpms")]),
         (asl_reg_wf, compt_qccbf_wf, [("outputnode.itk_t1_to_asl", "inputnode.t1_asl_xform")]),
         (inputnode, compt_qccbf_wf, [("t1w_mask", "inputnode.t1w_mask")]),
-        (compt_cbf_wf, compt_qccbf_wf, [("outputnode.out_mean", "inputnode.meancbf")]),
+        (compute_cbf_wf, compt_qccbf_wf, [("outputnode.out_mean", "inputnode.meancbf")]),
         (asl_confounds_wf, compt_qccbf_wf, [("outputnode.confounds_file", "inputnode.confmat")]),
         (compt_qccbf_wf, outputnode, [("outputnode.qc_file", "qc_file")]),
         (compt_qccbf_wf, asl_derivatives_wf, [("outputnode.qc_file", "inputnode.qc_file")]),
@@ -890,7 +892,7 @@ def init_asl_preproc_wf(asl_file):
     if scorescrub:
         # fmt:off
         workflow.connect([
-            (compt_cbf_wf, compt_qccbf_wf, [
+            (compute_cbf_wf, compt_qccbf_wf, [
                 ("outputnode.out_avgscore", "inputnode.avgscore"),
                 ("outputnode.out_scrub", "inputnode.scrub"),
             ]),
@@ -900,7 +902,7 @@ def init_asl_preproc_wf(asl_file):
     if basil:
         # fmt:off
         workflow.connect([
-            (compt_cbf_wf, compt_qccbf_wf, [
+            (compute_cbf_wf, compt_qccbf_wf, [
                 ("outputnode.out_cbfb", "inputnode.basil"),
                 ("outputnode.out_cbfpv", "inputnode.pv"),
             ]),
@@ -924,7 +926,7 @@ def init_asl_preproc_wf(asl_file):
     )
     # fmt:off
     workflow.connect([
-        (compt_cbf_wf, cbf_plot, [
+        (compute_cbf_wf, cbf_plot, [
             ("outputnode.out_mean", "inputnode.cbf"),
             ("outputnode.out_avgscore", "inputnode.score"),
             ("outputnode.out_scrub", "inputnode.scrub"),
@@ -985,7 +987,7 @@ def init_asl_preproc_wf(asl_file):
         (asl_asl_trans_wf, cbfroiqu, [("outputnode.asl_mask", "inputnode.aslmask")]),
         (inputnode, cbfroiqu, [("std2anat_xfm", "inputnode.std2anat_xfm")]),
         (asl_reg_wf, cbfroiqu, [("outputnode.itk_t1_to_asl", "inputnode.t1_asl_xform")]),
-        (compt_cbf_wf, cbfroiqu, [("outputnode.out_mean", "inputnode.cbf")]),
+        (compute_cbf_wf, cbfroiqu, [("outputnode.out_mean", "inputnode.cbf")]),
         (cbfroiqu, asl_derivatives_wf, [
             ("outputnode.cbf_hvoxf", "inputnode.cbf_hvoxf"),
             ("outputnode.cbf_sc207", "inputnode.cbf_sc207"),
@@ -999,7 +1001,7 @@ def init_asl_preproc_wf(asl_file):
     if scorescrub:
         # fmt:off
         workflow.connect([
-            (compt_cbf_wf, cbfroiqu, [
+            (compute_cbf_wf, cbfroiqu, [
                 ("outputnode.out_avgscore", "inputnode.score"),
                 ("outputnode.out_scrub", "inputnode.scrub"),
             ]),
@@ -1021,7 +1023,7 @@ def init_asl_preproc_wf(asl_file):
     if basil:
         # fmt:off
         workflow.connect([
-            (compt_cbf_wf, cbfroiqu, [
+            (compute_cbf_wf, cbfroiqu, [
                 ("outputnode.out_cbfb", "inputnode.basil"),
                 ("outputnode.out_cbfpv", "inputnode.pvc"),
             ]),
diff --git a/aslprep/workflows/asl/cbf.py b/aslprep/workflows/asl/cbf.py
index b7e863329..f407e3001 100644
--- a/aslprep/workflows/asl/cbf.py
+++ b/aslprep/workflows/asl/cbf.py
@@ -28,7 +28,7 @@
 )
 
 
-def init_cbf_compt_wf(
+def init_compute_cbf_wf(
     name_source,
     aslcontext,
     metadata,
@@ -37,7 +37,7 @@ def init_cbf_compt_wf(
     basil=False,
     M0Scale=1,
     smooth_kernel=5,
-    name="cbf_compt_wf",
+    name="compute_cbf_wf",
 ):
     """Create a workflow for :abbr:`CCBF (compute cbf)`.
 
@@ -46,9 +46,9 @@ def init_cbf_compt_wf(
             :graph2use: orig
             :simple_form: yes
 
-            from aslprep.workflows.asl.cbf import init_cbf_compt_wf
+            from aslprep.workflows.asl.cbf import init_compute_cbf_wf
 
-            wf = init_cbf_compt_wf(
+            wf = init_compute_cbf_wf(
                 name_source="",
                 metadata={},
                 dummy_vols=0,
@@ -68,7 +68,7 @@ def init_cbf_compt_wf(
     M0Scale
     smooth_kernel
     name : :obj:`str`
-        Name of workflow (default: ``cbf_compt_wf``)
+        Name of workflow (default: ``compute_cbf_wf``)
 
     Inputs
     ------
@@ -366,7 +366,7 @@ def _getfiledir(file):
     return workflow
 
 
-def init_gecbf_compt_wf(
+def init_compute_cbf_ge_wf(
     name_source,
     aslcontext,
     metadata,
@@ -374,7 +374,7 @@ def init_gecbf_compt_wf(
     M0Scale=1,
     scorescrub=False,
     basil=False,
-    name="cbf_compt_wf",
+    name="compute_cbf_wf",
 ):
     """Calculate CBF for GE data.
 
@@ -383,9 +383,9 @@ def init_gecbf_compt_wf(
             :graph2use: orig
             :simple_form: yes
 
-            from aslprep.workflows.asl.cbf import init_gecbf_compt_wf
+            from aslprep.workflows.asl.cbf import init_compute_cbf_ge_wf
 
-            wf = init_gecbf_compt_wf(
+            wf = init_compute_cbf_ge_wf(
                 name_source="",
                 metadata={},
                 mem_gb=0.1,
diff --git a/aslprep/workflows/asl/gecbf.py b/aslprep/workflows/asl/gecbf.py
index 19a6e8f34..840d0916b 100644
--- a/aslprep/workflows/asl/gecbf.py
+++ b/aslprep/workflows/asl/gecbf.py
@@ -12,7 +12,7 @@
 from aslprep.niworkflows.interfaces.nibabel import ApplyMask
 from aslprep.utils.bids import collect_run_data
 from aslprep.utils.misc import _create_mem_gb, _get_wf_name
-from aslprep.workflows.asl.cbf import init_cbfroiquant_wf, init_gecbf_compt_wf
+from aslprep.workflows.asl.cbf import init_cbfroiquant_wf, init_compute_cbf_ge_wf
 from aslprep.workflows.asl.ge_utils import (
     init_asl_geref_wf,
     init_asl_gereg_wf,
@@ -353,7 +353,7 @@ def init_asl_gepreproc_wf(asl_file):
     ])
     # fmt:on
 
-    cbf_compt_wf = init_gecbf_compt_wf(
+    compute_cbf_wf = init_compute_cbf_ge_wf(
         name_source=asl_file,
         aslcontext=run_data["aslcontext"],
         metadata=metadata,
@@ -361,22 +361,22 @@ def init_asl_gepreproc_wf(asl_file):
         basil=basil,
         M0Scale=mscale,
         mem_gb=mem_gb["filesize"],
-        name="cbf_compt_wf",
+        name="compute_cbf_wf",
     )
 
     # fmt:off
     workflow.connect([
-        (inputnode, cbf_compt_wf, [
+        (inputnode, compute_cbf_wf, [
             ("asl_file", "inputnode.asl_file"),
             ("t1w_tpms", "inputnode.t1w_tpms"),
             ("t1w_mask", "inputnode.t1w_mask"),
         ]),
-        (gen_ref_wf, cbf_compt_wf, [
+        (gen_ref_wf, compute_cbf_wf, [
             ("outputnode.asl_mask", "inputnode.asl_mask"),
             ("outputnode.m0_file", "inputnode.m0_file"),
             ("outputnode.m0tr", "inputnode.m0tr"),
         ]),
-        (asl_reg_wf, cbf_compt_wf, [
+        (asl_reg_wf, compute_cbf_wf, [
             ("outputnode.itk_asl_to_t1", "inputnode.itk_asl_to_t1"),
             ("outputnode.itk_t1_to_asl", "inputnode.t1_asl_xform"),
         ]),
@@ -414,7 +414,7 @@ def init_asl_gepreproc_wf(asl_file):
             ("outputnode.asl_t1", "asl_t1"),
             ("outputnode.asl_t1_ref", "asl_t1_ref"),
         ]),
-        (cbf_compt_wf, t1w_gereg_wf, [
+        (compute_cbf_wf, t1w_gereg_wf, [
             ("outputnode.out_cbf", "inputnode.cbf"),
             ("outputnode.out_mean", "inputnode.meancbf"),
         ]),
@@ -448,7 +448,7 @@ def init_asl_gepreproc_wf(asl_file):
 
     # fmt:off
     workflow.connect([
-        (cbf_compt_wf, cbf_plot, [("outputnode.out_mean", "inputnode.cbf")]),
+        (compute_cbf_wf, cbf_plot, [("outputnode.out_mean", "inputnode.cbf")]),
         (gen_ref_wf, cbf_plot, [("outputnode.ref_image_brain", "inputnode.asl_ref")]),
     ])
     # fmt:on
@@ -456,7 +456,7 @@ def init_asl_gepreproc_wf(asl_file):
     if scorescrub:
         # fmt:off
         workflow.connect([
-            (cbf_compt_wf, cbf_plot, [
+            (compute_cbf_wf, cbf_plot, [
                 ("outputnode.out_avgscore", "inputnode.score"),
                 ("outputnode.out_scrub", "inputnode.scrub"),
             ]),
@@ -466,7 +466,7 @@ def init_asl_gepreproc_wf(asl_file):
     if basil:
         # fmt:off
         workflow.connect([
-            (cbf_compt_wf, cbf_plot, [
+            (compute_cbf_wf, cbf_plot, [
                 ("outputnode.out_cbfb", "inputnode.basil"),
                 ("outputnode.out_cbfpv", "inputnode.pvc"),
             ]),
@@ -496,7 +496,7 @@ def init_asl_gepreproc_wf(asl_file):
         if scorescrub:
             # fmt:off
             workflow.connect([
-                (cbf_compt_wf, t1w_gereg_wf, [
+                (compute_cbf_wf, t1w_gereg_wf, [
                     ("outputnode.out_score", "inputnode.score"),
                     ("outputnode.out_avgscore", "inputnode.avgscore"),
                     ("outputnode.out_scrub", "inputnode.scrub"),
@@ -512,7 +512,7 @@ def init_asl_gepreproc_wf(asl_file):
         if basil:
             # fmt:off
             workflow.connect([
-                (cbf_compt_wf, t1w_gereg_wf, [
+                (compute_cbf_wf, t1w_gereg_wf, [
                     ("outputnode.out_cbfb", "inputnode.basil"),
                     ("outputnode.out_cbfpv", "inputnode.pv"),
                     ("outputnode.out_cbfpvwm", "inputnode.pvwm"),
@@ -542,7 +542,7 @@ def init_asl_gepreproc_wf(asl_file):
         ]),
         (refine_mask, compt_qccbf_wf, [("out_mask", "inputnode.asl_mask")]),
         (asl_reg_wf, compt_qccbf_wf, [("outputnode.itk_t1_to_asl", "inputnode.t1_asl_xform")]),
-        (cbf_compt_wf, compt_qccbf_wf, [("outputnode.out_mean", "inputnode.meancbf")]),
+        (compute_cbf_wf, compt_qccbf_wf, [("outputnode.out_mean", "inputnode.meancbf")]),
         (compt_qccbf_wf, outputnode, [("outputnode.qc_file", "qc_file")]),
         (compt_qccbf_wf, asl_derivatives_wf, [("outputnode.qc_file", "inputnode.qc_file")]),
         (compt_qccbf_wf, summary, [("outputnode.qc_file", "qc_file")]),
@@ -552,7 +552,7 @@ def init_asl_gepreproc_wf(asl_file):
     if scorescrub:
         # fmt:off
         workflow.connect([
-            (cbf_compt_wf, compt_qccbf_wf, [
+            (compute_cbf_wf, compt_qccbf_wf, [
                 ("outputnode.out_avgscore", "inputnode.avgscore"),
                 ("outputnode.out_scrub", "inputnode.scrub"),
             ]),
@@ -562,7 +562,7 @@ def init_asl_gepreproc_wf(asl_file):
     if basil:
         # fmt:off
         workflow.connect([
-            (cbf_compt_wf, compt_qccbf_wf, [
+            (compute_cbf_wf, compt_qccbf_wf, [
                 ("outputnode.out_cbfb", "inputnode.basil"),
                 ("outputnode.out_cbfpv", "inputnode.pv"),
             ]),
@@ -577,7 +577,7 @@ def init_asl_gepreproc_wf(asl_file):
                 ("outputnode.raw_ref_image", "inputnode.asl_native_ref"),
             ]),
             (refine_mask, asl_derivatives_wf, [("out_mask", "inputnode.asl_mask_native")]),
-            (cbf_compt_wf, asl_derivatives_wf, [
+            (compute_cbf_wf, asl_derivatives_wf, [
                 ("outputnode.out_cbf", "inputnode.cbf"),
                 ("outputnode.out_mean", "inputnode.meancbf"),
             ]),
@@ -587,7 +587,7 @@ def init_asl_gepreproc_wf(asl_file):
         if scorescrub:
             # fmt:off
             workflow.connect([
-                (cbf_compt_wf, asl_derivatives_wf, [
+                (compute_cbf_wf, asl_derivatives_wf, [
                     ("outputnode.out_score", "inputnode.score"),
                     ("outputnode.out_avgscore", "inputnode.avgscore"),
                     ("outputnode.out_scrub", "inputnode.scrub"),
@@ -598,7 +598,7 @@ def init_asl_gepreproc_wf(asl_file):
         if basil:
             # fmt:off
             workflow.connect([
-                (cbf_compt_wf, asl_derivatives_wf, [
+                (compute_cbf_wf, asl_derivatives_wf, [
                     ("outputnode.out_cbfb", "inputnode.basil"),
                     ("outputnode.out_cbfpv", "inputnode.pv"),
                     ("outputnode.out_cbfpvwm", "inputnode.pvwm"),
@@ -633,7 +633,7 @@ def init_asl_gepreproc_wf(asl_file):
                 ("outputnode.asl_std_ref", "asl_std_ref"),
                 ("outputnode.asl_mask_std", "asl_mask_std"),
             ]),
-            (cbf_compt_wf, std_gereg_wf, [
+            (compute_cbf_wf, std_gereg_wf, [
                 ("outputnode.out_cbf", "inputnode.cbf"),
                 ("outputnode.out_mean", "inputnode.meancbf"),
             ]),
@@ -643,7 +643,7 @@ def init_asl_gepreproc_wf(asl_file):
         if scorescrub:
             # fmt:off
             workflow.connect([
-                (cbf_compt_wf, std_gereg_wf, [
+                (compute_cbf_wf, std_gereg_wf, [
                     ("outputnode.out_score", "inputnode.score"),
                     ("outputnode.out_avgscore", "inputnode.avgscore"),
                     ("outputnode.out_scrub", "inputnode.scrub"),
@@ -654,7 +654,7 @@ def init_asl_gepreproc_wf(asl_file):
         if basil:
             # fmt:off
             workflow.connect([
-                (cbf_compt_wf, std_gereg_wf, [
+                (compute_cbf_wf, std_gereg_wf, [
                     ("outputnode.out_cbfb", "inputnode.basil"),
                     ("outputnode.out_cbfpv", "inputnode.pv"),
                     ("outputnode.out_cbfpvwm", "inputnode.pvwm"),
@@ -721,7 +721,7 @@ def init_asl_gepreproc_wf(asl_file):
         (refine_mask, cbfroiqu, [("out_mask", "inputnode.aslmask")]),
         (inputnode, cbfroiqu, [("std2anat_xfm", "inputnode.std2anat_xfm")]),
         (asl_reg_wf, cbfroiqu, [("outputnode.itk_t1_to_asl", "inputnode.t1_asl_xform")]),
-        (cbf_compt_wf, cbfroiqu, [("outputnode.out_mean", "inputnode.cbf")]),
+        (compute_cbf_wf, cbfroiqu, [("outputnode.out_mean", "inputnode.cbf")]),
         (cbfroiqu, asl_derivatives_wf, [
             ("outputnode.cbf_hvoxf", "inputnode.cbf_hvoxf"),
             ("outputnode.cbf_sc207", "inputnode.cbf_sc207"),
@@ -735,7 +735,7 @@ def init_asl_gepreproc_wf(asl_file):
     if scorescrub:
         # fmt:off
         workflow.connect([
-            (cbf_compt_wf, cbfroiqu, [
+            (compute_cbf_wf, cbfroiqu, [
                 ("outputnode.out_avgscore", "inputnode.score"),
                 ("outputnode.out_scrub", "inputnode.scrub"),
             ]),
@@ -757,7 +757,7 @@ def init_asl_gepreproc_wf(asl_file):
     if basil:
         # fmt:off
         workflow.connect([
-            (cbf_compt_wf, cbfroiqu, [
+            (compute_cbf_wf, cbfroiqu, [
                 ("outputnode.out_cbfb", "inputnode.basil"),
                 ("outputnode.out_cbfpv", "inputnode.pvc"),
             ]),

From 8acc889cddb178367e11a5332c555da58ff706c0 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Tue, 4 Apr 2023 12:16:14 -0400
Subject: [PATCH 03/27] Clean up SCORE code.

---
 aslprep/interfaces/cbf_computation.py | 148 +++++++++++++-------------
 aslprep/utils/misc.py                 | 118 ++++++++++++--------
 2 files changed, 146 insertions(+), 120 deletions(-)

diff --git a/aslprep/interfaces/cbf_computation.py b/aslprep/interfaces/cbf_computation.py
index aba45deeb..6483971e2 100644
--- a/aslprep/interfaces/cbf_computation.py
+++ b/aslprep/interfaces/cbf_computation.py
@@ -22,7 +22,7 @@
 
 from aslprep import config
 from aslprep.utils.misc import (
-    _getcbfscore,
+    _score_cbf,
     _scrubcbf,
     estimate_labeling_efficiency,
     parcellate_cbf,
@@ -597,32 +597,37 @@ def _run_interface(self, runtime):
 
 
 class _ScoreAndScrubCBFInputSpec(BaseInterfaceInputSpec):
-    in_file = File(exists=True, mandatory=True, desc="computed CBF from ComputeCBF")
-    in_greyM = File(exists=True, mandatory=True, desc="grey  matter")
-    in_whiteM = File(exists=True, mandatory=True, desc="white  matter")
-    in_mask = File(exists=True, mandatory=True, desc="mask")
-    in_csf = File(exists=True, mandatory=True, desc="csf")
-    in_thresh = traits.Float(
-        default_value=0.7, exists=True, mandatory=False, desc="threshold of propbaility matter"
+    cbf = File(exists=True, mandatory=True, desc="computed CBF from ComputeCBF")
+    gm_tpm = File(exists=True, mandatory=True, desc="Gray matter tissue probability map.")
+    wm_tpm = File(exists=True, mandatory=True, desc="White matter tissue probability map.")
+    csf_tpm = File(exists=True, mandatory=True, desc="Cerebrospinal fluid tissue probability map.")
+    brain_mask = File(exists=True, mandatory=True, desc="Brain mask.")
+    tpm_threshold = traits.Float(
+        default_value=0.7,
+        mandatory=False,
+        usedefault=True,
+        desc="Threshold for tissue probability maps.",
     )
-    in_wfun = traits.Str(
+    wavelet_function = traits.Str(
         exists=True,
         mandatory=False,
         default_value="huber",
         option=["bisquare", "andrews", "cauchy", "fair", "logistics", "ols", "talwar", "welsch"],
-        desc="wavelet fun ",
+        desc="Wavelet function used in SCRUB.",
     )
-    out_score = File(exists=False, desc="score timeseries data")
-    out_avgscore = File(exists=False, desc="average score")
-    out_scrub = File(exists=False, desc="average scrub")
-    out_scoreindex = File(exists=False, desc="index of volume remove or leave by score")
 
 
 class _ScoreAndScrubCBFOutputSpec(TraitedSpec):
-    out_score = File(exists=False, mandatory=False, desc="score timeseries data")
-    out_avgscore = File(exists=False, mandatory=False, desc="average score")
-    out_scrub = File(exists=False, mandatory=False, desc="average scrub")
-    out_scoreindex = File(exists=False, mandatory=False, desc="index of volume remove ")
+    cbf_ts_score = File(
+        exists=True,
+        desc="CBF time series after removing outlier volumes flagged by SCORE algorithm.",
+    )
+    score_outliers = File(exists=True, desc="Index of removed volumes, in a CSV file.")
+    mean_cbf_score = File(
+        exists=True,
+        desc="Mean CBF image calculated from SCORE-censored CBF time series.",
+    )
+    mean_cbf_scrub = File(exists=True, desc="average scrub")
 
 
 class ScoreAndScrubCBF(SimpleInterface):
@@ -644,84 +649,81 @@ class ScoreAndScrubCBF(SimpleInterface):
     output_spec = _ScoreAndScrubCBFOutputSpec
 
     def _run_interface(self, runtime):
-        cbf_ts = nb.load(self.inputs.in_file).get_fdata()
-        mask = nb.load(self.inputs.in_mask).get_fdata()
-        greym = nb.load(self.inputs.in_greyM).get_fdata()
-        whitem = nb.load(self.inputs.in_whiteM).get_fdata()
-        csf = nb.load(self.inputs.in_csf).get_fdata()
-        if cbf_ts.ndim > 3:
-            cbf_scorets, index_score = _getcbfscore(
-                cbfts=cbf_ts,
-                wm=whitem,
-                gm=greym,
-                csf=csf,
-                mask=mask,
-                thresh=self.inputs.in_thresh,
+        cbf_img = nb.load(self.inputs.cbf)
+        cbf_arr = cbf_img.get_fdata()
+
+        if cbf_arr.ndim > 3:
+            mask_arr = nb.load(self.inputs.brain_mask).get_fdata()
+            gm_tpm_arr = nb.load(self.inputs.gm_tpm).get_fdata()
+            wm_tpm_arr = nb.load(self.inputs.wm_tpm).get_fdata()
+            csf_tpm_arr = nb.load(self.inputs.csf_tpm).get_fdata()
+
+            cbf_ts_score, score_outliers_idx = _score_cbf(
+                cbf_ts=cbf_arr,
+                gm=gm_tpm_arr,
+                wm=wm_tpm_arr,
+                csf=csf_tpm_arr,
+                mask=mask_arr,
+                thresh=self.inputs.tpm_threshold,
             )
-            cbfscrub = _scrubcbf(
-                cbf_ts=cbf_scorets,
-                gm=greym,
-                wm=whitem,
-                csf=csf,
-                mask=mask,
-                wfun=self.inputs.in_wfun,
-                thresh=self.inputs.in_thresh,
+            mean_cbf_score = np.mean(cbf_ts_score, axis=3)
+            mean_cbf_scrub = _scrubcbf(
+                cbf_ts=cbf_ts_score,
+                gm=gm_tpm_arr,
+                wm=wm_tpm_arr,
+                csf=csf_tpm_arr,
+                mask=mask_arr,
+                wavelet_function=self.inputs.wavelet_function,
+                thresh=self.inputs.tpm_threshold,
             )
-            avgscore = np.mean(cbf_scorets, axis=3)
         else:
             config.loggers.interface.warning(
-                f"CBF time series is only {cbf_ts.ndim}D. Skipping SCORE and SCRUB."
+                f"CBF time series is only {cbf_arr.ndim}D. Skipping SCORE and SCRUB."
             )
-            cbf_scorets = cbf_ts
-            index_score = np.array([0])
-            cbfscrub = cbf_ts
-            avgscore = cbf_ts
+            cbf_ts_score = cbf_arr
+            mean_cbf_score = cbf_arr
+            score_outliers_idx = np.array([0])
+            mean_cbf_scrub = cbf_arr
 
-        self._results["out_score"] = fname_presuffix(
-            self.inputs.in_file,
+        self._results["cbf_ts_score"] = fname_presuffix(
+            self.inputs.cbf,
             suffix="_cbfscorets",
             newpath=runtime.cwd,
         )
-        self._results["out_avgscore"] = fname_presuffix(
-            self.inputs.in_file,
+        self._results["mean_cbf_score"] = fname_presuffix(
+            self.inputs.cbf,
             suffix="_meancbfscore",
             newpath=runtime.cwd,
         )
-        self._results["out_scrub"] = fname_presuffix(
-            self.inputs.in_file,
+        self._results["mean_cbf_scrub"] = fname_presuffix(
+            self.inputs.cbf,
             suffix="_cbfscrub",
             newpath=runtime.cwd,
         )
-        self._results["out_scoreindex"] = fname_presuffix(
-            self.inputs.in_file,
-            suffix="_scoreindex.txt",
+        self._results["score_outliers"] = fname_presuffix(
+            self.inputs.cbf,
+            suffix="_scoreindex.csv",
             newpath=runtime.cwd,
             use_ext=False,
         )
-        samplecbf = nb.load(self.inputs.in_mask)
 
         nb.Nifti1Image(
-            dataobj=cbf_scorets,
-            affine=samplecbf.affine,
-            header=samplecbf.header,
-        ).to_filename(self._results["out_score"])
+            dataobj=cbf_ts_score,
+            affine=cbf_img.affine,
+            header=cbf_img.header,
+        ).to_filename(self._results["cbf_ts_score"])
         nb.Nifti1Image(
-            dataobj=avgscore,
-            affine=samplecbf.affine,
-            header=samplecbf.header,
-        ).to_filename(self._results["out_avgscore"])
+            dataobj=mean_cbf_score,
+            affine=cbf_img.affine,
+            header=cbf_img.header,
+        ).to_filename(self._results["mean_cbf_score"])
+        np.savetxt(self._results["score_outliers"], score_outliers_idx, delimiter=",")
         nb.Nifti1Image(
-            dataobj=cbfscrub,
-            affine=samplecbf.affine,
-            header=samplecbf.header,
-        ).to_filename(self._results["out_scrub"])
-
-        np.savetxt(self._results["out_scoreindex"], index_score, delimiter=",")
+            dataobj=mean_cbf_scrub,
+            affine=cbf_img.affine,
+            header=cbf_img.header,
+        ).to_filename(self._results["mean_cbf_scrub"])
 
-        self.inputs.out_score = os.path.abspath(self._results["out_score"])
-        self.inputs.out_avgscore = os.path.abspath(self._results["out_avgscore"])
-        self.inputs.out_scrub = os.path.abspath(self._results["out_scrub"])
-        self.inputs.out_scoreindex = os.path.abspath(self._results["out_scoreindex"])
         return runtime
 
 
diff --git a/aslprep/utils/misc.py b/aslprep/utils/misc.py
index d49d4c382..5bb88a5de 100644
--- a/aslprep/utils/misc.py
+++ b/aslprep/utils/misc.py
@@ -518,67 +518,91 @@ def _getchisquare(n):
     return a[n - 1], b[n - 1]
 
 
-def _getcbfscore(cbfts, wm, gm, csf, mask, thresh=0.7):
+def _score_cbf(cbf_ts, wm, gm, csf, mask, thresh=0.7):
     """Apply SCORE algorithm to remove noisy CBF volumes.
 
     Parameters
     ----------
-    cbf_ts
-       nd array of 3D or 4D computed cbf
-    gm,wm,csf
-       numpy array of grey matter, whitematter, and csf
-    mask
-       numpy array of mask
+    cbf_ts : 4D numpy.ndarray
+        CBF time series.
+    gm, wm, csf : 3D numpy.ndarray
+        Gray matter, white matter, and CSF tissue probability maps.
+    mask : 3D numpy.ndarray
+        Binary brain mask.
+    thresh : float
+        Tissue probability threshold.
+        Voxels in gm, wm, and csf with values greater than or equal to the threshold are
+        considered to be of that tissue.
+
+    Returns
+    -------
+    cbf_ts_good_vols : 4D numpy.ndarray
+        CBF time series containing only non-outlier volumes. This time series is also brain-masked.
+    outlier_idx : 1D numpy.ndarray
+        An index of outlier volumes in the CBF time series. Values of 1 indicate outliers.
     """
+    n_volumes = cbf_ts.shape[3]
+
+    # Create binary tissue-type masks
     gm_bin = gm >= thresh
     wm_bin = wm >= thresh
     csf_bin = csf >= thresh
+    tissue_type_masks = [gm_bin, wm_bin, csf_bin]
 
-    # get the total number of voxle within csf,gm and wm
+    # Get the number of voxels within each tissue type, then subtract 1 to get Nk-1.
     n_gm_voxels = np.sum(gm_bin) - 1
     n_wm_voxels = np.sum(wm_bin) - 1
     n_csf_voxels = np.sum(csf_bin) - 1
-    mask1 = gm_bin + wm_bin + csf_bin
-
-    # mean  of times series cbf within greymatter
-    gm_cbf_ts = np.squeeze(np.mean(cbfts[gm_bin, :], axis=0))
-
-    # robust mean and median
-    median_gm_cbf = np.median(gm_cbf_ts)
-    mad_gm_cbf = median_abs_deviation(gm_cbf_ts) / 0.675
-    indx = 1 * (np.abs(gm_cbf_ts - median_gm_cbf) > (2.5 * mad_gm_cbf))
-    R = np.mean(cbfts[:, :, :, indx == 0], axis=3)
-    V = (
-        n_gm_voxels * np.var(R[gm == 1])
-        + n_wm_voxels * np.var(R[wm == 1])
-        + n_csf_voxels * np.var(R[csf == 1])
-    )
-    V1 = V + 1
-    while V < V1:
-        V1 = V
-        CC = np.zeros(cbfts.shape[3]) * (-2)
-        for s in range(cbfts.shape[3]):
-            if indx[s] != 0:
+    n_voxels = [n_gm_voxels, n_wm_voxels, n_csf_voxels]
+
+    # Create a mask only including GM, WM, and CSF
+    tissue_type_mask = (gm_bin + wm_bin + csf_bin).astype(bool)
+
+    # Perform initial outlier exclusion.
+    # Remove volumes with mean GM-CBF outside 2.5 SDs of their means.
+    mean_gb_cbf = np.squeeze(np.mean(cbf_ts[gm_bin, :], axis=0))
+    median_gm_cbf = np.median(mean_gb_cbf)
+    mad_gm_cbf = median_abs_deviation(mean_gb_cbf) / 0.675
+    outlier_idx = np.abs(mean_gb_cbf - median_gm_cbf) > (2.5 * mad_gm_cbf)
+
+    iter_mean_cbf = np.mean(cbf_ts[:, :, :, ~outlier_idx], axis=3)
+    iter_pooled_variance = [
+        (n_voxels[k] * np.var(iter_mean_cbf[tissue_type_masks[k]])) / n_voxels[k] for k in range(3)
+    ]
+    last_pooled_variance = iter_pooled_variance + 1  # add 1 to ensure it is >
+    while iter_pooled_variance < last_pooled_variance:
+        # If variance increases from the previous iteration, we interpret this as removing
+        # random noise instead of artifact.
+        last_pooled_variance = iter_pooled_variance
+        volume_wise_correlations = np.full(n_volumes, np.nan)
+        for i_vol in range(n_volumes):
+            if outlier_idx[i_vol]:
                 break
-            else:
-                tmp1 = cbfts[:, :, :, s]
-                CC[s] = np.corrcoef(R[mask1 > 0], tmp1[mask1 > 0])[0][1]
-
-        inx = np.argmax(CC)
-        indx[inx] = 2
-        R = np.mean(cbfts[:, :, :, indx == 0], axis=3)
-        V = (
-            (n_gm_voxels * np.var(R[gm == 1]))
-            + (n_wm_voxels * np.var(R[wm == 1]))
-            + (n_csf_voxels * np.var(R[csf == 1]))
-        )
-    cbfts_recon = cbfts[:, :, :, indx == 0]
-    cbfts_recon1 = np.zeros_like(cbfts_recon)
-    for i in range(cbfts_recon.shape[3]):
-        cbfts_recon1[:, :, :, i] = cbfts_recon[:, :, :, i] * mask
 
-    cbfts_recon1 = np.nan_to_num(cbfts_recon1)
-    return cbfts_recon1, indx
+            cbf_vol = cbf_ts[:, :, :, i_vol]
+            # Correlate mean CBF's masked values with volume CBF's masked values.
+            volume_wise_correlations[i_vol] = np.corrcoef(
+                iter_mean_cbf[tissue_type_mask],
+                cbf_vol[tissue_type_mask],
+            )[0, 1]
+
+        # The volume with the maximum correlation value is this iteration's outlier.
+        max_correlation_idx = np.nanargmax(volume_wise_correlations)
+        outlier_idx[max_correlation_idx] = True
+
+        # Recalculate mean CBF map and pooled variance.
+        iter_mean_cbf = np.mean(cbf_ts[:, :, :, ~outlier_idx], axis=3)
+        iter_pooled_variance = [
+            (n_voxels[k] * np.var(iter_mean_cbf[tissue_type_masks[k]])) / n_voxels[k]
+            for k in range(3)
+        ]
+
+    # Create CBF time series from just the non-outlier volumes and brain-mask it.
+    cbf_ts_good_vols = cbf_ts[:, :, :, ~outlier_idx]
+    cbf_ts_good_vols *= mask[..., None]
+    cbf_ts_good_vols = np.nan_to_num(cbf_ts_good_vols)
+
+    return cbf_ts_good_vols, outlier_idx.astype(int)
 
 
 def _robust_fit(

From 337274d4b7ea2137100c135f626aad3759e22032 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Tue, 4 Apr 2023 12:16:23 -0400
Subject: [PATCH 04/27] Fix misc. bugs.

---
 aslprep/interfaces/cbf_computation.py | 4 ++--
 docs/workflows.rst                    | 6 +++---
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/aslprep/interfaces/cbf_computation.py b/aslprep/interfaces/cbf_computation.py
index 6483971e2..389d3f696 100644
--- a/aslprep/interfaces/cbf_computation.py
+++ b/aslprep/interfaces/cbf_computation.py
@@ -99,11 +99,11 @@ class _ExtractCBFInputSpec(BaseInterfaceInputSpec):
     in_mask = File(exists=True, mandatory=True, desc="mask")
     dummy_vols = traits.Int(
         default_value=0,
-        use_default=True,
+        usedefault=True,
         mandatory=False,
         desc="remove first n volumes",
     )
-    fwhm = traits.Float(default_value=5, use_default=True, mandatory=False, desc="fwhm")
+    fwhm = traits.Float(default_value=5, usedefault=True, mandatory=False, desc="fwhm")
 
 
 class _ExtractCBFOutputSpec(TraitedSpec):
diff --git a/docs/workflows.rst b/docs/workflows.rst
index 4926a2d64..df26b98ab 100644
--- a/docs/workflows.rst
+++ b/docs/workflows.rst
@@ -276,7 +276,7 @@ Interpolation uses a Lanczos kernel.
 CBF Computation in native space
 *******************************
 
-:py:func:`~aslprep.workflows.asl.cbf.init_cbf_compt_wf`
+:py:func:`~aslprep.workflows.asl.cbf.init_compute_cbf_wf`
 
 .. workflow::
     :graph2use: orig
@@ -286,14 +286,14 @@ CBF Computation in native space
     from pathlib import Path
     from pkg_resources import resource_filename as pkgrf
 
-    from aslprep.workflows.asl.cbf import init_cbf_compt_wf
+    from aslprep.workflows.asl.cbf import init_compute_cbf_wf
 
     bids_dir = Path(pkgrf('aslprep', 'tests/data/ds000240')).absolute()
     metadata_file = bids_dir / 'sub-01' / 'perf'/ 'sub-01_asl.json'
     with open(metadata_file) as f:
         metadata = json.load(f)
 
-    wf = init_cbf_compt_wf(
+    wf = init_compute_cbf_wf(
         bids_dir=str(bids_dir),
         scorescrub=False,
         basil=False,

From adb7ace83b77e23ce9e46d2fae3c7184bf57a3d2 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Tue, 4 Apr 2023 12:17:04 -0400
Subject: [PATCH 05/27] Rename SCORE/SCRUB outputs.

---
 aslprep/workflows/asl/base.py |  12 +--
 aslprep/workflows/asl/cbf.py  | 172 +++++++++++++++++-----------------
 2 files changed, 93 insertions(+), 91 deletions(-)

diff --git a/aslprep/workflows/asl/base.py b/aslprep/workflows/asl/base.py
index f753986d0..14634e02c 100644
--- a/aslprep/workflows/asl/base.py
+++ b/aslprep/workflows/asl/base.py
@@ -687,15 +687,15 @@ def init_asl_preproc_wf(asl_file):
         if scorescrub:
             # fmt:off
             workflow.connect([
-                    ("outputnode.out_score", "inputnode.score"),
-                    ("outputnode.out_avgscore", "inputnode.avgscore"),
-                    ("outputnode.out_scrub", "inputnode.scrub"),
                 (compute_cbf_wf, asl_t1_trans_wf, [
+                    ("outputnode.cbf_ts_score", "inputnode.cbf_ts_score"),
+                    ("outputnode.mean_cbf_score", "inputnode.mean_cbf_score"),
+                    ("outputnode.mean_cbf_scrub", "inputnode.mean_cbf_scrub"),
                 ]),
                 (asl_t1_trans_wf, asl_derivatives_wf, [
-                    ("outputnode.scrub_t1", "inputnode.scrub_t1"),
-                    ("outputnode.score_t1", "inputnode.score_t1"),
-                    ("outputnode.avgscore_t1", "inputnode.avgscore_t1"),
+                    ("outputnode.cbf_ts_score_t1", "inputnode.cbf_ts_score_t1"),
+                    ("outputnode.mean_cbf_score_t1", "inputnode.mean_cbf_score_t1"),
+                    ("outputnode.mean_cbf_scrub_t1", "inputnode.mean_cbf_scrub_t1"),
                 ]),
             ])
             # fmt:on
diff --git a/aslprep/workflows/asl/cbf.py b/aslprep/workflows/asl/cbf.py
index f407e3001..ce099527c 100644
--- a/aslprep/workflows/asl/cbf.py
+++ b/aslprep/workflows/asl/cbf.py
@@ -101,8 +101,6 @@ def init_compute_cbf_wf(
 
 *ASLPrep* was configured to calculate cerebral blood flow (CBF) using the following methods.
 
-The cerebral blood flow (CBF) was quantified from preprocessed ASL data using a general kinetic
-model [@buxton1998general].
 """
     inputnode = pe.Node(
         niu.IdentityInterface(
@@ -123,14 +121,14 @@ def init_compute_cbf_wf(
     outputnode = pe.Node(
         niu.IdentityInterface(
             fields=[
-                "out_cbf",
-                "out_mean",
-                "out_score",
+                "cbf_ts",
+                "mean_cbf",
+                "cbf_ts_score",
+                "score_outliers",
+                "mean_cbf_score",
+                "mean_cbf_scrub",
                 "out_cbfpvwm",
-                "out_avgscore",
-                "out_scrub",
                 "out_cbfb",
-                "out_scoreindex",
                 "out_cbfpv",
                 "out_att",
             ]
@@ -256,67 +254,9 @@ def _getfiledir(file):
     ])
     # fmt:on
 
-    compute_cbf = pe.Node(
-        ComputeCBF(
-            cbf_only=cbf_only,
-            m0scale=M0Scale,
-        ),
-        mem_gb=0.2,
-        run_without_submitting=True,
-        name="compute_cbf",
-    )
-
-    # fmt:off
-    workflow.connect([
-        (refine_mask, compute_cbf, [("out_mask", "mask")]),
-        (extract_deltam, compute_cbf, [
-            ("out_file", "deltam"),
-            ("m0_file", "m0_file"),
-            ("metadata", "metadata"),
-        ]),
-        (compute_cbf, outputnode, [
-            ("cbf", "out_cbf"),
-            ("mean_cbf", "out_mean"),
-        ]),
-    ])
-    # fmt:on
-
-    if scorescrub:
-        score_and_scrub_cbf = pe.Node(
-            ScoreAndScrubCBF(in_thresh=0.7, in_wfun="huber"),
-            mem_gb=0.2,
-            name="score_and_scrub_cbf",
-            run_without_submitting=True,
-        )
-
-        workflow.__desc__ += """
-
-Structural Correlation based Outlier Rejection (SCORE) algorithm was applied to the CBF timeseries
-to discard CBF volumes with outlying values [@dolui2017structural] before computing the mean CBF.
-Following SCORE, the Structural Correlation with RobUst Bayesian (SCRUB) algorithm was applied to
-the CBF maps using structural tissue probability maps to reweight the mean CBF
-[@dolui2017structural;@dolui2016scrub].
-"""
-        # fmt:off
-        workflow.connect([
-            (refine_mask, score_and_scrub_cbf, [("out_mask", "in_mask")]),
-            (compute_cbf, score_and_scrub_cbf, [("cbf", "in_file")]),
-            (gm_tfm, score_and_scrub_cbf, [("output_image", "in_greyM")]),
-            (wm_tfm, score_and_scrub_cbf, [("output_image", "in_whiteM")]),
-            (csf_tfm, score_and_scrub_cbf, [("output_image", "in_csf")]),
-            (score_and_scrub_cbf, outputnode, [
-                ("out_score", "out_score"),
-                ("out_scoreindex", "out_scoreindex"),
-                ("out_avgscore", "out_avgscore"),
-                ("out_scrub", "out_scrub"),
-            ]),
-        ])
-        # fmt:on
-
     if basil:
         workflow.__desc__ += f"""
-
-CBF was also computed with Bayesian Inference for Arterial Spin Labeling (BASIL)
+CBF was computed with Bayesian Inference for Arterial Spin Labeling (BASIL)
 [@chappell2008variational], as implemented in *FSL* {Info.version().split(':')[0]}.
 BASIL computes CBF using a spatial regularization of the estimated perfusion image and
 additionally calculates a partial-volume corrected CBF image [@chappell_pvc].
@@ -355,14 +295,76 @@ def _getfiledir(file):
             (gm_tfm, basilcbf, [("output_image", "pvgm")]),
             (wm_tfm, basilcbf, [("output_image", "pvwm")]),
             (basilcbf, outputnode, [
-                ("out_cbfb", "out_cbfb"),
-                ("out_cbfpv", "out_cbfpv"),
-                ("out_cbfpvwm", "out_cbfpvwm"),
-                ("out_att", "out_att"),
+                ("out_cbfb", "cbf"),
+                ("out_cbfpv", "cbf_pv"),
+                ("out_cbfpvwm", "cbf_pv_wm"),
+                ("out_att", "att"),
+            ]),
+        ])
+        # fmt:on
+
+    else:
+        workflow.__desc__ += """
+The cerebral blood flow (CBF) was quantified from preprocessed ASL data using a general kinetic
+model [@buxton1998general].
+"""
+        compute_cbf = pe.Node(
+            ComputeCBF(
+                cbf_only=cbf_only,
+                m0scale=M0Scale,
+            ),
+            mem_gb=0.2,
+            run_without_submitting=True,
+            name="compute_cbf",
+        )
+
+        # fmt:off
+        workflow.connect([
+            (refine_mask, compute_cbf, [("out_mask", "mask")]),
+            (extract_deltam, compute_cbf, [
+                ("out_file", "deltam"),
+                ("m0_file", "m0_file"),
+                ("metadata", "metadata"),
+            ]),
+            (compute_cbf, outputnode, [
+                ("cbf", "out_cbf"),
+                ("mean_cbf", "out_mean"),
             ]),
         ])
         # fmt:on
 
+        if scorescrub:
+            score_and_scrub_cbf = pe.Node(
+                ScoreAndScrubCBF(tpm_threshold=0.7, wavelet_function="huber"),
+                mem_gb=0.2,
+                name="score_and_scrub_cbf",
+                run_without_submitting=True,
+            )
+
+            workflow.__desc__ += """
+
+Structural Correlation based Outlier Rejection (SCORE) algorithm was applied to the CBF timeseries
+to discard CBF volumes with outlying values [@dolui2017structural] before computing the mean CBF.
+Following SCORE, the Structural Correlation with RobUst Bayesian (SCRUB) algorithm was applied to
+the CBF maps using structural tissue probability maps to reweight the mean CBF
+[@dolui2017structural;@dolui2016scrub].
+"""
+            # fmt:off
+            workflow.connect([
+                (refine_mask, score_and_scrub_cbf, [("out_mask", "brain_mask")]),
+                (compute_cbf, score_and_scrub_cbf, [("cbf", "cbf")]),
+                (gm_tfm, score_and_scrub_cbf, [("output_image", "gm_tpm")]),
+                (wm_tfm, score_and_scrub_cbf, [("output_image", "wm_tpm")]),
+                (csf_tfm, score_and_scrub_cbf, [("output_image", "csf_tpm")]),
+                (score_and_scrub_cbf, outputnode, [
+                    ("cbf_ts_score", "cbf_ts_score"),
+                    ("score_outliers", "score_outliers"),
+                    ("mean_cbf_score", "mean_cbf_score"),
+                    ("mean_cbf_scrub", "mean_cbf_scrub"),
+                ]),
+            ])
+            # fmt:on
+
     return workflow
 
 
@@ -416,11 +418,11 @@ def init_compute_cbf_ge_wf(
             fields=[
                 "out_cbf",
                 "out_mean",
-                "out_score",
-                "out_avgscore",
-                "out_scrub",
+                "cbf_ts_score",
+                "score_outliers",
+                "mean_cbf_score",
+                "mean_cbf_scrub",
                 "out_cbfb",
-                "out_scoreindex",
                 "out_cbfpv",
                 "out_att",
                 "out_cbfpvwm",
@@ -621,7 +623,7 @@ def _getfiledir(file):
 """
 
         score_and_scrub_cbf = pe.Node(
-            ScoreAndScrubCBF(in_thresh=0.7, in_wfun="huber"),
+            ScoreAndScrubCBF(tpm_threshold=0.7, wavelet_function="huber"),
             mem_gb=mem_gb,
             name="scorescrub",
             run_without_submitting=True,
@@ -629,16 +631,16 @@ def _getfiledir(file):
 
         # fmt:off
         workflow.connect([
-            (refine_mask, score_and_scrub_cbf, [("out_mask", "in_mask")]),
-            (gm_tfm, score_and_scrub_cbf, [("output_image", "in_greyM")]),
-            (wm_tfm, score_and_scrub_cbf, [("output_image", "in_whiteM")]),
-            (csf_tfm, score_and_scrub_cbf, [("output_image", "in_csf")]),
-            (collect_cbf, score_and_scrub_cbf, [("cbf", "in_file")]),
+            (refine_mask, score_and_scrub_cbf, [("out_mask", "brain_mask")]),
+            (gm_tfm, score_and_scrub_cbf, [("output_image", "gm_tpm")]),
+            (wm_tfm, score_and_scrub_cbf, [("output_image", "wm_tpm")]),
+            (csf_tfm, score_and_scrub_cbf, [("output_image", "csf_tpm")]),
+            (collect_cbf, score_and_scrub_cbf, [("cbf", "cbf")]),
             (score_and_scrub_cbf, outputnode, [
-                ("out_score", "out_score"),
-                ("out_scoreindex", "out_scoreindex"),
-                ("out_avgscore", "out_avgscore"),
-                ("out_scrub", "out_scrub"),
+                ("cbf_ts_score", "cbf_ts_score"),
+                ("score_outliers", "score_outliers"),
+                ("mean_cbf_score", "mean_cbf_score"),
+                ("mean_cbf_scrub", "mean_cbf_scrub"),
             ]),
         ])
         # fmt:on

From 535814a8d04898015d0b98c4beecce03403e0a36 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Tue, 4 Apr 2023 12:19:34 -0400
Subject: [PATCH 06/27] Patch connections.

---
 aslprep/workflows/asl/base.py  | 22 +++++++++++-----------
 aslprep/workflows/asl/gecbf.py | 18 +++++++++---------
 2 files changed, 20 insertions(+), 20 deletions(-)

diff --git a/aslprep/workflows/asl/base.py b/aslprep/workflows/asl/base.py
index 14634e02c..c0989883c 100644
--- a/aslprep/workflows/asl/base.py
+++ b/aslprep/workflows/asl/base.py
@@ -739,9 +739,9 @@ def init_asl_preproc_wf(asl_file):
             # fmt:off
             workflow.connect([
                 (compute_cbf_wf, asl_derivatives_wf, [
-                    ("outputnode.out_score", "inputnode.score"),
-                    ("outputnode.out_avgscore", "inputnode.avgscore"),
-                    ("outputnode.out_scrub", "inputnode.scrub"),
+                    ("outputnode.cbf_ts_score", "inputnode.score"),
+                    ("outputnode.mean_cbf_score", "inputnode.avgscore"),
+                    ("outputnode.mean_cbf_scrub", "inputnode.scrub"),
                 ]),
             ])
             # fmt:on
@@ -795,9 +795,9 @@ def init_asl_preproc_wf(asl_file):
             # fmt:off
             workflow.connect([
                 (compute_cbf_wf, asl_std_trans_wf, [
-                    ("outputnode.out_score", "inputnode.score"),
-                    ("outputnode.out_avgscore", "inputnode.avgscore"),
-                    ("outputnode.out_scrub", "inputnode.scrub"),
+                    ("outputnode.cbf_ts_score", "inputnode.score"),
+                    ("outputnode.mean_cbf_score", "inputnode.avgscore"),
+                    ("outputnode.mean_cbf_scrub", "inputnode.scrub"),
                 ]),
             ])
             # fmt:on
@@ -927,14 +927,14 @@ def init_asl_preproc_wf(asl_file):
     # fmt:off
     workflow.connect([
         (compute_cbf_wf, cbf_plot, [
+            ("outputnode.out_cbf", "inputnode.cbf_ts"),
             ("outputnode.out_mean", "inputnode.cbf"),
-            ("outputnode.out_avgscore", "inputnode.score"),
-            ("outputnode.out_scrub", "inputnode.scrub"),
             ("outputnode.out_cbfb", "inputnode.basil"),
             ("outputnode.out_cbfpv", "inputnode.pvc"),
-            ("outputnode.out_score", "inputnode.score_ts"),
-            ("outputnode.out_cbf", "inputnode.cbf_ts"),
-            ("outputnode.out_scoreindex", "inputnode.scoreindex"),
+            ("outputnode.cbf_ts_score", "inputnode.score_ts"),
+            ("outputnode.mean_cbf_score", "inputnode.score"),
+            ("outputnode.mean_cbf_scrub", "inputnode.scrub"),
+            ("outputnode.score_outliers", "inputnode.scoreindex"),
         ]),
         (inputnode, cbf_plot, [("std2anat_xfm", "inputnode.std2anat_xfm")]),
         (asl_reg_wf, cbf_plot, [("outputnode.itk_t1_to_asl", "inputnode.t1_asl_xform")]),
diff --git a/aslprep/workflows/asl/gecbf.py b/aslprep/workflows/asl/gecbf.py
index 840d0916b..187e55086 100644
--- a/aslprep/workflows/asl/gecbf.py
+++ b/aslprep/workflows/asl/gecbf.py
@@ -497,9 +497,9 @@ def init_asl_gepreproc_wf(asl_file):
             # fmt:off
             workflow.connect([
                 (compute_cbf_wf, t1w_gereg_wf, [
-                    ("outputnode.out_score", "inputnode.score"),
-                    ("outputnode.out_avgscore", "inputnode.avgscore"),
-                    ("outputnode.out_scrub", "inputnode.scrub"),
+                    ("outputnode.cbf_ts_score", "inputnode.score"),
+                    ("outputnode.mean_cbf_score", "inputnode.avgscore"),
+                    ("outputnode.mean_cbf_scrub", "inputnode.scrub"),
                 ]),
                 (t1w_gereg_wf, asl_derivatives_wf, [
                     ("outputnode.score_t1", "inputnode.score_t1"),
@@ -588,9 +588,9 @@ def init_asl_gepreproc_wf(asl_file):
             # fmt:off
             workflow.connect([
                 (compute_cbf_wf, asl_derivatives_wf, [
-                    ("outputnode.out_score", "inputnode.score"),
-                    ("outputnode.out_avgscore", "inputnode.avgscore"),
-                    ("outputnode.out_scrub", "inputnode.scrub"),
+                    ("outputnode.cbf_ts_score", "inputnode.score"),
+                    ("outputnode.mean_cbf_score", "inputnode.avgscore"),
+                    ("outputnode.mean_cbf_scrub", "inputnode.scrub"),
                 ]),
             ])
             # fmt:on
@@ -644,9 +644,9 @@ def init_asl_gepreproc_wf(asl_file):
             # fmt:off
             workflow.connect([
                 (compute_cbf_wf, std_gereg_wf, [
-                    ("outputnode.out_score", "inputnode.score"),
-                    ("outputnode.out_avgscore", "inputnode.avgscore"),
-                    ("outputnode.out_scrub", "inputnode.scrub"),
+                    ("outputnode.cbf_ts_score", "inputnode.score"),
+                    ("outputnode.mean_cbf_score", "inputnode.avgscore"),
+                    ("outputnode.mean_cbf_scrub", "inputnode.scrub"),
                 ]),
             ])
             # fmt:on

From 0cd690a0b3e9688c08424c630ad25b8736584b92 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Tue, 4 Apr 2023 12:41:44 -0400
Subject: [PATCH 07/27] Patch more connections.

---
 aslprep/workflows/asl/base.py  | 40 +++++++++++++++++-----------------
 aslprep/workflows/asl/gecbf.py | 30 ++++++++++++-------------
 2 files changed, 35 insertions(+), 35 deletions(-)

diff --git a/aslprep/workflows/asl/base.py b/aslprep/workflows/asl/base.py
index c0989883c..bfe96ffcb 100644
--- a/aslprep/workflows/asl/base.py
+++ b/aslprep/workflows/asl/base.py
@@ -674,8 +674,8 @@ def init_asl_preproc_wf(asl_file):
             (refine_mask, aslmask_to_t1w, [("out_mask", "input_image")]),
             (aslmask_to_t1w, outputnode, [("output_image", "asl_mask_t1")]),
             (compute_cbf_wf, asl_t1_trans_wf, [
-                ("outputnode.out_cbf", "inputnode.cbf"),
-                ("outputnode.out_mean", "inputnode.meancbf"),
+                ("outputnode.cbf_ts", "inputnode.cbf"),
+                ("outputnode.mean_cbf", "inputnode.meancbf"),
             ]),
             (asl_t1_trans_wf, asl_derivatives_wf, [
                 ("outputnode.cbf_t1", "inputnode.cbf_t1"),
@@ -688,14 +688,14 @@ def init_asl_preproc_wf(asl_file):
             # fmt:off
             workflow.connect([
                 (compute_cbf_wf, asl_t1_trans_wf, [
-                    ("outputnode.cbf_ts_score", "inputnode.cbf_ts_score"),
-                    ("outputnode.mean_cbf_score", "inputnode.mean_cbf_score"),
-                    ("outputnode.mean_cbf_scrub", "inputnode.mean_cbf_scrub"),
+                    ("outputnode.cbf_ts_score", "inputnode.score"),
+                    ("outputnode.mean_cbf_score", "inputnode.avgscore"),
+                    ("outputnode.mean_cbf_scrub", "inputnode.scrub"),
                 ]),
                 (asl_t1_trans_wf, asl_derivatives_wf, [
-                    ("outputnode.cbf_ts_score_t1", "inputnode.cbf_ts_score_t1"),
-                    ("outputnode.mean_cbf_score_t1", "inputnode.mean_cbf_score_t1"),
-                    ("outputnode.mean_cbf_scrub_t1", "inputnode.mean_cbf_scrub_t1"),
+                    ("outputnode.score_t1", "inputnode.score_t1"),
+                    ("outputnode.avgscore_t1", "inputnode.avgscore_t1"),
+                    ("outputnode.scrub_t1", "inputnode.scrub_t1"),
                 ]),
             ])
             # fmt:on
@@ -729,8 +729,8 @@ def init_asl_preproc_wf(asl_file):
                 ("out_mask", "inputnode.asl_mask_native"),
             ]),
             (compute_cbf_wf, asl_derivatives_wf, [
-                ("outputnode.out_cbf", "inputnode.cbf"),
-                ("outputnode.out_mean", "inputnode.meancbf"),
+                ("outputnode.cbf_ts", "inputnode.cbf"),
+                ("outputnode.mean_cbf", "inputnode.meancbf"),
             ]),
         ])
         # fmt:on
@@ -785,8 +785,8 @@ def init_asl_preproc_wf(asl_file):
             (refine_mask, asl_std_trans_wf, [("out_mask", "inputnode.asl_mask")]),
             (asl_sdc_wf, asl_std_trans_wf, [("outputnode.out_warp", "inputnode.fieldwarp")]),
             (compute_cbf_wf, asl_std_trans_wf, [
-                ("outputnode.out_cbf", "inputnode.cbf"),
-                ("outputnode.out_mean", "inputnode.meancbf"),
+                ("outputnode.cbf_ts", "inputnode.cbf"),
+                ("outputnode.mean_cbf", "inputnode.meancbf"),
             ]),
         ])
         # fmt:on
@@ -880,7 +880,7 @@ def init_asl_preproc_wf(asl_file):
         (inputnode, compt_qccbf_wf, [("t1w_tpms", "inputnode.t1w_tpms")]),
         (asl_reg_wf, compt_qccbf_wf, [("outputnode.itk_t1_to_asl", "inputnode.t1_asl_xform")]),
         (inputnode, compt_qccbf_wf, [("t1w_mask", "inputnode.t1w_mask")]),
-        (compute_cbf_wf, compt_qccbf_wf, [("outputnode.out_mean", "inputnode.meancbf")]),
+        (compute_cbf_wf, compt_qccbf_wf, [("outputnode.mean_cbf", "inputnode.meancbf")]),
         (asl_confounds_wf, compt_qccbf_wf, [("outputnode.confounds_file", "inputnode.confmat")]),
         (compt_qccbf_wf, outputnode, [("outputnode.qc_file", "qc_file")]),
         (compt_qccbf_wf, asl_derivatives_wf, [("outputnode.qc_file", "inputnode.qc_file")]),
@@ -893,8 +893,8 @@ def init_asl_preproc_wf(asl_file):
         # fmt:off
         workflow.connect([
             (compute_cbf_wf, compt_qccbf_wf, [
-                ("outputnode.out_avgscore", "inputnode.avgscore"),
-                ("outputnode.out_scrub", "inputnode.scrub"),
+                ("outputnode.mean_cbf_score", "inputnode.avgscore"),
+                ("outputnode.mean_cbf_scrub", "inputnode.scrub"),
             ]),
         ])
         # fmt:on
@@ -927,8 +927,8 @@ def init_asl_preproc_wf(asl_file):
     # fmt:off
     workflow.connect([
         (compute_cbf_wf, cbf_plot, [
-            ("outputnode.out_cbf", "inputnode.cbf_ts"),
-            ("outputnode.out_mean", "inputnode.cbf"),
+            ("outputnode.cbf_ts", "inputnode.cbf_ts"),
+            ("outputnode.mean_cbf", "inputnode.cbf"),
             ("outputnode.out_cbfb", "inputnode.basil"),
             ("outputnode.out_cbfpv", "inputnode.pvc"),
             ("outputnode.cbf_ts_score", "inputnode.score_ts"),
@@ -987,7 +987,7 @@ def init_asl_preproc_wf(asl_file):
         (asl_asl_trans_wf, cbfroiqu, [("outputnode.asl_mask", "inputnode.aslmask")]),
         (inputnode, cbfroiqu, [("std2anat_xfm", "inputnode.std2anat_xfm")]),
         (asl_reg_wf, cbfroiqu, [("outputnode.itk_t1_to_asl", "inputnode.t1_asl_xform")]),
-        (compute_cbf_wf, cbfroiqu, [("outputnode.out_mean", "inputnode.cbf")]),
+        (compute_cbf_wf, cbfroiqu, [("outputnode.mean_cbf", "inputnode.cbf")]),
         (cbfroiqu, asl_derivatives_wf, [
             ("outputnode.cbf_hvoxf", "inputnode.cbf_hvoxf"),
             ("outputnode.cbf_sc207", "inputnode.cbf_sc207"),
@@ -1002,8 +1002,8 @@ def init_asl_preproc_wf(asl_file):
         # fmt:off
         workflow.connect([
             (compute_cbf_wf, cbfroiqu, [
-                ("outputnode.out_avgscore", "inputnode.score"),
-                ("outputnode.out_scrub", "inputnode.scrub"),
+                ("outputnode.mean_cbf_score", "inputnode.score"),
+                ("outputnode.mean_cbf_scrub", "inputnode.scrub"),
             ]),
             (cbfroiqu, asl_derivatives_wf, [
                 ("outputnode.score_hvoxf", "inputnode.score_hvoxf"),
diff --git a/aslprep/workflows/asl/gecbf.py b/aslprep/workflows/asl/gecbf.py
index 187e55086..a661ee653 100644
--- a/aslprep/workflows/asl/gecbf.py
+++ b/aslprep/workflows/asl/gecbf.py
@@ -415,8 +415,8 @@ def init_asl_gepreproc_wf(asl_file):
             ("outputnode.asl_t1_ref", "asl_t1_ref"),
         ]),
         (compute_cbf_wf, t1w_gereg_wf, [
-            ("outputnode.out_cbf", "inputnode.cbf"),
-            ("outputnode.out_mean", "inputnode.meancbf"),
+            ("outputnode.cbf_ts", "inputnode.cbf"),
+            ("outputnode.mean_cbf", "inputnode.meancbf"),
         ]),
         (t1w_gereg_wf, asl_derivatives_wf, [
             ("outputnode.cbf_t1", "inputnode.cbf_t1"),
@@ -448,7 +448,7 @@ def init_asl_gepreproc_wf(asl_file):
 
     # fmt:off
     workflow.connect([
-        (compute_cbf_wf, cbf_plot, [("outputnode.out_mean", "inputnode.cbf")]),
+        (compute_cbf_wf, cbf_plot, [("outputnode.mean_cbf", "inputnode.cbf")]),
         (gen_ref_wf, cbf_plot, [("outputnode.ref_image_brain", "inputnode.asl_ref")]),
     ])
     # fmt:on
@@ -457,8 +457,8 @@ def init_asl_gepreproc_wf(asl_file):
         # fmt:off
         workflow.connect([
             (compute_cbf_wf, cbf_plot, [
-                ("outputnode.out_avgscore", "inputnode.score"),
-                ("outputnode.out_scrub", "inputnode.scrub"),
+                ("outputnode.mean_cbf_score", "inputnode.score"),
+                ("outputnode.mean_cbf_scrub", "inputnode.scrub"),
             ]),
         ])
         # fmt:on
@@ -542,7 +542,7 @@ def init_asl_gepreproc_wf(asl_file):
         ]),
         (refine_mask, compt_qccbf_wf, [("out_mask", "inputnode.asl_mask")]),
         (asl_reg_wf, compt_qccbf_wf, [("outputnode.itk_t1_to_asl", "inputnode.t1_asl_xform")]),
-        (compute_cbf_wf, compt_qccbf_wf, [("outputnode.out_mean", "inputnode.meancbf")]),
+        (compute_cbf_wf, compt_qccbf_wf, [("outputnode.mean_cbf", "inputnode.meancbf")]),
         (compt_qccbf_wf, outputnode, [("outputnode.qc_file", "qc_file")]),
         (compt_qccbf_wf, asl_derivatives_wf, [("outputnode.qc_file", "inputnode.qc_file")]),
         (compt_qccbf_wf, summary, [("outputnode.qc_file", "qc_file")]),
@@ -553,8 +553,8 @@ def init_asl_gepreproc_wf(asl_file):
         # fmt:off
         workflow.connect([
             (compute_cbf_wf, compt_qccbf_wf, [
-                ("outputnode.out_avgscore", "inputnode.avgscore"),
-                ("outputnode.out_scrub", "inputnode.scrub"),
+                ("outputnode.mean_cbf_score", "inputnode.avgscore"),
+                ("outputnode.mean_cbf_scrub", "inputnode.scrub"),
             ]),
         ])
         # fmt:on
@@ -578,8 +578,8 @@ def init_asl_gepreproc_wf(asl_file):
             ]),
             (refine_mask, asl_derivatives_wf, [("out_mask", "inputnode.asl_mask_native")]),
             (compute_cbf_wf, asl_derivatives_wf, [
-                ("outputnode.out_cbf", "inputnode.cbf"),
-                ("outputnode.out_mean", "inputnode.meancbf"),
+                ("outputnode.cbf_ts", "inputnode.cbf"),
+                ("outputnode.mean_cbf", "inputnode.meancbf"),
             ]),
         ])
         # fmt:on
@@ -634,8 +634,8 @@ def init_asl_gepreproc_wf(asl_file):
                 ("outputnode.asl_mask_std", "asl_mask_std"),
             ]),
             (compute_cbf_wf, std_gereg_wf, [
-                ("outputnode.out_cbf", "inputnode.cbf"),
-                ("outputnode.out_mean", "inputnode.meancbf"),
+                ("outputnode.cbf_ts", "inputnode.cbf"),
+                ("outputnode.mean_cbf", "inputnode.meancbf"),
             ]),
         ])
         # fmt:on
@@ -721,7 +721,7 @@ def init_asl_gepreproc_wf(asl_file):
         (refine_mask, cbfroiqu, [("out_mask", "inputnode.aslmask")]),
         (inputnode, cbfroiqu, [("std2anat_xfm", "inputnode.std2anat_xfm")]),
         (asl_reg_wf, cbfroiqu, [("outputnode.itk_t1_to_asl", "inputnode.t1_asl_xform")]),
-        (compute_cbf_wf, cbfroiqu, [("outputnode.out_mean", "inputnode.cbf")]),
+        (compute_cbf_wf, cbfroiqu, [("outputnode.mean_cbf", "inputnode.cbf")]),
         (cbfroiqu, asl_derivatives_wf, [
             ("outputnode.cbf_hvoxf", "inputnode.cbf_hvoxf"),
             ("outputnode.cbf_sc207", "inputnode.cbf_sc207"),
@@ -736,8 +736,8 @@ def init_asl_gepreproc_wf(asl_file):
         # fmt:off
         workflow.connect([
             (compute_cbf_wf, cbfroiqu, [
-                ("outputnode.out_avgscore", "inputnode.score"),
-                ("outputnode.out_scrub", "inputnode.scrub"),
+                ("outputnode.mean_cbf_score", "inputnode.score"),
+                ("outputnode.mean_cbf_scrub", "inputnode.scrub"),
             ]),
             (cbfroiqu, asl_derivatives_wf, [
                 ("outputnode.score_hvoxf", "inputnode.score_hvoxf"),

From 292b0c02d43609b1a96850dc18e76fee5994817a Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Tue, 4 Apr 2023 13:26:07 -0400
Subject: [PATCH 08/27] Rename and document BASIL inputs and outputs.

---
 aslprep/interfaces/cbf_computation.py |  32 +++--
 aslprep/workflows/asl/base.py         |  36 +++---
 aslprep/workflows/asl/cbf.py          | 171 +++++++++++++-------------
 aslprep/workflows/asl/gecbf.py        |  36 +++---
 4 files changed, 145 insertions(+), 130 deletions(-)

diff --git a/aslprep/interfaces/cbf_computation.py b/aslprep/interfaces/cbf_computation.py
index 389d3f696..d3bcb3a37 100644
--- a/aslprep/interfaces/cbf_computation.py
+++ b/aslprep/interfaces/cbf_computation.py
@@ -797,6 +797,7 @@ class _BASILCBFInputSpec(FSLCommandInputSpec):
         mandatory=False,
         argstr="--pvcorr",
         default_value=True,
+        usedefault=True,
     )
     pvgm = File(
         exists=True,
@@ -821,12 +822,22 @@ class _BASILCBFInputSpec(FSLCommandInputSpec):
 
 
 class _BASILCBFOutputSpec(TraitedSpec):
-    out_cbfb = File(exists=False, desc="cbf with spatial correction")
-    out_cbfpv = File(exists=False, desc="cbf with spatial correction")
-    out_cbfpvwm = File(
-        exists=False, desc="cbf with spatial partial volume white matter correction"
+    mean_cbf = File(exists=True, desc="CBF map in absolute units (ml/100g/min).")
+    mean_cbf_gm = File(
+        exists=True,
+        desc=(
+            "CBF map with gray matter partial volume correction. "
+            "This means that the map contains 'pure' GM perfusion estimates, in absolute units."
+        ),
+    )
+    mean_cbf_wm = File(
+        exists=True,
+        desc=(
+            "CBF map with white matter partial volume correction. "
+            "This means that the map contains 'pure' WM perfusion estimates, in absolute units."
+        ),
     )
-    out_att = File(exists=False, desc="aretrial transist time")
+    att = File(exists=True, desc="Arterial transit time map.")
 
 
 class BASILCBF(FSLCommand):
@@ -859,10 +870,13 @@ def _list_outputs(self):
 
         outputs = self.output_spec().get()
 
-        outputs["out_cbfb"] = os.path.join(basename, "native_space/perfusion_calib.nii.gz")
-        outputs["out_att"] = os.path.join(basename, "native_space/arrival.nii.gz")
-        outputs["out_cbfpv"] = os.path.join(basename, "native_space/pvcorr/perfusion_calib.nii.gz")
-        outputs["out_cbfpvwm"] = os.path.join(
+        outputs["mean_cbf"] = os.path.join(basename, "native_space/perfusion_calib.nii.gz")
+        outputs["att"] = os.path.join(basename, "native_space/arrival.nii.gz")
+        outputs["mean_cbf_gm"] = os.path.join(
+            basename,
+            "native_space/pvcorr/perfusion_calib.nii.gz",
+        )
+        outputs["mean_cbf_wm"] = os.path.join(
             basename,
             "native_space/pvcorr/perfusion_wm_calib.nii.gz",
         )
diff --git a/aslprep/workflows/asl/base.py b/aslprep/workflows/asl/base.py
index bfe96ffcb..ae8dadde5 100644
--- a/aslprep/workflows/asl/base.py
+++ b/aslprep/workflows/asl/base.py
@@ -704,10 +704,10 @@ def init_asl_preproc_wf(asl_file):
             # fmt:off
             workflow.connect([
                 (compute_cbf_wf, asl_t1_trans_wf, [
-                    ("outputnode.out_cbfb", "inputnode.basil"),
-                    ("outputnode.out_cbfpv", "inputnode.pv"),
-                    ("outputnode.out_cbfpvwm", "inputnode.pvwm"),
-                    ("outputnode.out_att", "inputnode.att"),
+                    ("outputnode.mean_cbf_basil", "inputnode.basil"),
+                    ("outputnode.mean_cbf_gm_basil", "inputnode.pv"),
+                    ("outputnode.mean_cbf_wm_basil", "inputnode.pvwm"),
+                    ("outputnode.att_basil", "inputnode.att"),
                 ]),
                 (asl_t1_trans_wf, asl_derivatives_wf, [
                     ("outputnode.basil_t1", "inputnode.basil_t1"),
@@ -750,10 +750,10 @@ def init_asl_preproc_wf(asl_file):
             # fmt:off
             workflow.connect([
                 (compute_cbf_wf, asl_derivatives_wf, [
-                    ("outputnode.out_cbfb", "inputnode.basil"),
-                    ("outputnode.out_cbfpv", "inputnode.pv"),
-                    ("outputnode.out_cbfpvwm", "inputnode.pvwm"),
-                    ("outputnode.out_att", "inputnode.att"),
+                    ("outputnode.mean_cbf_basil", "inputnode.basil"),
+                    ("outputnode.mean_cbf_gm_basil", "inputnode.pv"),
+                    ("outputnode.mean_cbf_wm_basil", "inputnode.pvwm"),
+                    ("outputnode.att_basil", "inputnode.att"),
                 ]),
             ])
             # fmt:on
@@ -806,10 +806,10 @@ def init_asl_preproc_wf(asl_file):
             # fmt:off
             workflow.connect([
                 (compute_cbf_wf, asl_std_trans_wf, [
-                    ("outputnode.out_cbfb", "inputnode.basil"),
-                    ("outputnode.out_cbfpv", "inputnode.pv"),
-                    ("outputnode.out_cbfpvwm", "inputnode.pvwm"),
-                    ("outputnode.out_att", "inputnode.att"),
+                    ("outputnode.mean_cbf_basil", "inputnode.basil"),
+                    ("outputnode.mean_cbf_gm_basil", "inputnode.pv"),
+                    ("outputnode.mean_cbf_wm_basil", "inputnode.pvwm"),
+                    ("outputnode.att_basil", "inputnode.att"),
                 ]),
             ])
             # fmt:on
@@ -903,8 +903,8 @@ def init_asl_preproc_wf(asl_file):
         # fmt:off
         workflow.connect([
             (compute_cbf_wf, compt_qccbf_wf, [
-                ("outputnode.out_cbfb", "inputnode.basil"),
-                ("outputnode.out_cbfpv", "inputnode.pv"),
+                ("outputnode.mean_cbf_basil", "inputnode.basil"),
+                ("outputnode.mean_cbf_gm_basil", "inputnode.pv"),
             ]),
         ])
         # fmt:on
@@ -929,12 +929,12 @@ def init_asl_preproc_wf(asl_file):
         (compute_cbf_wf, cbf_plot, [
             ("outputnode.cbf_ts", "inputnode.cbf_ts"),
             ("outputnode.mean_cbf", "inputnode.cbf"),
-            ("outputnode.out_cbfb", "inputnode.basil"),
-            ("outputnode.out_cbfpv", "inputnode.pvc"),
             ("outputnode.cbf_ts_score", "inputnode.score_ts"),
             ("outputnode.mean_cbf_score", "inputnode.score"),
             ("outputnode.mean_cbf_scrub", "inputnode.scrub"),
             ("outputnode.score_outliers", "inputnode.scoreindex"),
+            ("outputnode.mean_cbf_basil", "inputnode.basil"),
+            ("outputnode.mean_cbf_gm_basil", "inputnode.pvc"),
         ]),
         (inputnode, cbf_plot, [("std2anat_xfm", "inputnode.std2anat_xfm")]),
         (asl_reg_wf, cbf_plot, [("outputnode.itk_t1_to_asl", "inputnode.t1_asl_xform")]),
@@ -1024,8 +1024,8 @@ def init_asl_preproc_wf(asl_file):
         # fmt:off
         workflow.connect([
             (compute_cbf_wf, cbfroiqu, [
-                ("outputnode.out_cbfb", "inputnode.basil"),
-                ("outputnode.out_cbfpv", "inputnode.pvc"),
+                ("outputnode.mean_cbf_basil", "inputnode.basil"),
+                ("outputnode.mean_cbf_gm_basil", "inputnode.pvc"),
             ]),
             (cbfroiqu, asl_derivatives_wf, [
                 ("outputnode.basil_hvoxf", "inputnode.basil_hvoxf"),
diff --git a/aslprep/workflows/asl/cbf.py b/aslprep/workflows/asl/cbf.py
index ce099527c..b36adc7ad 100644
--- a/aslprep/workflows/asl/cbf.py
+++ b/aslprep/workflows/asl/cbf.py
@@ -101,6 +101,8 @@ def init_compute_cbf_wf(
 
 *ASLPrep* was configured to calculate cerebral blood flow (CBF) using the following methods.
 
+The cerebral blood flow (CBF) was quantified from preprocessed ASL data using a general kinetic
+model [@buxton1998general].
 """
     inputnode = pe.Node(
         niu.IdentityInterface(
@@ -123,14 +125,16 @@ def init_compute_cbf_wf(
             fields=[
                 "cbf_ts",
                 "mean_cbf",
+                # SCORE/SCRUB outputs
                 "cbf_ts_score",
                 "score_outliers",
                 "mean_cbf_score",
                 "mean_cbf_scrub",
-                "out_cbfpvwm",
-                "out_cbfb",
-                "out_cbfpv",
-                "out_att",
+                # BASIL outputs
+                "mean_cbf_basil",
+                "mean_cbf_gm_basil",
+                "mean_cbf_wm_basil",
+                "att_basil",
             ]
         ),
         name="outputnode",
@@ -254,6 +258,63 @@ def _getfiledir(file):
     ])
     # fmt:on
 
+    compute_cbf = pe.Node(
+        ComputeCBF(
+            cbf_only=cbf_only,
+            m0scale=M0Scale,
+        ),
+        mem_gb=0.2,
+        run_without_submitting=True,
+        name="compute_cbf",
+    )
+
+    # fmt:off
+    workflow.connect([
+        (refine_mask, compute_cbf, [("out_mask", "mask")]),
+        (extract_deltam, compute_cbf, [
+            ("out_file", "deltam"),
+            ("m0_file", "m0_file"),
+            ("metadata", "metadata"),
+        ]),
+        (compute_cbf, outputnode, [
+            ("cbf", "cbf_ts"),
+            ("mean_cbf", "mean_cbf"),
+        ]),
+    ])
+    # fmt:on
+
+    if scorescrub:
+        score_and_scrub_cbf = pe.Node(
+            ScoreAndScrubCBF(tpm_threshold=0.7, wavelet_function="huber"),
+            mem_gb=0.2,
+            name="score_and_scrub_cbf",
+            run_without_submitting=True,
+        )
+
+        workflow.__desc__ += """
+
+Structural Correlation based Outlier Rejection (SCORE) algorithm was applied to the CBF timeseries
+to discard CBF volumes with outlying values [@dolui2017structural] before computing the mean CBF.
+Following SCORE, the Structural Correlation with RobUst Bayesian (SCRUB) algorithm was applied to
+the CBF maps using structural tissue probability maps to reweight the mean CBF
+[@dolui2017structural;@dolui2016scrub].
+"""
+        # fmt:off
+        workflow.connect([
+            (refine_mask, score_and_scrub_cbf, [("out_mask", "brain_mask")]),
+            (compute_cbf, score_and_scrub_cbf, [("cbf", "cbf")]),
+            (gm_tfm, score_and_scrub_cbf, [("output_image", "gm_tpm")]),
+            (wm_tfm, score_and_scrub_cbf, [("output_image", "wm_tpm")]),
+            (csf_tfm, score_and_scrub_cbf, [("output_image", "csf_tpm")]),
+            (score_and_scrub_cbf, outputnode, [
+                ("cbf_ts_score", "cbf_ts_score"),
+                ("score_outliers", "score_outliers"),
+                ("mean_cbf_score", "mean_cbf_score"),
+                ("mean_cbf_scrub", "mean_cbf_scrub"),
+            ]),
+        ])
+        # fmt:on
+
     if basil:
         workflow.__desc__ += f"""
 CBF was computed with Bayesian Inference for Arterial Spin Labeling (BASIL)
@@ -295,76 +356,14 @@ def _getfiledir(file):
             (gm_tfm, basilcbf, [("output_image", "pvgm")]),
             (wm_tfm, basilcbf, [("output_image", "pvwm")]),
             (basilcbf, outputnode, [
-                ("out_cbfb", "cbf"),
-                ("out_cbfpv", "cbf_pv"),
-                ("out_cbfpvwm", "cbf_pv_wm"),
-                ("out_att", "att"),
-            ]),
-        ])
-        # fmt:on
-
-    else:
-        workflow.__desc__ += """
-The cerebral blood flow (CBF) was quantified from preprocessed ASL data using a general kinetic
-model [@buxton1998general].
-"""
-        compute_cbf = pe.Node(
-            ComputeCBF(
-                cbf_only=cbf_only,
-                m0scale=M0Scale,
-            ),
-            mem_gb=0.2,
-            run_without_submitting=True,
-            name="compute_cbf",
-        )
-
-        # fmt:off
-        workflow.connect([
-            (refine_mask, compute_cbf, [("out_mask", "mask")]),
-            (extract_deltam, compute_cbf, [
-                ("out_file", "deltam"),
-                ("m0_file", "m0_file"),
-                ("metadata", "metadata"),
-            ]),
-            (compute_cbf, outputnode, [
-                ("cbf", "out_cbf"),
-                ("mean_cbf", "out_mean"),
+                ("mean_cbf", "mean_cbf_basil"),
+                ("mean_cbf_gm", "mean_cbf_gm_basil"),
+                ("mean_cbf_wm", "mean_cbf_wm_basil"),
+                ("att", "att_basil"),
             ]),
         ])
         # fmt:on
 
-        if scorescrub:
-            score_and_scrub_cbf = pe.Node(
-                ScoreAndScrubCBF(tpm_threshold=0.7, wavelet_function="huber"),
-                mem_gb=0.2,
-                name="score_and_scrub_cbf",
-                run_without_submitting=True,
-            )
-
-            workflow.__desc__ += """
-
-Structural Correlation based Outlier Rejection (SCORE) algorithm was applied to the CBF timeseries
-to discard CBF volumes with outlying values [@dolui2017structural] before computing the mean CBF.
-Following SCORE, the Structural Correlation with RobUst Bayesian (SCRUB) algorithm was applied to
-the CBF maps using structural tissue probability maps to reweight the mean CBF
-[@dolui2017structural;@dolui2016scrub].
-"""
-            # fmt:off
-            workflow.connect([
-                (refine_mask, score_and_scrub_cbf, [("out_mask", "brain_mask")]),
-                (compute_cbf, score_and_scrub_cbf, [("cbf", "cbf")]),
-                (gm_tfm, score_and_scrub_cbf, [("output_image", "gm_tpm")]),
-                (wm_tfm, score_and_scrub_cbf, [("output_image", "wm_tpm")]),
-                (csf_tfm, score_and_scrub_cbf, [("output_image", "csf_tpm")]),
-                (score_and_scrub_cbf, outputnode, [
-                    ("cbf_ts_score", "cbf_ts_score"),
-                    ("score_outliers", "score_outliers"),
-                    ("mean_cbf_score", "mean_cbf_score"),
-                    ("mean_cbf_scrub", "mean_cbf_scrub"),
-                ]),
-            ])
-            # fmt:on
-
     return workflow
 
 
@@ -416,16 +415,18 @@ def init_compute_cbf_ge_wf(
     outputnode = pe.Node(
         niu.IdentityInterface(
             fields=[
-                "out_cbf",
-                "out_mean",
+                "cbf_ts",
+                "mean_cbf",
+                # SCORE/SCRUB outputs
                 "cbf_ts_score",
                 "score_outliers",
                 "mean_cbf_score",
                 "mean_cbf_scrub",
-                "out_cbfb",
-                "out_cbfpv",
-                "out_att",
-                "out_cbfpvwm",
+                # BASIL outputs
+                "mean_cbf_basil",
+                "mean_cbf_gm_basil",
+                "mean_cbf_wm_basil",
+                "att_basil",
             ]
         ),
         name="outputnode",
@@ -572,8 +573,8 @@ def _getfiledir(file):
             (refine_mask, compute_cbf, [("out_mask", "mask")]),
             (compute_cbf, collect_cbf, [("cbf", "cbf")]),
             (compute_cbf, outputnode, [
-                ("cbf", "out_cbf"),
-                ("mean_cbf", "out_mean"),
+                ("cbf", "cbf_ts"),
+                ("mean_cbf", "mean_cbf"),
             ]),
         ])
         # fmt:on
@@ -608,8 +609,8 @@ def _getfiledir(file):
             (extract_cbf, mask_cbf, [("out_file", "operand_files")]),
             (mask_cbf, collect_cbf, [("out_file", "cbf")]),
             (mask_cbf, outputnode, [
-                ("out_file", "out_cbf"),
-                ("out_file", "out_mean"),
+                ("out_file", "cbf_ts"),
+                ("out_file", "mean_cbf"),
             ]),
         ])
         # fmt:on
@@ -646,7 +647,7 @@ def _getfiledir(file):
         # fmt:on
 
     if basil:
-        workflow.__desc__ += """\
+        workflow.__desc__ += """
 In addition, CBF was also computed by Bayesian Inference for Arterial Spin Labeling
 (BASIL) as implemented in FSL. BASIL is based on Bayesian inference principles
 [@chappell2008variational], and computed CBF from ASL by incorporating natural
@@ -689,10 +690,10 @@ def _getfiledir(file):
             (wm_tfm, basilcbf, [("output_image", "pvwm")]),
             (refine_mask, basilcbf, [("out_mask", "mask")]),
             (basilcbf, outputnode, [
-                ("out_cbfb", "out_cbfb"),
-                ("out_cbfpv", "out_cbfpv"),
-                ("out_cbfpvwm", "out_cbfpvwm"),
-                ("out_att", "out_att"),
+                ("mean_cbf", "mean_cbf_basil"),
+                ("mean_cbf_gm", "mean_cbf_gm_basil"),
+                ("mean_cbf_wm", "mean_cbf_wm_basil"),
+                ("att", "att_basil"),
             ]),
         ])
         # fmt:on
diff --git a/aslprep/workflows/asl/gecbf.py b/aslprep/workflows/asl/gecbf.py
index a661ee653..4417b14a1 100644
--- a/aslprep/workflows/asl/gecbf.py
+++ b/aslprep/workflows/asl/gecbf.py
@@ -467,8 +467,8 @@ def init_asl_gepreproc_wf(asl_file):
         # fmt:off
         workflow.connect([
             (compute_cbf_wf, cbf_plot, [
-                ("outputnode.out_cbfb", "inputnode.basil"),
-                ("outputnode.out_cbfpv", "inputnode.pvc"),
+                ("outputnode.mean_cbf_basil", "inputnode.basil"),
+                ("outputnode.mean_cbf_gm_basil", "inputnode.pvc"),
             ]),
         ])
         # fmt:on
@@ -513,10 +513,10 @@ def init_asl_gepreproc_wf(asl_file):
             # fmt:off
             workflow.connect([
                 (compute_cbf_wf, t1w_gereg_wf, [
-                    ("outputnode.out_cbfb", "inputnode.basil"),
-                    ("outputnode.out_cbfpv", "inputnode.pv"),
-                    ("outputnode.out_cbfpvwm", "inputnode.pvwm"),
-                    ("outputnode.out_att", "inputnode.att"),
+                    ("outputnode.mean_cbf_basil", "inputnode.basil"),
+                    ("outputnode.mean_cbf_gm_basil", "inputnode.pv"),
+                    ("outputnode.mean_cbf_wm_basil", "inputnode.pvwm"),
+                    ("outputnode.att_basil", "inputnode.att"),
                 ]),
                 (t1w_gereg_wf, asl_derivatives_wf, [
                     ("outputnode.basil_t1", "inputnode.basil_t1"),
@@ -563,8 +563,8 @@ def init_asl_gepreproc_wf(asl_file):
         # fmt:off
         workflow.connect([
             (compute_cbf_wf, compt_qccbf_wf, [
-                ("outputnode.out_cbfb", "inputnode.basil"),
-                ("outputnode.out_cbfpv", "inputnode.pv"),
+                ("outputnode.mean_cbf_basil", "inputnode.basil"),
+                ("outputnode.mean_cbf_gm_basil", "inputnode.pv"),
             ]),
         ])
         # fmt:on
@@ -599,10 +599,10 @@ def init_asl_gepreproc_wf(asl_file):
             # fmt:off
             workflow.connect([
                 (compute_cbf_wf, asl_derivatives_wf, [
-                    ("outputnode.out_cbfb", "inputnode.basil"),
-                    ("outputnode.out_cbfpv", "inputnode.pv"),
-                    ("outputnode.out_cbfpvwm", "inputnode.pvwm"),
-                    ("outputnode.out_att", "inputnode.att"),
+                    ("outputnode.mean_cbf_basil", "inputnode.basil"),
+                    ("outputnode.mean_cbf_gm_basil", "inputnode.pv"),
+                    ("outputnode.mean_cbf_wm_basil", "inputnode.pvwm"),
+                    ("outputnode.att_basil", "inputnode.att"),
                 ]),
             ])
             # fmt:on
@@ -655,10 +655,10 @@ def init_asl_gepreproc_wf(asl_file):
             # fmt:off
             workflow.connect([
                 (compute_cbf_wf, std_gereg_wf, [
-                    ("outputnode.out_cbfb", "inputnode.basil"),
-                    ("outputnode.out_cbfpv", "inputnode.pv"),
-                    ("outputnode.out_cbfpvwm", "inputnode.pvwm"),
-                    ("outputnode.out_att", "inputnode.att"),
+                    ("outputnode.mean_cbf_basil", "inputnode.basil"),
+                    ("outputnode.mean_cbf_gm_basil", "inputnode.pv"),
+                    ("outputnode.mean_cbf_wm_basil", "inputnode.pvwm"),
+                    ("outputnode.att_basil", "inputnode.att"),
                 ]),
             ])
             # fmt:on
@@ -758,8 +758,8 @@ def init_asl_gepreproc_wf(asl_file):
         # fmt:off
         workflow.connect([
             (compute_cbf_wf, cbfroiqu, [
-                ("outputnode.out_cbfb", "inputnode.basil"),
-                ("outputnode.out_cbfpv", "inputnode.pvc"),
+                ("outputnode.mean_cbf_basil", "inputnode.basil"),
+                ("outputnode.mean_cbf_gm_basil", "inputnode.pvc"),
             ]),
             (cbfroiqu, asl_derivatives_wf, [
                 ("outputnode.basil_hvoxf", "inputnode.basil_hvoxf"),

From 3e35cf1d1caeac5e904031bbb66d4ea83105dd07 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Tue, 4 Apr 2023 13:26:40 -0400
Subject: [PATCH 09/27] Drop mutual exclusion for this PR.

---
 aslprep/cli/parser.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/aslprep/cli/parser.py b/aslprep/cli/parser.py
index 3635916c8..3ba367d1b 100644
--- a/aslprep/cli/parser.py
+++ b/aslprep/cli/parser.py
@@ -329,8 +329,7 @@ def _bids_filter(value):
     )
 
     # CBF calculation settings
-    cbf_group = g_conf.add_mutually_exclusive_group()
-    cbf_group.add_argument(
+    g_conf.add_argument(
         "--scorescrub",
         action="store_true",
         default=False,
@@ -339,7 +338,7 @@ def _bids_filter(value):
             "These algorithms will be applied to the standard CBF results."
         ),
     )
-    cbf_group.add_argument(
+    g_conf.add_argument(
         "--basil",
         action="store_true",
         default=False,

From 4602ef00a189ce45c24015168a7787bfed93c5af Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Tue, 4 Apr 2023 13:28:24 -0400
Subject: [PATCH 10/27] Update parser.py

---
 aslprep/cli/parser.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/aslprep/cli/parser.py b/aslprep/cli/parser.py
index 3ba367d1b..4756ce019 100644
--- a/aslprep/cli/parser.py
+++ b/aslprep/cli/parser.py
@@ -343,8 +343,7 @@ def _bids_filter(value):
         action="store_true",
         default=False,
         help=(
-            "Use FSL's CBF computation with spatial regularization and partial volume correction, "
-            "rather than ASLPrep's standard CBF estimation method."
+            "Use FSL's CBF computation with spatial regularization and partial volume correction."
         ),
     )
 

From 28caafac92a6a80926c545680c7a3f816bbd3b94 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Tue, 4 Apr 2023 13:33:25 -0400
Subject: [PATCH 11/27] Split off QC interfaces.

---
 aslprep/interfaces/cbf_computation.py | 376 -------------------------
 aslprep/interfaces/qc.py              | 389 ++++++++++++++++++++++++++
 aslprep/workflows/asl/base.py         |   2 +-
 aslprep/workflows/asl/qc.py           |   2 +-
 4 files changed, 391 insertions(+), 378 deletions(-)
 create mode 100644 aslprep/interfaces/qc.py

diff --git a/aslprep/interfaces/cbf_computation.py b/aslprep/interfaces/cbf_computation.py
index d3bcb3a37..f6e859628 100644
--- a/aslprep/interfaces/cbf_computation.py
+++ b/aslprep/interfaces/cbf_computation.py
@@ -28,15 +28,6 @@
     parcellate_cbf,
     pcasl_or_pasl,
 )
-from aslprep.utils.qc import (
-    cbf_qei,
-    coverage,
-    crosscorr,
-    dice,
-    globalcbf,
-    jaccard,
-    negativevoxel,
-)
 
 
 class _RefineMaskInputSpec(BaseInterfaceInputSpec):
@@ -884,373 +875,6 @@ def _list_outputs(self):
         return outputs
 
 
-class _ComputeCBFQCInputSpec(BaseInterfaceInputSpec):
-    in_file = File(exists=True, mandatory=True, desc="original asl_file")
-    in_meancbf = File(exists=True, mandatory=True, desc="cbf img")
-    in_avgscore = File(exists=True, mandatory=False, desc="cbf img")
-    in_scrub = File(exists=True, mandatory=False, desc="cbf img")
-    in_basil = File(exists=True, mandatory=False, desc="cbf img")
-    in_pvc = File(exists=True, mandatory=False, desc="cbf img")
-    in_greyM = File(exists=True, mandatory=True, desc="grey  matter")
-    in_whiteM = File(exists=True, mandatory=True, desc="white  matter")
-    in_csf = File(exists=True, mandatory=True, desc="csf")
-    in_confmat = File(exists=True, mandatory=False, desc=" cofnound matrix")
-    in_aslmask = File(exists=True, mandatory=True, desc="asl mask in native space")
-    in_t1mask = File(exists=True, mandatory=True, desc="t1wmask in native space ")
-    in_aslmaskstd = File(exists=True, mandatory=False, desc="asl mask in native space")
-    in_templatemask = File(exists=True, mandatory=False, desc="template mask or image")
-    rmsd_file = File(exists=True, mandatory=True, desc="rmsd file")
-
-
-class _ComputeCBFQCOutputSpec(TraitedSpec):
-    qc_file = File(exists=False, desc="qc file ")
-
-
-class ComputeCBFQC(SimpleInterface):
-    """Calculate a series of CBF quality control metrics for non-GE data.
-
-    compute qc from confound regressors
-    and cbf maps,
-    coregistration and regsitration indexes
-    """
-
-    input_spec = _ComputeCBFQCInputSpec
-    output_spec = _ComputeCBFQCOutputSpec
-
-    def _run_interface(self, runtime):
-        time1 = pd.read_table(self.inputs.in_confmat)
-        time1.fillna(0, inplace=True)
-        fd = np.mean(time1["framewise_displacement"])
-        # rms = time1[['rot_x', 'rot_y', 'rot_z']]
-        # rms1 = rms.pow(2)
-        # rms = np.mean(np.sqrt(rms1.sum(axis=1)/3))
-        print(self.inputs.rmsd_file)
-        rms = pd.read_csv(self.inputs.rmsd_file, header=None).mean().values[0]
-        regDC = dice(self.inputs.in_aslmask, self.inputs.in_t1mask)
-        regJC = jaccard(self.inputs.in_aslmask, self.inputs.in_t1mask)
-        regCC = crosscorr(self.inputs.in_aslmask, self.inputs.in_t1mask)
-        regCov = coverage(self.inputs.in_aslmask, self.inputs.in_t1mask)
-
-        if self.inputs.in_aslmaskstd and self.inputs.in_templatemask:
-            normDC = dice(self.inputs.in_aslmaskstd, self.inputs.in_templatemask)
-            normJC = jaccard(self.inputs.in_aslmaskstd, self.inputs.in_templatemask)
-            normCC = crosscorr(self.inputs.in_aslmaskstd, self.inputs.in_templatemask)
-            normCov = coverage(self.inputs.in_aslmaskstd, self.inputs.in_templatemask)
-
-        meancbf_qei = cbf_qei(
-            gm=self.inputs.in_greyM,
-            wm=self.inputs.in_whiteM,
-            csf=self.inputs.in_csf,
-            img=self.inputs.in_meancbf,
-            thresh=0.7,
-        )
-        meancbf = globalcbf(
-            gm=self.inputs.in_greyM,
-            wm=self.inputs.in_whiteM,
-            csf=self.inputs.in_csf,
-            cbf=self.inputs.in_meancbf,
-            thresh=0.7,
-        )
-
-        if self.inputs.in_avgscore:
-            scorecbf_qei = cbf_qei(
-                gm=self.inputs.in_greyM,
-                wm=self.inputs.in_whiteM,
-                csf=self.inputs.in_csf,
-                img=self.inputs.in_avgscore,
-                thresh=0.7,
-            )
-            scrub_qei = cbf_qei(
-                gm=self.inputs.in_greyM,
-                wm=self.inputs.in_whiteM,
-                csf=self.inputs.in_csf,
-                img=self.inputs.in_scrub,
-                thresh=0.7,
-            )
-            negscore = negativevoxel(
-                cbf=self.inputs.in_avgscore, gm=self.inputs.in_greyM, thresh=0.7
-            )
-            negscrub = negativevoxel(cbf=self.inputs.in_scrub, gm=self.inputs.in_greyM, thresh=0.7)
-        else:
-            print("no score inputs, setting to np.nan")
-            scorecbf_qei = np.nan
-            scrub_qei = np.nan
-            negscore = np.nan
-            negscrub = np.nan
-
-        if self.inputs.in_basil:
-            basilcbf_qei = cbf_qei(
-                gm=self.inputs.in_greyM,
-                wm=self.inputs.in_whiteM,
-                csf=self.inputs.in_csf,
-                img=self.inputs.in_basil,
-                thresh=0.7,
-            )
-            pvcbf_qei = cbf_qei(
-                gm=self.inputs.in_greyM,
-                wm=self.inputs.in_whiteM,
-                csf=self.inputs.in_csf,
-                img=self.inputs.in_pvc,
-                thresh=0.7,
-            )
-            negbasil = negativevoxel(cbf=self.inputs.in_basil, gm=self.inputs.in_greyM, thresh=0.7)
-            negpvc = negativevoxel(cbf=self.inputs.in_pvc, gm=self.inputs.in_greyM, thresh=0.7)
-        else:
-            print("no basil inputs, setting to np.nan")
-            basilcbf_qei = np.nan
-            pvcbf_qei = np.nan
-            negbasil = np.nan
-            negpvc = np.nan
-
-        gwratio = np.divide(meancbf[0], meancbf[1])
-        negcbf = negativevoxel(cbf=self.inputs.in_meancbf, gm=self.inputs.in_greyM, thresh=0.7)
-
-        if self.inputs.in_aslmaskstd and self.inputs.in_templatemask:
-            dict1 = {
-                "FD": [fd],
-                "rmsd": [rms],
-                "coregDC": [regDC],
-                "coregJC": [regJC],
-                "coregCC": [regCC],
-                "coregCOV": [regCov],
-                "normDC": [normDC],
-                "normJC": [normJC],
-                "normCC": [normCC],
-                "normCOV": [normCov],
-                "cbfQEI": [meancbf_qei],
-                "scoreQEI": [scorecbf_qei],
-                "scrubQEI": [scrub_qei],
-                "basilQEI": [basilcbf_qei],
-                "pvcQEI": [pvcbf_qei],
-                "GMmeanCBF": [meancbf[0]],
-                "WMmeanCBF": [meancbf[1]],
-                "Gm_Wm_CBF_ratio": [gwratio],
-                "NEG_CBF_PERC": [negcbf],
-                "NEG_SCORE_PERC": [negscore],
-                "NEG_SCRUB_PERC": [negscrub],
-                "NEG_BASIL_PERC": [negbasil],
-                "NEG_PVC_PERC": [negpvc],
-            }
-        else:
-            dict1 = {
-                "FD": [fd],
-                "rmsd": [rms],
-                "coregDC": [regDC],
-                "coregJC": [regJC],
-                "coregCC": [regCC],
-                "coregCOV": [regCov],
-                "cbfQEI": [meancbf_qei],
-                "scoreQEI": [scorecbf_qei],
-                "scrubQEI": [scrub_qei],
-                "basilQEI": [basilcbf_qei],
-                "pvcQEI": [pvcbf_qei],
-                "GMmeanCBF": [meancbf[0]],
-                "WMmeanCBF": [meancbf[1]],
-                "Gm_Wm_CBF_ratio": [gwratio],
-                "NEG_CBF_PERC": [negcbf],
-                "NEG_SCORE_PERC": [negscore],
-                "NEG_SCRUB_PERC": [negscrub],
-                "NEG_BASIL_PERC": [negbasil],
-                "NEG_PVC_PERC": [negpvc],
-            }
-        _, file1 = os.path.split(self.inputs.in_file)
-        bb = file1.split("_")
-        dict2 = {}
-        for i in range(len(bb) - 1):
-            dict2.update({bb[i].split("-")[0]: bb[i].split("-")[1]})
-        dict2.update(dict1)
-
-        df = pd.DataFrame(dict2)
-
-        self._results["qc_file"] = fname_presuffix(
-            self.inputs.in_meancbf,
-            suffix="qc_cbf.csv",
-            newpath=runtime.cwd,
-            use_ext=False,
-        )
-        df.to_csv(self._results["qc_file"], index=False, header=True)
-
-        return runtime
-
-
-class _ComputeCBFQCforGEInputSpec(BaseInterfaceInputSpec):
-    in_file = File(exists=True, mandatory=True, desc="original asl_file")
-    in_meancbf = File(exists=True, mandatory=True, desc="cbf img")
-    in_avgscore = File(exists=True, mandatory=False, desc="cbf img")
-    in_scrub = File(exists=True, mandatory=False, desc="cbf img")
-    in_basil = File(exists=True, mandatory=False, desc="cbf img")
-    in_pvc = File(exists=True, mandatory=False, desc="cbf img")
-    in_greyM = File(exists=True, mandatory=True, desc="grey  matter")
-    in_whiteM = File(exists=True, mandatory=True, desc="white  matter")
-    in_csf = File(exists=True, mandatory=True, desc="csf")
-    in_aslmask = File(exists=True, mandatory=True, desc="asl mask in native space")
-    in_t1mask = File(exists=True, mandatory=True, desc="t1wmask in native space ")
-    in_aslmaskstd = File(exists=True, mandatory=False, desc="asl mask in native space")
-    in_templatemask = File(exists=True, mandatory=False, desc="template mask or image")
-
-
-class _ComputeCBFQCforGEOutputSpec(TraitedSpec):
-    qc_file = File(exists=False, desc="qc file ")
-
-
-class ComputeCBFQCforGE(SimpleInterface):
-    """Calculate a series of CBF quality control metrics for GE data.
-
-    compute qc from confound regressors
-    and cbf maps,
-    coregistration and regsitration indexes
-    """
-
-    input_spec = _ComputeCBFQCforGEInputSpec
-    output_spec = _ComputeCBFQCforGEOutputSpec
-
-    def _run_interface(self, runtime):
-        regDC = dice(self.inputs.in_aslmask, self.inputs.in_t1mask)
-        regJC = jaccard(self.inputs.in_aslmask, self.inputs.in_t1mask)
-        regCC = crosscorr(self.inputs.in_aslmask, self.inputs.in_t1mask)
-        regCov = coverage(self.inputs.in_aslmask, self.inputs.in_t1mask)
-
-        if self.inputs.in_aslmaskstd and self.inputs.in_templatemask:
-            normDC = dice(self.inputs.in_aslmaskstd, self.inputs.in_templatemask)
-            normJC = jaccard(self.inputs.in_aslmaskstd, self.inputs.in_templatemask)
-            normCC = crosscorr(self.inputs.in_aslmaskstd, self.inputs.in_templatemask)
-            normCov = coverage(self.inputs.in_aslmaskstd, self.inputs.in_templatemask)
-
-        meancbf_qei = cbf_qei(
-            gm=self.inputs.in_greyM,
-            wm=self.inputs.in_whiteM,
-            csf=self.inputs.in_csf,
-            img=self.inputs.in_meancbf,
-            thresh=0.8,
-        )
-        meancbf = globalcbf(
-            gm=self.inputs.in_greyM,
-            wm=self.inputs.in_whiteM,
-            csf=self.inputs.in_csf,
-            cbf=self.inputs.in_meancbf,
-            thresh=0.8,
-        )
-
-        if self.inputs.in_avgscore:
-            scorecbf_qei = cbf_qei(
-                gm=self.inputs.in_greyM,
-                wm=self.inputs.in_whiteM,
-                csf=self.inputs.in_csf,
-                img=self.inputs.in_avgscore,
-                thresh=0.8,
-            )
-            scrub_qei = cbf_qei(
-                gm=self.inputs.in_greyM,
-                wm=self.inputs.in_whiteM,
-                csf=self.inputs.in_csf,
-                img=self.inputs.in_scrub,
-                thresh=0.8,
-            )
-            negscore = negativevoxel(
-                cbf=self.inputs.in_avgscore, gm=self.inputs.in_greyM, thresh=0.8
-            )
-            negscrub = negativevoxel(cbf=self.inputs.in_scrub, gm=self.inputs.in_greyM, thresh=0.8)
-        else:
-            scorecbf_qei = 0
-            scrub_qei = 0
-            negscore = 0
-            negscrub = 0
-
-        if self.inputs.in_basil:
-            basilcbf_qei = cbf_qei(
-                gm=self.inputs.in_greyM,
-                wm=self.inputs.in_whiteM,
-                csf=self.inputs.in_csf,
-                img=self.inputs.in_basil,
-                thresh=0.8,
-            )
-            pvcbf_qei = cbf_qei(
-                gm=self.inputs.in_greyM,
-                wm=self.inputs.in_whiteM,
-                csf=self.inputs.in_csf,
-                img=self.inputs.in_pvc,
-                thresh=0.8,
-            )
-            negbasil = negativevoxel(cbf=self.inputs.in_basil, gm=self.inputs.in_greyM, thresh=0.8)
-            negpvc = negativevoxel(cbf=self.inputs.in_pvc, gm=self.inputs.in_greyM, thresh=0.8)
-        else:
-            basilcbf_qei = 0
-            pvcbf_qei = 0
-            negbasil = 0
-            negpvc = 0
-        gwratio = np.divide(meancbf[0], meancbf[1])
-        negcbf = negativevoxel(cbf=self.inputs.in_meancbf, gm=self.inputs.in_greyM, thresh=0.8)
-
-        if self.inputs.in_aslmaskstd and self.inputs.in_templatemask:
-            dict1 = {
-                "FD": 0,
-                "relRMS": 0,
-                "rmsd": np.nan,
-                "coregDC": [regDC],
-                "coregJC": [regJC],
-                "coregCC": [regCC],
-                "coregCOV": [regCov],
-                "normDC": [normDC],
-                "normJC": [normJC],
-                "normCC": [normCC],
-                "normCOV": [normCov],
-                "cbfQEI": [meancbf_qei],
-                "scoreQEI": [scorecbf_qei],
-                "scrubQEI": [scrub_qei],
-                "basilQEI": [basilcbf_qei],
-                "pvcQEI": [pvcbf_qei],
-                "GMmeanCBF": [meancbf[0]],
-                "WMmeanCBF": [meancbf[1]],
-                "Gm_Wm_CBF_ratio": [gwratio],
-                "NEG_CBF_PERC": [negcbf],
-                "NEG_SCORE_PERC": [negscore],
-                "NEG_SCRUB_PERC": [negscrub],
-                "NEG_BASIL_PERC": [negbasil],
-                "NEG_PVC_PERC": [negpvc],
-            }
-        else:
-            dict1 = {
-                "FD": 0,
-                "relRMS": 0,
-                "rmsd": np.nan,
-                "coregDC": [regDC],
-                "coregJC": [regJC],
-                "coregCC": [regCC],
-                "coregCOV": [regCov],
-                "cbfQEI": [meancbf_qei],
-                "scoreQEI": [scorecbf_qei],
-                "scrubQEI": [scrub_qei],
-                "basilQEI": [basilcbf_qei],
-                "pvcQEI": [pvcbf_qei],
-                "GMmeanCBF": [meancbf[0]],
-                "WMmeanCBF": [meancbf[1]],
-                "Gm_Wm_CBF_ratio": [gwratio],
-                "NEG_CBF_PERC": [negcbf],
-                "NEG_SCORE_PERC": [negscore],
-                "NEG_SCRUB_PERC": [negscrub],
-                "NEG_BASIL_PERC": [negbasil],
-                "NEG_PVC_PERC": [negpvc],
-            }
-        _, file1 = os.path.split(self.inputs.in_file)
-        bb = file1.split("_")
-        dict2 = {}
-        for i in range(len(bb) - 1):
-            dict2.update({bb[i].split("-")[0]: bb[i].split("-")[1]})
-        dict2.update(dict1)
-
-        df = pd.DataFrame(dict2)
-
-        self._results["qc_file"] = fname_presuffix(
-            self.inputs.in_meancbf,
-            suffix="qc_cbf.csv",
-            newpath=runtime.cwd,
-            use_ext=False,
-        )
-        df.to_csv(self._results["qc_file"], index=False, header=True)
-
-        return runtime
-
-
 class _ParcellateCBFInputSpec(BaseInterfaceInputSpec):
     in_cbf = File(exists=True, mandatory=True, desc="cbf img")
     atlasfile = File(exists=True, mandatory=True, desc="data")
diff --git a/aslprep/interfaces/qc.py b/aslprep/interfaces/qc.py
new file mode 100644
index 000000000..6027f6d6b
--- /dev/null
+++ b/aslprep/interfaces/qc.py
@@ -0,0 +1,389 @@
+"""Interfaces for calculating QC measures."""
+import os
+
+import numpy as np
+import pandas as pd
+from nipype.interfaces.base import (
+    BaseInterfaceInputSpec,
+    File,
+    SimpleInterface,
+    TraitedSpec,
+)
+from nipype.utils.filemanip import fname_presuffix
+
+from aslprep.utils.qc import (
+    cbf_qei,
+    coverage,
+    crosscorr,
+    dice,
+    globalcbf,
+    jaccard,
+    negativevoxel,
+)
+
+
+class _ComputeCBFQCInputSpec(BaseInterfaceInputSpec):
+    in_file = File(exists=True, mandatory=True, desc="original asl_file")
+    in_meancbf = File(exists=True, mandatory=True, desc="cbf img")
+    in_avgscore = File(exists=True, mandatory=False, desc="cbf img")
+    in_scrub = File(exists=True, mandatory=False, desc="cbf img")
+    in_basil = File(exists=True, mandatory=False, desc="cbf img")
+    in_pvc = File(exists=True, mandatory=False, desc="cbf img")
+    in_greyM = File(exists=True, mandatory=True, desc="grey  matter")
+    in_whiteM = File(exists=True, mandatory=True, desc="white  matter")
+    in_csf = File(exists=True, mandatory=True, desc="csf")
+    in_confmat = File(exists=True, mandatory=False, desc=" cofnound matrix")
+    in_aslmask = File(exists=True, mandatory=True, desc="asl mask in native space")
+    in_t1mask = File(exists=True, mandatory=True, desc="t1wmask in native space ")
+    in_aslmaskstd = File(exists=True, mandatory=False, desc="asl mask in native space")
+    in_templatemask = File(exists=True, mandatory=False, desc="template mask or image")
+    rmsd_file = File(exists=True, mandatory=True, desc="rmsd file")
+
+
+class _ComputeCBFQCOutputSpec(TraitedSpec):
+    qc_file = File(exists=False, desc="qc file ")
+
+
+class ComputeCBFQC(SimpleInterface):
+    """Calculate a series of CBF quality control metrics for non-GE data.
+
+    compute qc from confound regressors
+    and cbf maps,
+    coregistration and regsitration indexes
+    """
+
+    input_spec = _ComputeCBFQCInputSpec
+    output_spec = _ComputeCBFQCOutputSpec
+
+    def _run_interface(self, runtime):
+        time1 = pd.read_table(self.inputs.in_confmat)
+        time1.fillna(0, inplace=True)
+        fd = np.mean(time1["framewise_displacement"])
+        # rms = time1[['rot_x', 'rot_y', 'rot_z']]
+        # rms1 = rms.pow(2)
+        # rms = np.mean(np.sqrt(rms1.sum(axis=1)/3))
+        print(self.inputs.rmsd_file)
+        rms = pd.read_csv(self.inputs.rmsd_file, header=None).mean().values[0]
+        regDC = dice(self.inputs.in_aslmask, self.inputs.in_t1mask)
+        regJC = jaccard(self.inputs.in_aslmask, self.inputs.in_t1mask)
+        regCC = crosscorr(self.inputs.in_aslmask, self.inputs.in_t1mask)
+        regCov = coverage(self.inputs.in_aslmask, self.inputs.in_t1mask)
+
+        if self.inputs.in_aslmaskstd and self.inputs.in_templatemask:
+            normDC = dice(self.inputs.in_aslmaskstd, self.inputs.in_templatemask)
+            normJC = jaccard(self.inputs.in_aslmaskstd, self.inputs.in_templatemask)
+            normCC = crosscorr(self.inputs.in_aslmaskstd, self.inputs.in_templatemask)
+            normCov = coverage(self.inputs.in_aslmaskstd, self.inputs.in_templatemask)
+
+        meancbf_qei = cbf_qei(
+            gm=self.inputs.in_greyM,
+            wm=self.inputs.in_whiteM,
+            csf=self.inputs.in_csf,
+            img=self.inputs.in_meancbf,
+            thresh=0.7,
+        )
+        meancbf = globalcbf(
+            gm=self.inputs.in_greyM,
+            wm=self.inputs.in_whiteM,
+            csf=self.inputs.in_csf,
+            cbf=self.inputs.in_meancbf,
+            thresh=0.7,
+        )
+
+        if self.inputs.in_avgscore:
+            scorecbf_qei = cbf_qei(
+                gm=self.inputs.in_greyM,
+                wm=self.inputs.in_whiteM,
+                csf=self.inputs.in_csf,
+                img=self.inputs.in_avgscore,
+                thresh=0.7,
+            )
+            scrub_qei = cbf_qei(
+                gm=self.inputs.in_greyM,
+                wm=self.inputs.in_whiteM,
+                csf=self.inputs.in_csf,
+                img=self.inputs.in_scrub,
+                thresh=0.7,
+            )
+            negscore = negativevoxel(
+                cbf=self.inputs.in_avgscore, gm=self.inputs.in_greyM, thresh=0.7
+            )
+            negscrub = negativevoxel(cbf=self.inputs.in_scrub, gm=self.inputs.in_greyM, thresh=0.7)
+        else:
+            print("no score inputs, setting to np.nan")
+            scorecbf_qei = np.nan
+            scrub_qei = np.nan
+            negscore = np.nan
+            negscrub = np.nan
+
+        if self.inputs.in_basil:
+            basilcbf_qei = cbf_qei(
+                gm=self.inputs.in_greyM,
+                wm=self.inputs.in_whiteM,
+                csf=self.inputs.in_csf,
+                img=self.inputs.in_basil,
+                thresh=0.7,
+            )
+            pvcbf_qei = cbf_qei(
+                gm=self.inputs.in_greyM,
+                wm=self.inputs.in_whiteM,
+                csf=self.inputs.in_csf,
+                img=self.inputs.in_pvc,
+                thresh=0.7,
+            )
+            negbasil = negativevoxel(cbf=self.inputs.in_basil, gm=self.inputs.in_greyM, thresh=0.7)
+            negpvc = negativevoxel(cbf=self.inputs.in_pvc, gm=self.inputs.in_greyM, thresh=0.7)
+        else:
+            print("no basil inputs, setting to np.nan")
+            basilcbf_qei = np.nan
+            pvcbf_qei = np.nan
+            negbasil = np.nan
+            negpvc = np.nan
+
+        gwratio = np.divide(meancbf[0], meancbf[1])
+        negcbf = negativevoxel(cbf=self.inputs.in_meancbf, gm=self.inputs.in_greyM, thresh=0.7)
+
+        if self.inputs.in_aslmaskstd and self.inputs.in_templatemask:
+            dict1 = {
+                "FD": [fd],
+                "rmsd": [rms],
+                "coregDC": [regDC],
+                "coregJC": [regJC],
+                "coregCC": [regCC],
+                "coregCOV": [regCov],
+                "normDC": [normDC],
+                "normJC": [normJC],
+                "normCC": [normCC],
+                "normCOV": [normCov],
+                "cbfQEI": [meancbf_qei],
+                "scoreQEI": [scorecbf_qei],
+                "scrubQEI": [scrub_qei],
+                "basilQEI": [basilcbf_qei],
+                "pvcQEI": [pvcbf_qei],
+                "GMmeanCBF": [meancbf[0]],
+                "WMmeanCBF": [meancbf[1]],
+                "Gm_Wm_CBF_ratio": [gwratio],
+                "NEG_CBF_PERC": [negcbf],
+                "NEG_SCORE_PERC": [negscore],
+                "NEG_SCRUB_PERC": [negscrub],
+                "NEG_BASIL_PERC": [negbasil],
+                "NEG_PVC_PERC": [negpvc],
+            }
+        else:
+            dict1 = {
+                "FD": [fd],
+                "rmsd": [rms],
+                "coregDC": [regDC],
+                "coregJC": [regJC],
+                "coregCC": [regCC],
+                "coregCOV": [regCov],
+                "cbfQEI": [meancbf_qei],
+                "scoreQEI": [scorecbf_qei],
+                "scrubQEI": [scrub_qei],
+                "basilQEI": [basilcbf_qei],
+                "pvcQEI": [pvcbf_qei],
+                "GMmeanCBF": [meancbf[0]],
+                "WMmeanCBF": [meancbf[1]],
+                "Gm_Wm_CBF_ratio": [gwratio],
+                "NEG_CBF_PERC": [negcbf],
+                "NEG_SCORE_PERC": [negscore],
+                "NEG_SCRUB_PERC": [negscrub],
+                "NEG_BASIL_PERC": [negbasil],
+                "NEG_PVC_PERC": [negpvc],
+            }
+        _, file1 = os.path.split(self.inputs.in_file)
+        bb = file1.split("_")
+        dict2 = {}
+        for i in range(len(bb) - 1):
+            dict2.update({bb[i].split("-")[0]: bb[i].split("-")[1]})
+        dict2.update(dict1)
+
+        df = pd.DataFrame(dict2)
+
+        self._results["qc_file"] = fname_presuffix(
+            self.inputs.in_meancbf,
+            suffix="qc_cbf.csv",
+            newpath=runtime.cwd,
+            use_ext=False,
+        )
+        df.to_csv(self._results["qc_file"], index=False, header=True)
+
+        return runtime
+
+
+class _ComputeCBFQCforGEInputSpec(BaseInterfaceInputSpec):
+    in_file = File(exists=True, mandatory=True, desc="original asl_file")
+    in_meancbf = File(exists=True, mandatory=True, desc="cbf img")
+    in_avgscore = File(exists=True, mandatory=False, desc="cbf img")
+    in_scrub = File(exists=True, mandatory=False, desc="cbf img")
+    in_basil = File(exists=True, mandatory=False, desc="cbf img")
+    in_pvc = File(exists=True, mandatory=False, desc="cbf img")
+    in_greyM = File(exists=True, mandatory=True, desc="grey  matter")
+    in_whiteM = File(exists=True, mandatory=True, desc="white  matter")
+    in_csf = File(exists=True, mandatory=True, desc="csf")
+    in_aslmask = File(exists=True, mandatory=True, desc="asl mask in native space")
+    in_t1mask = File(exists=True, mandatory=True, desc="t1wmask in native space ")
+    in_aslmaskstd = File(exists=True, mandatory=False, desc="asl mask in native space")
+    in_templatemask = File(exists=True, mandatory=False, desc="template mask or image")
+
+
+class _ComputeCBFQCforGEOutputSpec(TraitedSpec):
+    qc_file = File(exists=False, desc="qc file ")
+
+
+class ComputeCBFQCforGE(SimpleInterface):
+    """Calculate a series of CBF quality control metrics for GE data.
+
+    compute qc from confound regressors
+    and cbf maps,
+    coregistration and regsitration indexes
+    """
+
+    input_spec = _ComputeCBFQCforGEInputSpec
+    output_spec = _ComputeCBFQCforGEOutputSpec
+
+    def _run_interface(self, runtime):
+        regDC = dice(self.inputs.in_aslmask, self.inputs.in_t1mask)
+        regJC = jaccard(self.inputs.in_aslmask, self.inputs.in_t1mask)
+        regCC = crosscorr(self.inputs.in_aslmask, self.inputs.in_t1mask)
+        regCov = coverage(self.inputs.in_aslmask, self.inputs.in_t1mask)
+
+        if self.inputs.in_aslmaskstd and self.inputs.in_templatemask:
+            normDC = dice(self.inputs.in_aslmaskstd, self.inputs.in_templatemask)
+            normJC = jaccard(self.inputs.in_aslmaskstd, self.inputs.in_templatemask)
+            normCC = crosscorr(self.inputs.in_aslmaskstd, self.inputs.in_templatemask)
+            normCov = coverage(self.inputs.in_aslmaskstd, self.inputs.in_templatemask)
+
+        meancbf_qei = cbf_qei(
+            gm=self.inputs.in_greyM,
+            wm=self.inputs.in_whiteM,
+            csf=self.inputs.in_csf,
+            img=self.inputs.in_meancbf,
+            thresh=0.8,
+        )
+        meancbf = globalcbf(
+            gm=self.inputs.in_greyM,
+            wm=self.inputs.in_whiteM,
+            csf=self.inputs.in_csf,
+            cbf=self.inputs.in_meancbf,
+            thresh=0.8,
+        )
+
+        if self.inputs.in_avgscore:
+            scorecbf_qei = cbf_qei(
+                gm=self.inputs.in_greyM,
+                wm=self.inputs.in_whiteM,
+                csf=self.inputs.in_csf,
+                img=self.inputs.in_avgscore,
+                thresh=0.8,
+            )
+            scrub_qei = cbf_qei(
+                gm=self.inputs.in_greyM,
+                wm=self.inputs.in_whiteM,
+                csf=self.inputs.in_csf,
+                img=self.inputs.in_scrub,
+                thresh=0.8,
+            )
+            negscore = negativevoxel(
+                cbf=self.inputs.in_avgscore, gm=self.inputs.in_greyM, thresh=0.8
+            )
+            negscrub = negativevoxel(cbf=self.inputs.in_scrub, gm=self.inputs.in_greyM, thresh=0.8)
+        else:
+            scorecbf_qei = 0
+            scrub_qei = 0
+            negscore = 0
+            negscrub = 0
+
+        if self.inputs.in_basil:
+            basilcbf_qei = cbf_qei(
+                gm=self.inputs.in_greyM,
+                wm=self.inputs.in_whiteM,
+                csf=self.inputs.in_csf,
+                img=self.inputs.in_basil,
+                thresh=0.8,
+            )
+            pvcbf_qei = cbf_qei(
+                gm=self.inputs.in_greyM,
+                wm=self.inputs.in_whiteM,
+                csf=self.inputs.in_csf,
+                img=self.inputs.in_pvc,
+                thresh=0.8,
+            )
+            negbasil = negativevoxel(cbf=self.inputs.in_basil, gm=self.inputs.in_greyM, thresh=0.8)
+            negpvc = negativevoxel(cbf=self.inputs.in_pvc, gm=self.inputs.in_greyM, thresh=0.8)
+        else:
+            basilcbf_qei = 0
+            pvcbf_qei = 0
+            negbasil = 0
+            negpvc = 0
+        gwratio = np.divide(meancbf[0], meancbf[1])
+        negcbf = negativevoxel(cbf=self.inputs.in_meancbf, gm=self.inputs.in_greyM, thresh=0.8)
+
+        if self.inputs.in_aslmaskstd and self.inputs.in_templatemask:
+            dict1 = {
+                "FD": 0,
+                "relRMS": 0,
+                "rmsd": np.nan,
+                "coregDC": [regDC],
+                "coregJC": [regJC],
+                "coregCC": [regCC],
+                "coregCOV": [regCov],
+                "normDC": [normDC],
+                "normJC": [normJC],
+                "normCC": [normCC],
+                "normCOV": [normCov],
+                "cbfQEI": [meancbf_qei],
+                "scoreQEI": [scorecbf_qei],
+                "scrubQEI": [scrub_qei],
+                "basilQEI": [basilcbf_qei],
+                "pvcQEI": [pvcbf_qei],
+                "GMmeanCBF": [meancbf[0]],
+                "WMmeanCBF": [meancbf[1]],
+                "Gm_Wm_CBF_ratio": [gwratio],
+                "NEG_CBF_PERC": [negcbf],
+                "NEG_SCORE_PERC": [negscore],
+                "NEG_SCRUB_PERC": [negscrub],
+                "NEG_BASIL_PERC": [negbasil],
+                "NEG_PVC_PERC": [negpvc],
+            }
+        else:
+            dict1 = {
+                "FD": 0,
+                "relRMS": 0,
+                "rmsd": np.nan,
+                "coregDC": [regDC],
+                "coregJC": [regJC],
+                "coregCC": [regCC],
+                "coregCOV": [regCov],
+                "cbfQEI": [meancbf_qei],
+                "scoreQEI": [scorecbf_qei],
+                "scrubQEI": [scrub_qei],
+                "basilQEI": [basilcbf_qei],
+                "pvcQEI": [pvcbf_qei],
+                "GMmeanCBF": [meancbf[0]],
+                "WMmeanCBF": [meancbf[1]],
+                "Gm_Wm_CBF_ratio": [gwratio],
+                "NEG_CBF_PERC": [negcbf],
+                "NEG_SCORE_PERC": [negscore],
+                "NEG_SCRUB_PERC": [negscrub],
+                "NEG_BASIL_PERC": [negbasil],
+                "NEG_PVC_PERC": [negpvc],
+            }
+        _, file1 = os.path.split(self.inputs.in_file)
+        bb = file1.split("_")
+        dict2 = {}
+        for i in range(len(bb) - 1):
+            dict2.update({bb[i].split("-")[0]: bb[i].split("-")[1]})
+        dict2.update(dict1)
+
+        df = pd.DataFrame(dict2)
+
+        self._results["qc_file"] = fname_presuffix(
+            self.inputs.in_meancbf,
+            suffix="qc_cbf.csv",
+            newpath=runtime.cwd,
+            use_ext=False,
+        )
+        df.to_csv(self._results["qc_file"], index=False, header=True)
+
+        return runtime
diff --git a/aslprep/workflows/asl/base.py b/aslprep/workflows/asl/base.py
index ae8dadde5..ce5d0fa0f 100644
--- a/aslprep/workflows/asl/base.py
+++ b/aslprep/workflows/asl/base.py
@@ -17,7 +17,7 @@
 from aslprep.utils.bids import collect_run_data
 from aslprep.utils.meepi import combine_meepi_source
 from aslprep.utils.misc import _create_mem_gb, _get_series_len, _get_wf_name
-from aslprep.workflows.asl.cbf import init_compute_cbf_wf, init_cbfroiquant_wf
+from aslprep.workflows.asl.cbf import init_cbfroiquant_wf, init_compute_cbf_wf
 from aslprep.workflows.asl.confounds import init_asl_confs_wf, init_carpetplot_wf
 from aslprep.workflows.asl.hmc import init_asl_hmc_wf
 from aslprep.workflows.asl.outputs import init_asl_derivatives_wf
diff --git a/aslprep/workflows/asl/qc.py b/aslprep/workflows/asl/qc.py
index 618bd4203..ad59d53cb 100644
--- a/aslprep/workflows/asl/qc.py
+++ b/aslprep/workflows/asl/qc.py
@@ -6,7 +6,7 @@
 from nipype.pipeline import engine as pe
 from templateflow.api import get as get_template
 
-from aslprep.interfaces.cbf_computation import ComputeCBFQC, ComputeCBFQCforGE
+from aslprep.interfaces.qc import ComputeCBFQC, ComputeCBFQCforGE
 from aslprep.niworkflows.engine.workflows import LiterateWorkflow as Workflow
 from aslprep.niworkflows.interfaces.fixes import (
     FixHeaderApplyTransforms as ApplyTransforms,

From 051d12590a7fae1e23f4521bf9b8bb79ac8520d1 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Tue, 4 Apr 2023 14:00:10 -0400
Subject: [PATCH 12/27] Start refactoring _scrub_csf.

---
 aslprep/interfaces/cbf_computation.py |   4 +-
 aslprep/utils/misc.py                 | 122 +++++++++++++-------------
 2 files changed, 63 insertions(+), 63 deletions(-)

diff --git a/aslprep/interfaces/cbf_computation.py b/aslprep/interfaces/cbf_computation.py
index 84cce13c1..45ea55d99 100644
--- a/aslprep/interfaces/cbf_computation.py
+++ b/aslprep/interfaces/cbf_computation.py
@@ -23,7 +23,7 @@
 from aslprep import config
 from aslprep.utils.misc import (
     _score_cbf,
-    _scrubcbf,
+    _scrub_cbf,
     estimate_labeling_efficiency,
     parcellate_cbf,
     pcasl_or_pasl,
@@ -658,7 +658,7 @@ def _run_interface(self, runtime):
                 thresh=self.inputs.tpm_threshold,
             )
             mean_cbf_score = np.mean(cbf_ts_score, axis=3)
-            mean_cbf_scrub = _scrubcbf(
+            mean_cbf_scrub = _scrub_cbf(
                 cbf_ts=cbf_ts_score,
                 gm=gm_tpm_arr,
                 wm=wm_tpm_arr,
diff --git a/aslprep/utils/misc.py b/aslprep/utils/misc.py
index 5bb88a5de..a0ca97100 100644
--- a/aslprep/utils/misc.py
+++ b/aslprep/utils/misc.py
@@ -608,7 +608,7 @@ def _score_cbf(cbf_ts, wm, gm, csf, mask, thresh=0.7):
 def _robust_fit(
     Y,
     mu,
-    Globalprior,
+    global_prior,
     modrobprior,
     lmd=0,
     localprior=0,
@@ -622,7 +622,7 @@ def _robust_fit(
     priow = np.ones([dimcbf[0], dimcbf[1]])
     sw = 1
     X = priow
-    b = (np.sum(X * Y, axis=0) + mu * Globalprior + lmd * localprior) / (
+    b = (np.sum(X * Y, axis=0) + mu * global_prior + lmd * localprior) / (
         np.sum(X * X, axis=0) + mu + lmd
     )
     b0 = np.repeat(0, len(b))
@@ -658,96 +658,96 @@ def _robust_fit(
         z = np.sqrt(w)
         x = X * z
         yz = Y * z
-        b = (np.sum(x * yz, axis=0) + mu * Globalprior + lmd * localprior) / (
+        b = (np.sum(x * yz, axis=0) + mu * global_prior + lmd * localprior) / (
             np.sum(x * x, axis=0) + mu + lmd
         )
         b = np.nan_to_num(b)
     return b
 
 
-def _scrubcbf(cbf_ts, gm, wm, csf, mask, wfun="huber", thresh=0.7):
+def _scrub_cbf(cbf_ts, gm, wm, csf, mask, wavelet_function="huber", thresh=0.7):
     """Apply SCRUB algorithm to CBF data.
 
     Parameters
     ----------
-    cbf_ts
-       nd array of 3D or 4D computed cbf
-    gm,wm,csf
-       numpy array of grey matter, whitematter, and csf
-    mask
-       numpy array of mask
-    wf
-      wave function
-    """
-    gm = mask * gm
-    gmidx = gm[mask == 1]
-    gmidx[gmidx < thresh] = 0
-    gmidx[gmidx > 0] = 1
+    cbf_ts : 4D numpy.ndarray
+        CBF time series. Generally comes after SCORE censoring.
+    gm, wm, csf : 3D numpy.ndarray
+        Gray matter, white matter, and CSF tissue probability maps.
+    mask : 3D numpy.ndarray
+        Binary brain mask.
+    wavelet_function : str
+        Wavelet function
+    thresh : float
+        Tissue probability threshold.
+        Voxels in gm, wm, and csf with values greater than or equal to the threshold are
+        considered to be of that tissue.
 
-    wm = mask * wm
-    wmidx = wm[mask == 1]
-    wmidx[wmidx < thresh] = 0
-    wmidx[wmidx > 0] = 1
+    Returns
+    -------
+    newcbf : 3D numpy.ndarray
+    """
+    n_volumes = cbf_ts.shape[3]
 
+    # Mask out non-brain voxels from tissue probability maps
+    mask = mask.astype(bool)
+    gm = gm * mask
+    wm = wm * mask
     csf = csf * mask
-    csfidx = csf[mask == 1]
-    csfidx[csfidx < thresh] = 0
-    csfidx[csfidx > 0] = 1
-    # midx = mask[mask==1]
-
-    meancbf = np.mean(cbf_ts, axis=3)
-    y = np.transpose(
-        cbf_ts[
-            mask == 1,
-            :,
-        ]
-    )
-    VV = np.var(y, axis=0)
-    thresh1, thresh3 = _getchisquare(y.shape[0])
-    mu1 = VV / (np.median(VV[gmidx == 1]) * thresh3)
+
+    gm_bin = gm >= thresh
+    gm_idx = gm_bin[mask]
+
+    mean_cbf = np.mean(cbf_ts, axis=3)
+    masked_cbf_ts = cbf_ts[mask, :].T  # TxS array
+    masked_var_map = np.var(masked_cbf_ts, axis=0)  # S array
+    thresh1, thresh3 = _getchisquare(n_volumes)
+    mu1 = masked_var_map / (np.median(masked_var_map[gm_idx == 1]) * thresh3)
     mu = (
         ((mu1 > thresh1) & (mu1 < 10 * thresh1)) * (mu1 - thresh1)
         + (mu1 >= 10 * thresh1) * (1 / (2 * thresh1 * 10) * np.power(mu1, 2))
         + (thresh1 * 10 / 2 - thresh1)
     )
-    M = meancbf * mask
-    M[mask == 1] = mu
+    # Create 3D array from mu
+    mu_map = np.zeros_like(mean_cbf)
+    mu_map[mask] = mu
+
     modrobprior = mu / 10
-    gmidx2 = (
-        1 * ([gm.flatten() > thresh] and [M.flatten() == 0] and [wm.flatten() > csf.flatten()])[0]
-    )
-    wmidx2 = (
-        1 * ([wm.flatten() > thresh] and [M.flatten() == 0] and [gm.flatten() > csf.flatten()])[0]
-    )
-    if np.sum(gmidx2) == 0 or np.sum(wmidx2) == 0:
-        gmidx2 = 1 * (gm.flatten() > thresh)
-        wmidx2 = 1 * (wm.flatten() > thresh)
-    idxx = gmidx2 + wmidx2
+    gm_idx2 = (
+        [gm.flatten() > thresh] and [mu_map.flatten() == 0] and [wm.flatten() > csf.flatten()]
+    )[0].astype(int)
+    wm_idx2 = (
+        [wm.flatten() > thresh] and [mu_map.flatten() == 0] and [gm.flatten() > csf.flatten()]
+    )[0].astype(int)
+    if np.sum(gm_idx2) == 0 or np.sum(wm_idx2) == 0:
+        gm_idx2 = (gm.flatten() > thresh).astype(int)
+        wm_idx2 = (wm.flatten() > thresh).astype(int)
+
+    idxx = gm_idx2 + wm_idx2
     idxx[idxx > 0] = 1
     X = np.zeros([len(idxx), 2])
     X[:, 0] = gm.flatten()[gm.flatten() >= (0)] * idxx
     X[:, 1] = wm.flatten()[wm.flatten() >= (0)] * idxx
-    A = (meancbf.flatten()[idxx >= 0]) * idxx
+    A = (mean_cbf.flatten()[idxx >= 0]) * idxx
     c = np.linalg.lstsq(X, A)[0]
-    Globalpriorfull = c[0] * gm.flatten() + c[1] * wm.flatten()
-    Globalprior = Globalpriorfull[mask.flatten() == 1]
-    localprior = 0
-    lmd = 0
-    tune = _tune(wfun=wfun)
+    global_prior_full = c[0] * gm.flatten() + c[1] * wm.flatten()
+    global_prior = global_prior_full[mask.flatten() == 1]
+
+    tune = _tune(wfun=wavelet_function)
     bb = _robust_fit(
-        Y=y,
+        Y=masked_cbf_ts,
         mu=mu,
-        Globalprior=Globalprior,
+        global_prior=global_prior,
         modrobprior=modrobprior,
-        lmd=lmd,
-        localprior=localprior,
-        wfun=wfun,
+        lmd=0,
+        localprior=0,
+        wfun=wavelet_function,
         tune=tune,
         flagstd=1,
         flagmodrobust=1,
     )
-    newcbf = meancbf * mask
-    newcbf[mask == 1] = bb
+    newcbf = np.zeros_like(mean_cbf)
+    newcbf[mask] = bb
     newcbf = np.nan_to_num(newcbf)
     return newcbf
 

From f44257cd1f9a95c2b85bb3b9bbd7a65625eae367 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Tue, 4 Apr 2023 15:45:35 -0400
Subject: [PATCH 13/27] Fix bug I introduced in SCORE.

---
 aslprep/interfaces/cbf_computation.py |  2 +-
 aslprep/utils/misc.py                 | 43 ++++++++++++---------------
 2 files changed, 20 insertions(+), 25 deletions(-)

diff --git a/aslprep/interfaces/cbf_computation.py b/aslprep/interfaces/cbf_computation.py
index 45ea55d99..f8f7a1177 100644
--- a/aslprep/interfaces/cbf_computation.py
+++ b/aslprep/interfaces/cbf_computation.py
@@ -953,7 +953,7 @@ def refine_ref_mask(t1w_mask, ref_asl_mask, t12ref_transform, tmp_mask, refined_
 
     mat1 = MultiImageMaths()
     mat1.inputs.in_file = tmp_mask
-    mat1.inputs.op_string = " -mul  %s -bin"
+    mat1.inputs.op_string = " -mul %s -bin"
     mat1.inputs.operand_files = ref_asl_mask
     mat1.inputs.out_file = refined_mask
     mat1.run()
diff --git a/aslprep/utils/misc.py b/aslprep/utils/misc.py
index a0ca97100..99635fd39 100644
--- a/aslprep/utils/misc.py
+++ b/aslprep/utils/misc.py
@@ -287,23 +287,16 @@ def _tune(wfun="huber"):
 
     But wait, you might say, the docstring makes no sense! Correct.
     """
-    if wfun == "andrews":
-        tuner = 1.339
-    elif wfun == "bisquare":
-        tuner = 4.685
-    elif wfun == "cauchy":
-        tuner = 2.385
-    elif wfun == "logistic":
-        tuner = 1.205
-    elif wfun == "ols":
-        tuner = 1
-    elif wfun == "talwar":
-        tuner = 2.795
-    elif wfun == "welsch":
-        tuner = 2.985
-    else:
-        tuner = 1.345
-    return tuner
+    weights = {
+        "andrews": 1.339,
+        "bisquare": 4.685,
+        "cauchy": 2.385,
+        "logistic": 1.205,
+        "ols": 1,
+        "talwar": 2.795,
+        "welsch": 2.985,
+    }
+    return weights.get(wfun, 1.345)
 
 
 def _getchisquare(n):
@@ -566,9 +559,11 @@ def _score_cbf(cbf_ts, wm, gm, csf, mask, thresh=0.7):
     outlier_idx = np.abs(mean_gb_cbf - median_gm_cbf) > (2.5 * mad_gm_cbf)
 
     iter_mean_cbf = np.mean(cbf_ts[:, :, :, ~outlier_idx], axis=3)
-    iter_pooled_variance = [
-        (n_voxels[k] * np.var(iter_mean_cbf[tissue_type_masks[k]])) / n_voxels[k] for k in range(3)
-    ]
+    n_voxels_total = np.sum(n_voxels)
+    iter_pooled_variance = (
+        np.sum([n_voxels[k] * np.var(iter_mean_cbf[tissue_type_masks[k]]) for k in range(3)])
+        / n_voxels_total
+    )
     last_pooled_variance = iter_pooled_variance + 1  # add 1 to ensure it is >
     while iter_pooled_variance < last_pooled_variance:
         # If variance increases from the previous iteration, we interpret this as removing
@@ -592,10 +587,10 @@ def _score_cbf(cbf_ts, wm, gm, csf, mask, thresh=0.7):
 
         # Recalculate mean CBF map and pooled variance.
         iter_mean_cbf = np.mean(cbf_ts[:, :, :, ~outlier_idx], axis=3)
-        iter_pooled_variance = [
-            (n_voxels[k] * np.var(iter_mean_cbf[tissue_type_masks[k]])) / n_voxels[k]
-            for k in range(3)
-        ]
+        iter_pooled_variance = (
+            np.sum([n_voxels[k] * np.var(iter_mean_cbf[tissue_type_masks[k]]) for k in range(3)])
+            / n_voxels_total
+        )
 
     # Create CBF time series from just the non-outlier volumes and brain-mask it.
     cbf_ts_good_vols = cbf_ts[:, :, :, ~outlier_idx]

From 7588d344d8261f48054f33d3d34b7a48a12b291e Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Wed, 5 Apr 2023 16:18:34 -0400
Subject: [PATCH 14/27] Update .gitignore

---
 .gitignore | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.gitignore b/.gitignore
index fd7acdd46..a97882433 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,4 @@
+aslprep/tests/data/test_data
 .circleci/local_aslprep_path.txt
 
 .DS_Store

From 984bcbcfb89c5f215d38988fe6af849daa8a333c Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Wed, 5 Apr 2023 16:56:28 -0400
Subject: [PATCH 15/27] Simplify _get_wfun_weight.

---
 aslprep/utils/misc.py | 73 ++++++++++++++++++-------------------------
 1 file changed, 31 insertions(+), 42 deletions(-)

diff --git a/aslprep/utils/misc.py b/aslprep/utils/misc.py
index 99635fd39..a21b45990 100644
--- a/aslprep/utils/misc.py
+++ b/aslprep/utils/misc.py
@@ -253,36 +253,21 @@ def get_tis(metadata: "dict[str, Any]") -> list:
         return np.array(metadata["PostLabelingDelay"]).tolist()
 
 
-def _weightfun(x, wfun="huber"):
-    """Get weight fun and tuner."""
-    if wfun == "andrews":
-        tuner = 1.339
-        weight = (np.abs(x) < np.pi) * np.sin(x)
-    elif wfun == "bisquare":
-        tuner = 4.685
-        weight = (np.abs(x) < 1) * np.power((1 - np.power(x, 2)), 2)
-    elif wfun == "cauchy":
-        tuner = 2.385
-        weight = 1 / (1 + np.power(x, 2))
-    elif wfun == "logistic":
-        tuner = 1.205
-        weight = np.tanh(x) / x
-    elif wfun == "ols":
-        tuner = 1
-        weight = np.repeat(1, len(x))
-    elif wfun == "talwar":
-        tuner = 2.795
-        weight = 1 * (np.abs(x) < 1)
-    elif wfun == "welsch":
-        tuner = 2.985
-        weight = np.exp(-(np.power(x, 2)))
-    else:
-        tuner = 1.345
-        weight = 1 / np.abs(x)
-    return weight, tuner
+def _get_wfun_weight(x, wfun="huber"):
+    """Get wavelet function weight."""
+    weights = {
+        "andrews": (np.abs(x) < np.pi) * np.sin(x),
+        "bisquare": (np.abs(x) < 1) * np.power((1 - np.power(x, 2)), 2),
+        "cauchy": 1 / (1 + np.power(x, 2)),
+        "logistic": np.tanh(x) / x,
+        "ols": np.repeat(1, len(x)),
+        "talwar": 1 * (np.abs(x) < 1),
+        "welsch": np.exp(-(np.power(x, 2))),
+    }
+    return weights.get(wfun, 1 / np.abs(x))
 
 
-def _tune(wfun="huber"):
+def _get_wfun_tuner(wfun="huber"):
     """Get weight fun and tuner.
 
     But wait, you might say, the docstring makes no sense! Correct.
@@ -628,16 +613,17 @@ def _robust_fit(
     tiny_s = (1e-6) * (np.std(h, axis=0))
     tiny_s[tiny_s == 0] = 1
     D = np.sqrt(np.finfo(float).eps)
-    iter_num = 0
-    interlim = 10
-    while iter_num < interlim:
-        print("iteration  ", iter_num, "\n")
-        iter_num = iter_num + 1
+
+    iter_num, max_iters = 0, 10
+    while iter_num < max_iters:
+        print(f"iteration {iter_num}")
+        iter_num += 1
         check1 = np.subtract(np.abs(b - b0), (D * np.maximum(np.abs(b), np.abs(b0))))
         check1[check1 > 0] = 0
         if any(check1):
-            print(" \n converged after ", iter_num, "iterations\n")
+            print(f"converged after {iter_num} iterations")
             break
+
         r = Y - X * (np.tile(b, (dimcbf[0], 1)))
         radj = r * adjfactor / sw
         if flagstd == 1:
@@ -645,10 +631,11 @@ def _robust_fit(
         else:
             rs = np.sort(np.abs(radj), axis=0)
             s = np.median(rs, axis=0) / 0.6745
+
         rx1 = radj * (1 - flagmodrobust * np.exp(-np.tile(modrobprior, (dimcbf[0], 1))))
         rx2 = np.tile(np.maximum(s, tiny_s) * tune, (dimcbf[0], 1))
         r1 = rx1 / rx2
-        w, _ = _weightfun(r1, wfun)
+        w = _get_wfun_weight(r1, wfun)
         b0 = b
         z = np.sqrt(w)
         x = X * z
@@ -657,6 +644,7 @@ def _robust_fit(
             np.sum(x * x, axis=0) + mu + lmd
         )
         b = np.nan_to_num(b)
+
     return b
 
 
@@ -680,7 +668,8 @@ def _scrub_cbf(cbf_ts, gm, wm, csf, mask, wavelet_function="huber", thresh=0.7):
 
     Returns
     -------
-    newcbf : 3D numpy.ndarray
+    mean_cbf_scrub : 3D numpy.ndarray
+        Mean CBF after denoising with SCRUB.
     """
     n_volumes = cbf_ts.shape[3]
 
@@ -728,8 +717,8 @@ def _scrub_cbf(cbf_ts, gm, wm, csf, mask, wavelet_function="huber", thresh=0.7):
     global_prior_full = c[0] * gm.flatten() + c[1] * wm.flatten()
     global_prior = global_prior_full[mask.flatten() == 1]
 
-    tune = _tune(wfun=wavelet_function)
-    bb = _robust_fit(
+    tune = _get_wfun_tuner(wfun=wavelet_function)
+    masked_mean_cbf_scrub = _robust_fit(
         Y=masked_cbf_ts,
         mu=mu,
         global_prior=global_prior,
@@ -741,10 +730,10 @@ def _scrub_cbf(cbf_ts, gm, wm, csf, mask, wavelet_function="huber", thresh=0.7):
         flagstd=1,
         flagmodrobust=1,
     )
-    newcbf = np.zeros_like(mean_cbf)
-    newcbf[mask] = bb
-    newcbf = np.nan_to_num(newcbf)
-    return newcbf
+    mean_cbf_scrub = np.zeros_like(mean_cbf)
+    mean_cbf_scrub[mask] = masked_mean_cbf_scrub
+    mean_cbf_scrub = np.nan_to_num(mean_cbf_scrub)
+    return mean_cbf_scrub
 
 
 def get_atlas(atlasname):

From 6fb627294b3a78a56cb2f4e47c57b4123a2639d5 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Thu, 6 Apr 2023 15:13:50 -0400
Subject: [PATCH 16/27] Continue, not break.

---
 aslprep/utils/misc.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/aslprep/utils/misc.py b/aslprep/utils/misc.py
index a21b45990..eff4abb0a 100644
--- a/aslprep/utils/misc.py
+++ b/aslprep/utils/misc.py
@@ -557,7 +557,7 @@ def _score_cbf(cbf_ts, wm, gm, csf, mask, thresh=0.7):
         volume_wise_correlations = np.full(n_volumes, np.nan)
         for i_vol in range(n_volumes):
             if outlier_idx[i_vol]:
-                break
+                continue
 
             cbf_vol = cbf_ts[:, :, :, i_vol]
             # Correlate mean CBF's masked values with volume CBF's masked values.

From 66ee1fc9c7649f67911bfadfb76a0560c1248157 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Thu, 25 May 2023 15:32:36 -0400
Subject: [PATCH 17/27] Reimplement refactor of SCRUB/SCORE code.

---
 aslprep/utils/cbf.py | 326 ++++++++++++++++++++++---------------------
 1 file changed, 166 insertions(+), 160 deletions(-)

diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py
index 993ce47b9..5f8fdb6f1 100644
--- a/aslprep/utils/cbf.py
+++ b/aslprep/utils/cbf.py
@@ -8,57 +8,35 @@
 from aslprep import config
 
 
-def _weightfun(x, wfun="huber"):
-    """Get weight fun and tuner."""
-    if wfun == "andrews":
-        tuner = 1.339
-        weight = (np.abs(x) < np.pi) * np.sin(x)
-    elif wfun == "bisquare":
-        tuner = 4.685
-        weight = (np.abs(x) < 1) * np.power((1 - np.power(x, 2)), 2)
-    elif wfun == "cauchy":
-        tuner = 2.385
-        weight = 1 / (1 + np.power(x, 2))
-    elif wfun == "logistic":
-        tuner = 1.205
-        weight = np.tanh(x) / x
-    elif wfun == "ols":
-        tuner = 1
-        weight = np.repeat(1, len(x))
-    elif wfun == "talwar":
-        tuner = 2.795
-        weight = 1 * (np.abs(x) < 1)
-    elif wfun == "welsch":
-        tuner = 2.985
-        weight = np.exp(-(np.power(x, 2)))
-    else:
-        tuner = 1.345
-        weight = 1 / np.abs(x)
-    return weight, tuner
+def _get_wfun_weight(x, wfun="huber"):
+    """Get wavelet function weight."""
+    weights = {
+        "andrews": (np.abs(x) < np.pi) * np.sin(x),
+        "bisquare": (np.abs(x) < 1) * np.power((1 - np.power(x, 2)), 2),
+        "cauchy": 1 / (1 + np.power(x, 2)),
+        "logistic": np.tanh(x) / x,
+        "ols": np.repeat(1, len(x)),
+        "talwar": 1 * (np.abs(x) < 1),
+        "welsch": np.exp(-(np.power(x, 2))),
+    }
+    return weights.get(wfun, 1 / np.abs(x))
 
 
-def _tune(wfun="huber"):
+def _get_wfun_tuner(wfun="huber"):
     """Get weight fun and tuner.
 
     But wait, you might say, the docstring makes no sense! Correct.
     """
-    if wfun == "andrews":
-        tuner = 1.339
-    elif wfun == "bisquare":
-        tuner = 4.685
-    elif wfun == "cauchy":
-        tuner = 2.385
-    elif wfun == "logistic":
-        tuner = 1.205
-    elif wfun == "ols":
-        tuner = 1
-    elif wfun == "talwar":
-        tuner = 2.795
-    elif wfun == "welsch":
-        tuner = 2.985
-    else:
-        tuner = 1.345
-    return tuner
+    weights = {
+        "andrews": 1.339,
+        "bisquare": 4.685,
+        "cauchy": 2.385,
+        "logistic": 1.205,
+        "ols": 1,
+        "talwar": 2.795,
+        "welsch": 2.985,
+    }
+    return weights.get(wfun, 1.345)
 
 
 def _getchisquare(n):
@@ -273,75 +251,99 @@ def _getchisquare(n):
     return a[n - 1], b[n - 1]
 
 
-def _getcbfscore(cbfts, wm, gm, csf, mask, thresh=0.7):
+def _score_cbf(cbf_ts, wm, gm, csf, mask, thresh=0.7):
     """Apply SCORE algorithm to remove noisy CBF volumes.
 
     Parameters
     ----------
-    cbf_ts
-       nd array of 3D or 4D computed cbf
-    gm,wm,csf
-       numpy array of grey matter, whitematter, and csf
-    mask
-       numpy array of mask
+    cbf_ts : 4D numpy.ndarray
+        CBF time series.
+    gm, wm, csf : 3D numpy.ndarray
+        Gray matter, white matter, and CSF tissue probability maps.
+    mask : 3D numpy.ndarray
+        Binary brain mask.
+    thresh : float
+        Tissue probability threshold.
+        Voxels in gm, wm, and csf with values greater than or equal to the threshold are
+        considered to be of that tissue.
+
+    Returns
+    -------
+    cbf_ts_good_vols : 4D numpy.ndarray
+        CBF time series containing only non-outlier volumes. This time series is also brain-masked.
+    outlier_idx : 1D numpy.ndarray
+        An index of outlier volumes in the CBF time series. Values of 1 indicate outliers.
     """
+    n_volumes = cbf_ts.shape[3]
+
+    # Create binary tissue-type masks
     gm_bin = gm >= thresh
     wm_bin = wm >= thresh
     csf_bin = csf >= thresh
+    tissue_type_masks = [gm_bin, wm_bin, csf_bin]
 
-    # get the total number of voxle within csf,gm and wm
+    # Get the number of voxels within each tissue type, then subtract 1 to get Nk-1.
     n_gm_voxels = np.sum(gm_bin) - 1
     n_wm_voxels = np.sum(wm_bin) - 1
     n_csf_voxels = np.sum(csf_bin) - 1
-    mask1 = gm_bin + wm_bin + csf_bin
-
-    # mean  of times series cbf within greymatter
-    gm_cbf_ts = np.squeeze(np.mean(cbfts[gm_bin, :], axis=0))
-
-    # robust mean and median
-    median_gm_cbf = np.median(gm_cbf_ts)
-    mad_gm_cbf = median_abs_deviation(gm_cbf_ts) / 0.675
-    indx = 1 * (np.abs(gm_cbf_ts - median_gm_cbf) > (2.5 * mad_gm_cbf))
-    R = np.mean(cbfts[:, :, :, indx == 0], axis=3)
-    V = (
-        n_gm_voxels * np.var(R[gm == 1])
-        + n_wm_voxels * np.var(R[wm == 1])
-        + n_csf_voxels * np.var(R[csf == 1])
+    n_voxels = [n_gm_voxels, n_wm_voxels, n_csf_voxels]
+
+    # Create a mask only including GM, WM, and CSF
+    tissue_type_mask = (gm_bin + wm_bin + csf_bin).astype(bool)
+
+    # Perform initial outlier exclusion.
+    # Remove volumes with mean GM-CBF outside 2.5 SDs of their means.
+    mean_gb_cbf = np.squeeze(np.mean(cbf_ts[gm_bin, :], axis=0))
+    median_gm_cbf = np.median(mean_gb_cbf)
+    mad_gm_cbf = median_abs_deviation(mean_gb_cbf) / 0.675
+    outlier_idx = np.abs(mean_gb_cbf - median_gm_cbf) > (2.5 * mad_gm_cbf)
+
+    iter_mean_cbf = np.mean(cbf_ts[:, :, :, ~outlier_idx], axis=3)
+    n_voxels_total = np.sum(n_voxels)
+    iter_pooled_variance = (
+        np.sum([n_voxels[k] * np.var(iter_mean_cbf[tissue_type_masks[k]]) for k in range(3)])
+        / n_voxels_total
     )
-    V1 = V + 1
-    while V < V1:
-        V1 = V
-        CC = np.zeros(cbfts.shape[3]) * (-2)
-        for s in range(cbfts.shape[3]):
-            if indx[s] != 0:
-                break
-            else:
-                tmp1 = cbfts[:, :, :, s]
-                CC[s] = np.corrcoef(R[mask1 > 0], tmp1[mask1 > 0])[0][1]
-
-        inx = np.argmax(CC)
-        indx[inx] = 2
-        R = np.mean(cbfts[:, :, :, indx == 0], axis=3)
-        V = (
-            (n_gm_voxels * np.var(R[gm == 1]))
-            + (n_wm_voxels * np.var(R[wm == 1]))
-            + (n_csf_voxels * np.var(R[csf == 1]))
+    last_pooled_variance = iter_pooled_variance + 1  # add 1 to ensure it is >
+    while iter_pooled_variance < last_pooled_variance:
+        # If variance increases from the previous iteration, we interpret this as removing
+        # random noise instead of artifact.
+        last_pooled_variance = iter_pooled_variance
+        volume_wise_correlations = np.full(n_volumes, np.nan)
+        for i_vol in range(n_volumes):
+            if outlier_idx[i_vol]:
+                continue
+
+            cbf_vol = cbf_ts[:, :, :, i_vol]
+            # Correlate mean CBF's masked values with volume CBF's masked values.
+            volume_wise_correlations[i_vol] = np.corrcoef(
+                iter_mean_cbf[tissue_type_mask],
+                cbf_vol[tissue_type_mask],
+            )[0, 1]
+
+        # The volume with the maximum correlation value is this iteration's outlier.
+        max_correlation_idx = np.nanargmax(volume_wise_correlations)
+        outlier_idx[max_correlation_idx] = True
+
+        # Recalculate mean CBF map and pooled variance.
+        iter_mean_cbf = np.mean(cbf_ts[:, :, :, ~outlier_idx], axis=3)
+        iter_pooled_variance = (
+            np.sum([n_voxels[k] * np.var(iter_mean_cbf[tissue_type_masks[k]]) for k in range(3)])
+            / n_voxels_total
         )
 
-    config.loggers.utils.warning(f"SCORE retains {np.sum(indx == 0)}/{indx.size} volumes")
-    cbfts_recon = cbfts[:, :, :, indx == 0]
-    cbfts_recon1 = np.zeros_like(cbfts_recon)
-    for i in range(cbfts_recon.shape[3]):
-        cbfts_recon1[:, :, :, i] = cbfts_recon[:, :, :, i] * mask
+    # Create CBF time series from just the non-outlier volumes and brain-mask it.
+    cbf_ts_good_vols = cbf_ts[:, :, :, ~outlier_idx]
+    cbf_ts_good_vols *= mask[..., None]
+    cbf_ts_good_vols = np.nan_to_num(cbf_ts_good_vols)
 
-    cbfts_recon1 = np.nan_to_num(cbfts_recon1)
-    return cbfts_recon1, indx
+    return cbf_ts_good_vols, outlier_idx.astype(int)
 
 
 def _robust_fit(
     Y,
     mu,
-    Globalprior,
+    global_prior,
     modrobprior,
     lmd=0,
     localprior=0,
@@ -355,7 +357,7 @@ def _robust_fit(
     priow = np.ones([dimcbf[0], dimcbf[1]])
     sw = 1
     X = priow
-    b = (np.sum(X * Y, axis=0) + mu * Globalprior + lmd * localprior) / (
+    b = (np.sum(X * Y, axis=0) + mu * global_prior + lmd * localprior) / (
         np.sum(X * X, axis=0) + mu + lmd
     )
     b0 = np.repeat(0, len(b))
@@ -366,16 +368,17 @@ def _robust_fit(
     tiny_s = (1e-6) * (np.std(h, axis=0))
     tiny_s[tiny_s == 0] = 1
     D = np.sqrt(np.finfo(float).eps)
-    iter_num = 0
-    interlim = 10
-    while iter_num < interlim:
-        print("iteration  ", iter_num, "\n")
-        iter_num = iter_num + 1
+
+    iter_num, max_iters = 0, 10
+    while iter_num < max_iters:
+        print(f"iteration {iter_num}")
+        iter_num += 1
         check1 = np.subtract(np.abs(b - b0), (D * np.maximum(np.abs(b), np.abs(b0))))
         check1[check1 > 0] = 0
         if any(check1):
-            print(" \n converged after ", iter_num, "iterations\n")
+            print(f"converged after {iter_num} iterations")
             break
+
         r = Y - X * (np.tile(b, (dimcbf[0], 1)))
         radj = r * adjfactor / sw
         if flagstd == 1:
@@ -383,106 +386,109 @@ def _robust_fit(
         else:
             rs = np.sort(np.abs(radj), axis=0)
             s = np.median(rs, axis=0) / 0.6745
+
         rx1 = radj * (1 - flagmodrobust * np.exp(-np.tile(modrobprior, (dimcbf[0], 1))))
         rx2 = np.tile(np.maximum(s, tiny_s) * tune, (dimcbf[0], 1))
         r1 = rx1 / rx2
-        w, _ = _weightfun(r1, wfun)
+        w = _get_wfun_weight(r1, wfun)
         b0 = b
         z = np.sqrt(w)
         x = X * z
         yz = Y * z
-        b = (np.sum(x * yz, axis=0) + mu * Globalprior + lmd * localprior) / (
+        b = (np.sum(x * yz, axis=0) + mu * global_prior + lmd * localprior) / (
             np.sum(x * x, axis=0) + mu + lmd
         )
         b = np.nan_to_num(b)
+
     return b
 
 
-def _scrubcbf(cbf_ts, gm, wm, csf, mask, wfun="huber", thresh=0.7):
+def _scrub_cbf(cbf_ts, gm, wm, csf, mask, wavelet_function="huber", thresh=0.7):
     """Apply SCRUB algorithm to CBF data.
 
     Parameters
     ----------
-    cbf_ts
-       nd array of 3D or 4D computed cbf
-    gm,wm,csf
-       numpy array of grey matter, whitematter, and csf
-    mask
-       numpy array of mask
-    wf
-      wave function
-    """
-    gm = mask * gm
-    gmidx = gm[mask == 1]
-    gmidx[gmidx < thresh] = 0
-    gmidx[gmidx > 0] = 1
+    cbf_ts : 4D numpy.ndarray
+        CBF time series. Generally comes after SCORE censoring.
+    gm, wm, csf : 3D numpy.ndarray
+        Gray matter, white matter, and CSF tissue probability maps.
+    mask : 3D numpy.ndarray
+        Binary brain mask.
+    wavelet_function : str
+        Wavelet function
+    thresh : float
+        Tissue probability threshold.
+        Voxels in gm, wm, and csf with values greater than or equal to the threshold are
+        considered to be of that tissue.
 
-    wm = mask * wm
-    wmidx = wm[mask == 1]
-    wmidx[wmidx < thresh] = 0
-    wmidx[wmidx > 0] = 1
+    Returns
+    -------
+    mean_cbf_scrub : 3D numpy.ndarray
+        Mean CBF after denoising with SCRUB.
+    """
+    n_volumes = cbf_ts.shape[3]
 
+    # Mask out non-brain voxels from tissue probability maps
+    mask = mask.astype(bool)
+    gm = gm * mask
+    wm = wm * mask
     csf = csf * mask
-    csfidx = csf[mask == 1]
-    csfidx[csfidx < thresh] = 0
-    csfidx[csfidx > 0] = 1
-    # midx = mask[mask==1]
+
+    gm_bin = gm >= thresh
+    gm_idx = gm_bin[mask]
 
     mean_cbf = np.mean(cbf_ts, axis=3)
-    y = np.transpose(
-        cbf_ts[
-            mask == 1,
-            :,
-        ]
-    )
-    VV = np.var(y, axis=0)
-    thresh1, thresh3 = _getchisquare(y.shape[0])
-    mu1 = VV / (np.median(VV[gmidx == 1]) * thresh3)
+    masked_cbf_ts = cbf_ts[mask, :].T  # TxS array
+    masked_var_map = np.var(masked_cbf_ts, axis=0)  # S array
+    thresh1, thresh3 = _getchisquare(n_volumes)
+    mu1 = masked_var_map / (np.median(masked_var_map[gm_idx == 1]) * thresh3)
     mu = (
         ((mu1 > thresh1) & (mu1 < 10 * thresh1)) * (mu1 - thresh1)
         + (mu1 >= 10 * thresh1) * (1 / (2 * thresh1 * 10) * np.power(mu1, 2))
         + (thresh1 * 10 / 2 - thresh1)
     )
-    M = mean_cbf * mask
-    M[mask == 1] = mu
+    # Create 3D array from mu
+    mu_map = np.zeros_like(mean_cbf)
+    mu_map[mask] = mu
+
     modrobprior = mu / 10
-    gmidx2 = (
-        1 * ([gm.flatten() > thresh] and [M.flatten() == 0] and [wm.flatten() > csf.flatten()])[0]
-    )
-    wmidx2 = (
-        1 * ([wm.flatten() > thresh] and [M.flatten() == 0] and [gm.flatten() > csf.flatten()])[0]
-    )
-    if np.sum(gmidx2) == 0 or np.sum(wmidx2) == 0:
-        gmidx2 = 1 * (gm.flatten() > thresh)
-        wmidx2 = 1 * (wm.flatten() > thresh)
-    idxx = gmidx2 + wmidx2
+    gm_idx2 = (
+        [gm.flatten() > thresh] and [mu_map.flatten() == 0] and [wm.flatten() > csf.flatten()]
+    )[0].astype(int)
+    wm_idx2 = (
+        [wm.flatten() > thresh] and [mu_map.flatten() == 0] and [gm.flatten() > csf.flatten()]
+    )[0].astype(int)
+    if np.sum(gm_idx2) == 0 or np.sum(wm_idx2) == 0:
+        gm_idx2 = (gm.flatten() > thresh).astype(int)
+        wm_idx2 = (wm.flatten() > thresh).astype(int)
+
+    idxx = gm_idx2 + wm_idx2
     idxx[idxx > 0] = 1
     X = np.zeros([len(idxx), 2])
     X[:, 0] = gm.flatten()[gm.flatten() >= (0)] * idxx
     X[:, 1] = wm.flatten()[wm.flatten() >= (0)] * idxx
     A = (mean_cbf.flatten()[idxx >= 0]) * idxx
     c = np.linalg.lstsq(X, A)[0]
-    Globalpriorfull = c[0] * gm.flatten() + c[1] * wm.flatten()
-    Globalprior = Globalpriorfull[mask.flatten() == 1]
-    localprior = 0
-    lmd = 0
-    tune = _tune(wfun=wfun)
-    bb = _robust_fit(
-        Y=y,
+    global_prior_full = c[0] * gm.flatten() + c[1] * wm.flatten()
+    global_prior = global_prior_full[mask.flatten() == 1]
+
+    tune = _get_wfun_tuner(wfun=wavelet_function)
+    masked_mean_cbf_scrub = _robust_fit(
+        Y=masked_cbf_ts,
         mu=mu,
-        Globalprior=Globalprior,
+        global_prior=global_prior,
         modrobprior=modrobprior,
-        lmd=lmd,
-        localprior=localprior,
-        wfun=wfun,
+        lmd=0,
+        localprior=0,
+        wfun=wavelet_function,
         tune=tune,
         flagstd=1,
         flagmodrobust=1,
     )
-    newcbf = mean_cbf * mask
-    newcbf[mask == 1] = bb
-    newcbf = np.nan_to_num(newcbf)
-    return newcbf
+    mean_cbf_scrub = np.zeros_like(mean_cbf)
+    mean_cbf_scrub[mask] = masked_mean_cbf_scrub
+    mean_cbf_scrub = np.nan_to_num(mean_cbf_scrub)
+    return mean_cbf_scrub
 
 
 def estimate_att_pcasl(deltam_arr, plds, lds, t1blood, t1tissue):

From bb4c600f33680ce2ae1f5957df7b540cb3f35d9b Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Thu, 25 May 2023 15:34:48 -0400
Subject: [PATCH 18/27] Rename the functions.

---
 aslprep/interfaces/cbf.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/aslprep/interfaces/cbf.py b/aslprep/interfaces/cbf.py
index b8d5e6bfd..8753fae39 100644
--- a/aslprep/interfaces/cbf.py
+++ b/aslprep/interfaces/cbf.py
@@ -27,8 +27,8 @@
     pcasl_or_pasl,
 )
 from aslprep.utils.cbf import (
-    _getcbfscore,
-    _scrubcbf,
+    _score_cbf,
+    _scrub_cbf,
     estimate_cbf_pcasl_multipld,
     estimate_t1,
 )
@@ -724,7 +724,7 @@ def _run_interface(self, runtime):
         whitem = nb.load(self.inputs.wm_tpm).get_fdata()
         csf = nb.load(self.inputs.csf_tpm).get_fdata()
         if cbf_ts.ndim > 3:
-            cbf_scorets, index_score = _getcbfscore(
+            cbf_scorets, index_score = _score_cbf(
                 cbfts=cbf_ts,
                 wm=whitem,
                 gm=greym,
@@ -732,7 +732,7 @@ def _run_interface(self, runtime):
                 mask=mask,
                 thresh=self.inputs.tpm_threshold,
             )
-            cbfscrub = _scrubcbf(
+            cbfscrub = _scrub_cbf(
                 cbf_ts=cbf_scorets,
                 gm=greym,
                 wm=whitem,

From ed3a569feba78fbcf7cb95fda6cf5e11ecdfe981 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Thu, 25 May 2023 15:38:20 -0400
Subject: [PATCH 19/27] Update cbf.py

---
 aslprep/interfaces/cbf.py | 168 ++++++++++++++++++++------------------
 1 file changed, 90 insertions(+), 78 deletions(-)

diff --git a/aslprep/interfaces/cbf.py b/aslprep/interfaces/cbf.py
index 8753fae39..176750bd5 100644
--- a/aslprep/interfaces/cbf.py
+++ b/aslprep/interfaces/cbf.py
@@ -94,11 +94,11 @@ class _ExtractCBFInputSpec(BaseInterfaceInputSpec):
     in_mask = File(exists=True, mandatory=True, desc="mask")
     dummy_vols = traits.Int(
         default_value=0,
-        use_default=True,
+        usedefault=True,
         mandatory=False,
         desc="remove first n volumes",
     )
-    fwhm = traits.Float(default_value=5, use_default=True, mandatory=False, desc="fwhm")
+    fwhm = traits.Float(default_value=5, usedefault=True, mandatory=False, desc="fwhm")
 
 
 class _ExtractCBFOutputSpec(TraitedSpec):
@@ -672,31 +672,37 @@ def _run_interface(self, runtime):
 
 
 class _ScoreAndScrubCBFInputSpec(BaseInterfaceInputSpec):
-    cbf_ts = File(exists=True, mandatory=True, desc="Computed CBF from ComputeCBF.")
-    mask = File(exists=True, mandatory=True, desc="mask")
+    cbf = File(exists=True, mandatory=True, desc="computed CBF from ComputeCBF")
     gm_tpm = File(exists=True, mandatory=True, desc="Gray matter tissue probability map.")
     wm_tpm = File(exists=True, mandatory=True, desc="White matter tissue probability map.")
-    csf_tpm = File(exists=True, mandatory=True, desc="CSF tissue probability map.")
+    csf_tpm = File(exists=True, mandatory=True, desc="Cerebrospinal fluid tissue probability map.")
+    brain_mask = File(exists=True, mandatory=True, desc="Brain mask.")
     tpm_threshold = traits.Float(
         default_value=0.7,
-        usedefault=True,
         mandatory=False,
-        desc="Tissue probability threshold for binarizing GM, WM, and CSF masks.",
+        usedefault=True,
+        desc="Threshold for tissue probability maps.",
     )
     wavelet_function = traits.Str(
-        default_value="huber",
-        usedefault=True,
+        exists=True,
         mandatory=False,
+        default_value="huber",
         option=["bisquare", "andrews", "cauchy", "fair", "logistics", "ols", "talwar", "welsch"],
-        desc="Wavelet function",
+        desc="Wavelet function used in SCRUB.",
     )
 
 
 class _ScoreAndScrubCBFOutputSpec(TraitedSpec):
-    cbf_ts_score = File(exists=False, mandatory=False, desc="score timeseries data")
-    mean_cbf_score = File(exists=False, mandatory=False, desc="average score")
-    mean_cbf_scrub = File(exists=False, mandatory=False, desc="average scrub")
-    score_outlier_index = File(exists=False, mandatory=False, desc="index of volume remove ")
+    cbf_ts_score = File(
+        exists=True,
+        desc="CBF time series after removing outlier volumes flagged by SCORE algorithm.",
+    )
+    score_outliers = File(exists=True, desc="Index of removed volumes, in a CSV file.")
+    mean_cbf_score = File(
+        exists=True,
+        desc="Mean CBF image calculated from SCORE-censored CBF time series.",
+    )
+    mean_cbf_scrub = File(exists=True, desc="average scrub")
 
 
 class ScoreAndScrubCBF(SimpleInterface):
@@ -718,87 +724,88 @@ class ScoreAndScrubCBF(SimpleInterface):
     output_spec = _ScoreAndScrubCBFOutputSpec
 
     def _run_interface(self, runtime):
-        cbf_ts = nb.load(self.inputs.cbf_ts).get_fdata()
-        mask = nb.load(self.inputs.mask).get_fdata()
-        greym = nb.load(self.inputs.gm_tpm).get_fdata()
-        whitem = nb.load(self.inputs.wm_tpm).get_fdata()
-        csf = nb.load(self.inputs.csf_tpm).get_fdata()
-        if cbf_ts.ndim > 3:
-            cbf_scorets, index_score = _score_cbf(
-                cbfts=cbf_ts,
-                wm=whitem,
-                gm=greym,
-                csf=csf,
-                mask=mask,
+        cbf_img = nb.load(self.inputs.cbf)
+        cbf_arr = cbf_img.get_fdata()
+
+        if cbf_arr.ndim > 3:
+            mask_arr = nb.load(self.inputs.brain_mask).get_fdata()
+            gm_tpm_arr = nb.load(self.inputs.gm_tpm).get_fdata()
+            wm_tpm_arr = nb.load(self.inputs.wm_tpm).get_fdata()
+            csf_tpm_arr = nb.load(self.inputs.csf_tpm).get_fdata()
+
+            cbf_ts_score, score_outliers_idx = _score_cbf(
+                cbf_ts=cbf_arr,
+                gm=gm_tpm_arr,
+                wm=wm_tpm_arr,
+                csf=csf_tpm_arr,
+                mask=mask_arr,
                 thresh=self.inputs.tpm_threshold,
             )
-            cbfscrub = _scrub_cbf(
-                cbf_ts=cbf_scorets,
-                gm=greym,
-                wm=whitem,
-                csf=csf,
-                mask=mask,
-                wfun=self.inputs.wavelet_function,
+            mean_cbf_score = np.mean(cbf_ts_score, axis=3)
+            mean_cbf_scrub = _scrub_cbf(
+                cbf_ts=cbf_ts_score,
+                gm=gm_tpm_arr,
+                wm=wm_tpm_arr,
+                csf=csf_tpm_arr,
+                mask=mask_arr,
+                wavelet_function=self.inputs.wavelet_function,
                 thresh=self.inputs.tpm_threshold,
             )
-            mean_cbf_score = np.mean(cbf_scorets, axis=3)
         else:
             config.loggers.interface.warning(
-                f"CBF time series is only {cbf_ts.ndim}D. Skipping SCORE and SCRUB."
+                f"CBF time series is only {cbf_arr.ndim}D. Skipping SCORE and SCRUB."
             )
-            cbf_scorets = cbf_ts
-            index_score = np.array([0])
-            cbfscrub = cbf_ts
-            mean_cbf_score = cbf_ts
+            cbf_ts_score = cbf_arr
+            mean_cbf_score = cbf_arr
+            score_outliers_idx = np.array([0])
+            mean_cbf_scrub = cbf_arr
 
         self._results["cbf_ts_score"] = fname_presuffix(
-            self.inputs.cbf_ts,
+            self.inputs.cbf,
             suffix="_cbfscorets",
             newpath=runtime.cwd,
         )
         self._results["mean_cbf_score"] = fname_presuffix(
-            self.inputs.cbf_ts,
+            self.inputs.cbf,
             suffix="_meancbfscore",
             newpath=runtime.cwd,
         )
         self._results["mean_cbf_scrub"] = fname_presuffix(
-            self.inputs.cbf_ts,
+            self.inputs.cbf,
             suffix="_cbfscrub",
             newpath=runtime.cwd,
         )
-        self._results["score_outlier_index"] = fname_presuffix(
-            self.inputs.cbf_ts,
-            suffix="_scoreindex.txt",
+        self._results["score_outliers"] = fname_presuffix(
+            self.inputs.cbf,
+            suffix="_scoreindex.csv",
             newpath=runtime.cwd,
             use_ext=False,
         )
-        samplecbf = nb.load(self.inputs.mask)
 
         nb.Nifti1Image(
-            dataobj=cbf_scorets,
-            affine=samplecbf.affine,
-            header=samplecbf.header,
+            dataobj=cbf_ts_score,
+            affine=cbf_img.affine,
+            header=cbf_img.header,
         ).to_filename(self._results["cbf_ts_score"])
         nb.Nifti1Image(
             dataobj=mean_cbf_score,
-            affine=samplecbf.affine,
-            header=samplecbf.header,
+            affine=cbf_img.affine,
+            header=cbf_img.header,
         ).to_filename(self._results["mean_cbf_score"])
+        np.savetxt(self._results["score_outliers"], score_outliers_idx, delimiter=",")
         nb.Nifti1Image(
-            dataobj=cbfscrub,
-            affine=samplecbf.affine,
-            header=samplecbf.header,
+            dataobj=mean_cbf_scrub,
+            affine=cbf_img.affine,
+            header=cbf_img.header,
         ).to_filename(self._results["mean_cbf_scrub"])
 
-        np.savetxt(self._results["score_outlier_index"], index_score, delimiter=",")
-
         return runtime
 
 
 class _BASILCBFInputSpec(FSLCommandInputSpec):
     # We use position args here as list indices - so a negative number
     # will put something on the end
-    deltam = File(
+    in_file = File(
         exists=True,
         desc=(
             "ASL data after subtracting tag-control or control-tag. "
@@ -811,11 +818,11 @@ class _BASILCBFInputSpec(FSLCommandInputSpec):
     mask = File(
         exists=True,
         argstr="-m %s",
-        desc="mask in the same space as deltam",
+        desc="mask in the same space as in_file",
         mandatory=True,
     )
     mzero = File(exists=True, argstr="-c %s", desc="m0 scan", mandatory=False)
-    m0_scale = traits.Float(desc="calibration of asl", argstr="--cgain %.2f", mandatory=True)
+    m0scale = traits.Float(desc="calibration of asl", argstr="--cgain %.2f", mandatory=True)
     m0tr = traits.Float(
         desc="The repetition time for the calibration image (the M0 scan).",
         argstr="--tr %.2f",
@@ -860,24 +867,20 @@ class _BASILCBFInputSpec(FSLCommandInputSpec):
         mandatory=True,
         sep=",",
     )
-    slice_spacing = traits.Float(
-        desc="Slice times",
-        argstr="--slicedt %s",
-        mandatory=False,
-    )
     pvc = traits.Bool(
         desc="Do partial volume correction.",
         mandatory=False,
         argstr="--pvcorr",
         default_value=True,
+        usedefault=True,
     )
-    gm_tpm = File(
+    pvgm = File(
         exists=True,
         mandatory=False,
         desc="Partial volume estimates for GM. This is just a GM tissue probability map.",
         argstr="--pvgm %s",
     )
-    wm_tpm = File(
+    pvwm = File(
         exists=True,
         mandatory=False,
         desc="Partial volume estimates for WM. This is just a WM tissue probability map.",
@@ -886,7 +889,7 @@ class _BASILCBFInputSpec(FSLCommandInputSpec):
     alpha = traits.Float(
         desc=(
             "Inversion efficiency - [default: 0.98 (pASL); 0.85 (cASL)]. "
-            "This is equivalent to the BIDS metadata field 'LabelingEfficiency'."
+            "This is equivalent the BIDS metadata field 'LabelingEfficiency'."
         ),
         argstr="--alpha %.2f",
     )
@@ -894,13 +897,22 @@ class _BASILCBFInputSpec(FSLCommandInputSpec):
 
 
 class _BASILCBFOutputSpec(TraitedSpec):
-    mean_cbf_basil = File(exists=True, desc="cbf with spatial correction")
-    mean_cbf_gm_basil = File(exists=True, desc="cbf with spatial correction")
-    mean_cbf_wm_basil = File(
+    mean_cbf = File(exists=True, desc="CBF map in absolute units (ml/100g/min).")
+    mean_cbf_gm = File(
         exists=True,
-        desc="cbf with spatial partial volume white matter correction",
+        desc=(
+            "CBF map with gray matter partial volume correction. "
+            "This means that the map contains 'pure' GM perfusion estimates, in absolute units."
+        ),
+    )
+    mean_cbf_wm = File(
+        exists=True,
+        desc=(
+            "CBF map with white matter partial volume correction. "
+            "This means that the map contains 'pure' WM perfusion estimates, in absolute units."
+        ),
     )
-    att_basil = File(exists=True, desc="arterial transit time")
+    att = File(exists=True, desc="Arterial transit time map.")
 
 
 class BASILCBF(FSLCommand):
@@ -924,8 +936,8 @@ def _run_interface(self, runtime):
         return runtime
 
     def _gen_outfilename(self, suffix):
-        if isdefined(self.inputs.deltam):
-            out_file = self._gen_fname(self.inputs.deltam, suffix=suffix)
+        if isdefined(self.inputs.in_file):
+            out_file = self._gen_fname(self.inputs.in_file, suffix=suffix)
         return os.path.abspath(out_file)
 
     def _list_outputs(self):
@@ -933,13 +945,13 @@ def _list_outputs(self):
 
         outputs = self.output_spec().get()
 
-        outputs["mean_cbf_basil"] = os.path.join(basename, "native_space/perfusion_calib.nii.gz")
-        outputs["att_basil"] = os.path.join(basename, "native_space/arrival.nii.gz")
-        outputs["mean_cbf_gm_basil"] = os.path.join(
+        outputs["mean_cbf"] = os.path.join(basename, "native_space/perfusion_calib.nii.gz")
+        outputs["att"] = os.path.join(basename, "native_space/arrival.nii.gz")
+        outputs["mean_cbf_gm"] = os.path.join(
             basename,
             "native_space/pvcorr/perfusion_calib.nii.gz",
         )
-        outputs["mean_cbf_wm_basil"] = os.path.join(
+        outputs["mean_cbf_wm"] = os.path.join(
             basename,
             "native_space/pvcorr/perfusion_wm_calib.nii.gz",
         )

From 658e1d8269b212eb3fdb0b8f62990c96171cca5c Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Thu, 25 May 2023 15:46:19 -0400
Subject: [PATCH 20/27] Fix connection.

---
 aslprep/interfaces/cbf.py      | 16 +++++++---------
 aslprep/interfaces/plotting.py |  4 ++--
 2 files changed, 9 insertions(+), 11 deletions(-)

diff --git a/aslprep/interfaces/cbf.py b/aslprep/interfaces/cbf.py
index 176750bd5..8d285704d 100644
--- a/aslprep/interfaces/cbf.py
+++ b/aslprep/interfaces/cbf.py
@@ -380,12 +380,10 @@ class _ComputeCBFInputSpec(BaseInterfaceInputSpec):
         ),
     )
     metadata = traits.Dict(
-        exists=True,
         mandatory=True,
         desc="Metadata for the raw CBF file, taken from the raw ASL data's sidecar JSON file.",
     )
     m0_scale = traits.Float(
-        exists=True,
         mandatory=True,
         desc="Relative scale between ASL and M0.",
     )
@@ -672,7 +670,7 @@ def _run_interface(self, runtime):
 
 
 class _ScoreAndScrubCBFInputSpec(BaseInterfaceInputSpec):
-    cbf = File(exists=True, mandatory=True, desc="computed CBF from ComputeCBF")
+    cbf_ts = File(exists=True, mandatory=True, desc="computed CBF from ComputeCBF")
     gm_tpm = File(exists=True, mandatory=True, desc="Gray matter tissue probability map.")
     wm_tpm = File(exists=True, mandatory=True, desc="White matter tissue probability map.")
     csf_tpm = File(exists=True, mandatory=True, desc="Cerebrospinal fluid tissue probability map.")
@@ -684,11 +682,11 @@ class _ScoreAndScrubCBFInputSpec(BaseInterfaceInputSpec):
         desc="Threshold for tissue probability maps.",
     )
     wavelet_function = traits.Str(
-        exists=True,
         mandatory=False,
         default_value="huber",
         option=["bisquare", "andrews", "cauchy", "fair", "logistics", "ols", "talwar", "welsch"],
         desc="Wavelet function used in SCRUB.",
+        usedefault=True,
     )
 
 
@@ -724,7 +722,7 @@ class ScoreAndScrubCBF(SimpleInterface):
     output_spec = _ScoreAndScrubCBFOutputSpec
 
     def _run_interface(self, runtime):
-        cbf_img = nb.load(self.inputs.cbf)
+        cbf_img = nb.load(self.inputs.cbf_ts)
         cbf_arr = cbf_img.get_fdata()
 
         if cbf_arr.ndim > 3:
@@ -761,22 +759,22 @@ def _run_interface(self, runtime):
             mean_cbf_scrub = cbf_arr
 
         self._results["cbf_ts_score"] = fname_presuffix(
-            self.inputs.cbf,
+            self.inputs.cbf_ts,
             suffix="_cbfscorets",
             newpath=runtime.cwd,
         )
         self._results["mean_cbf_score"] = fname_presuffix(
-            self.inputs.cbf,
+            self.inputs.cbf_ts,
             suffix="_meancbfscore",
             newpath=runtime.cwd,
         )
         self._results["mean_cbf_scrub"] = fname_presuffix(
-            self.inputs.cbf,
+            self.inputs.cbf_ts,
             suffix="_cbfscrub",
             newpath=runtime.cwd,
         )
         self._results["score_outliers"] = fname_presuffix(
-            self.inputs.cbf,
+            self.inputs.cbf_ts,
             suffix="_scoreindex.csv",
             newpath=runtime.cwd,
             use_ext=False,
diff --git a/aslprep/interfaces/plotting.py b/aslprep/interfaces/plotting.py
index 17cffa56f..b62ffecf3 100644
--- a/aslprep/interfaces/plotting.py
+++ b/aslprep/interfaces/plotting.py
@@ -104,8 +104,8 @@ def _run_interface(self, runtime):
 
 class _CBFSummaryInputSpec(BaseInterfaceInputSpec):
     cbf = File(exists=True, mandatory=True, desc="")
-    label = traits.Str(exists=True, mandatory=True, desc="label")
-    vmax = traits.Int(exists=True, default_value=90, mandatory=True, desc="max value of asl")
+    label = traits.Str(mandatory=True, desc="label")
+    vmax = traits.Int(default_value=90, mandatory=True, desc="max value of asl")
     ref_vol = File(exists=True, mandatory=True, desc="")
 
 

From 5316259cc0cf8af742b2fa79a87b5781c2b8e1fb Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Thu, 25 May 2023 16:05:38 -0400
Subject: [PATCH 21/27] Update cbf.py

---
 aslprep/interfaces/cbf.py | 38 ++++++++++++++++++++++----------------
 1 file changed, 22 insertions(+), 16 deletions(-)

diff --git a/aslprep/interfaces/cbf.py b/aslprep/interfaces/cbf.py
index 8d285704d..e222fd423 100644
--- a/aslprep/interfaces/cbf.py
+++ b/aslprep/interfaces/cbf.py
@@ -803,7 +803,7 @@ def _run_interface(self, runtime):
 class _BASILCBFInputSpec(FSLCommandInputSpec):
     # We use position args here as list indices - so a negative number
     # will put something on the end
-    in_file = File(
+    deltam = File(
         exists=True,
         desc=(
             "ASL data after subtracting tag-control or control-tag. "
@@ -816,11 +816,11 @@ class _BASILCBFInputSpec(FSLCommandInputSpec):
     mask = File(
         exists=True,
         argstr="-m %s",
-        desc="mask in the same space as in_file",
+        desc="mask in the same space as deltam",
         mandatory=True,
     )
     mzero = File(exists=True, argstr="-c %s", desc="m0 scan", mandatory=False)
-    m0scale = traits.Float(desc="calibration of asl", argstr="--cgain %.2f", mandatory=True)
+    m0_scale = traits.Float(desc="calibration of asl", argstr="--cgain %.2f", mandatory=True)
     m0tr = traits.Float(
         desc="The repetition time for the calibration image (the M0 scan).",
         argstr="--tr %.2f",
@@ -865,6 +865,11 @@ class _BASILCBFInputSpec(FSLCommandInputSpec):
         mandatory=True,
         sep=",",
     )
+    slice_spacing = traits.Float(
+        desc="Slice times",
+        argstr="--slicedt %s",
+        mandatory=False,
+    )
     pvc = traits.Bool(
         desc="Do partial volume correction.",
         mandatory=False,
@@ -872,13 +877,13 @@ class _BASILCBFInputSpec(FSLCommandInputSpec):
         default_value=True,
         usedefault=True,
     )
-    pvgm = File(
+    gm_tpm = File(
         exists=True,
         mandatory=False,
         desc="Partial volume estimates for GM. This is just a GM tissue probability map.",
         argstr="--pvgm %s",
     )
-    pvwm = File(
+    wm_tpm = File(
         exists=True,
         mandatory=False,
         desc="Partial volume estimates for WM. This is just a WM tissue probability map.",
@@ -887,7 +892,7 @@ class _BASILCBFInputSpec(FSLCommandInputSpec):
     alpha = traits.Float(
         desc=(
             "Inversion efficiency - [default: 0.98 (pASL); 0.85 (cASL)]. "
-            "This is equivalent the BIDS metadata field 'LabelingEfficiency'."
+            "This is equivalent to the BIDS metadata field 'LabelingEfficiency'."
         ),
         argstr="--alpha %.2f",
     )
@@ -895,22 +900,22 @@ class _BASILCBFInputSpec(FSLCommandInputSpec):
 
 
 class _BASILCBFOutputSpec(TraitedSpec):
-    mean_cbf = File(exists=True, desc="CBF map in absolute units (ml/100g/min).")
-    mean_cbf_gm = File(
+    mean_cbf_basil = File(exists=True, desc="CBF map in absolute units (ml/100g/min).")
+    mean_cbf_gm_basil = File(
         exists=True,
         desc=(
             "CBF map with gray matter partial volume correction. "
             "This means that the map contains 'pure' GM perfusion estimates, in absolute units."
         ),
     )
-    mean_cbf_wm = File(
+    mean_cbf_wm_basil = File(
         exists=True,
         desc=(
             "CBF map with white matter partial volume correction. "
             "This means that the map contains 'pure' WM perfusion estimates, in absolute units."
         ),
     )
-    att = File(exists=True, desc="Arterial transit time map.")
+    att_basil = File(exists=True, desc="Arterial transit time map.")
 
 
 class BASILCBF(FSLCommand):
@@ -934,8 +939,9 @@ def _run_interface(self, runtime):
         return runtime
 
     def _gen_outfilename(self, suffix):
-        if isdefined(self.inputs.in_file):
-            out_file = self._gen_fname(self.inputs.in_file, suffix=suffix)
+        if isdefined(self.inputs.deltam):
+            out_file = self._gen_fname(self.inputs.deltam, suffix=suffix)
+
         return os.path.abspath(out_file)
 
     def _list_outputs(self):
@@ -943,13 +949,13 @@ def _list_outputs(self):
 
         outputs = self.output_spec().get()
 
-        outputs["mean_cbf"] = os.path.join(basename, "native_space/perfusion_calib.nii.gz")
-        outputs["att"] = os.path.join(basename, "native_space/arrival.nii.gz")
-        outputs["mean_cbf_gm"] = os.path.join(
+        outputs["mean_cbf_basil"] = os.path.join(basename, "native_space/perfusion_calib.nii.gz")
+        outputs["att_basil"] = os.path.join(basename, "native_space/arrival.nii.gz")
+        outputs["mean_cbf_gm_basil"] = os.path.join(
             basename,
             "native_space/pvcorr/perfusion_calib.nii.gz",
         )
-        outputs["mean_cbf_wm"] = os.path.join(
+        outputs["mean_cbf_wm_basil"] = os.path.join(
             basename,
             "native_space/pvcorr/perfusion_wm_calib.nii.gz",
         )

From 643e3521414a158fb1f9f6c975b2a00dc01d0a8f Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Thu, 25 May 2023 16:07:48 -0400
Subject: [PATCH 22/27] Update.

---
 aslprep/interfaces/cbf.py    | 6 +++---
 aslprep/workflows/asl/cbf.py | 4 ++--
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/aslprep/interfaces/cbf.py b/aslprep/interfaces/cbf.py
index e222fd423..6ec4dbd9d 100644
--- a/aslprep/interfaces/cbf.py
+++ b/aslprep/interfaces/cbf.py
@@ -695,7 +695,7 @@ class _ScoreAndScrubCBFOutputSpec(TraitedSpec):
         exists=True,
         desc="CBF time series after removing outlier volumes flagged by SCORE algorithm.",
     )
-    score_outliers = File(exists=True, desc="Index of removed volumes, in a CSV file.")
+    score_outlier_index = File(exists=True, desc="Index of removed volumes, in a CSV file.")
     mean_cbf_score = File(
         exists=True,
         desc="Mean CBF image calculated from SCORE-censored CBF time series.",
@@ -773,7 +773,7 @@ def _run_interface(self, runtime):
             suffix="_cbfscrub",
             newpath=runtime.cwd,
         )
-        self._results["score_outliers"] = fname_presuffix(
+        self._results["score_outlier_index"] = fname_presuffix(
             self.inputs.cbf_ts,
             suffix="_scoreindex.csv",
             newpath=runtime.cwd,
@@ -790,7 +790,7 @@ def _run_interface(self, runtime):
             affine=cbf_img.affine,
             header=cbf_img.header,
         ).to_filename(self._results["mean_cbf_score"])
-        np.savetxt(self._results["score_outliers"], score_outliers_idx, delimiter=",")
+        np.savetxt(self._results["score_outlier_index"], score_outliers_idx, delimiter=",")
         nb.Nifti1Image(
             dataobj=mean_cbf_scrub,
             affine=cbf_img.affine,
diff --git a/aslprep/workflows/asl/cbf.py b/aslprep/workflows/asl/cbf.py
index 84ca9543d..1d331eacf 100644
--- a/aslprep/workflows/asl/cbf.py
+++ b/aslprep/workflows/asl/cbf.py
@@ -381,7 +381,7 @@ def _getfiledir(file):
 """
         # fmt:off
         workflow.connect([
-            (refine_mask, score_and_scrub_cbf, [("out_mask", "mask")]),
+            (refine_mask, score_and_scrub_cbf, [("out_mask", "brain_mask")]),
             (compute_cbf, score_and_scrub_cbf, [("cbf_ts", "cbf_ts")]),
             (gm_tfm, score_and_scrub_cbf, [("output_image", "gm_tpm")]),
             (wm_tfm, score_and_scrub_cbf, [("output_image", "wm_tpm")]),
@@ -774,7 +774,7 @@ def _getfiledir(file):
 
         # fmt:off
         workflow.connect([
-            (refine_mask, score_and_scrub_cbf, [("out_mask", "mask")]),
+            (refine_mask, score_and_scrub_cbf, [("out_mask", "brain_mask")]),
             (gm_tfm, score_and_scrub_cbf, [("output_image", "gm_tpm")]),
             (wm_tfm, score_and_scrub_cbf, [("output_image", "wm_tpm")]),
             (csf_tfm, score_and_scrub_cbf, [("output_image", "csf_tpm")]),

From 6e1456bb77dff6426298c986e34aab1d4d7ae5be Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Fri, 26 May 2023 12:20:01 -0400
Subject: [PATCH 23/27] Update cbf.py

---
 aslprep/utils/cbf.py | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py
index 5f8fdb6f1..750b86073 100644
--- a/aslprep/utils/cbf.py
+++ b/aslprep/utils/cbf.py
@@ -3,7 +3,7 @@
 
 import numpy as np
 from scipy.interpolate import interp1d
-from scipy.stats import median_abs_deviation
+from scipy.stats import chi2, median_abs_deviation
 
 from aslprep import config
 
@@ -440,6 +440,11 @@ def _scrub_cbf(cbf_ts, gm, wm, csf, mask, wavelet_function="huber", thresh=0.7):
     mean_cbf = np.mean(cbf_ts, axis=3)
     masked_cbf_ts = cbf_ts[mask, :].T  # TxS array
     masked_var_map = np.var(masked_cbf_ts, axis=0)  # S array
+    n_voxels = masked_cbf_ts.shape[1]
+    # "Since μ was designed based on whole brain, we designed the algorithm to only penalize values
+    # whose normalized temporal variance is outside 99.99 percentile value (p99.99) of the chi
+    # square distribution computed with (Number of brain voxels – 1) degrees of freedom."
+    p_thresh_0001 = chi2.ppf(0.0001, df=n_voxels - 1)
     thresh1, thresh3 = _getchisquare(n_volumes)
     mu1 = masked_var_map / (np.median(masked_var_map[gm_idx == 1]) * thresh3)
     mu = (

From 4fa33baaf3ca0d3ebe378f50e3173cee8d867245 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Fri, 26 May 2023 15:10:27 -0400
Subject: [PATCH 24/27] Rename wavelet_function to cost_function.

---
 aslprep/interfaces/cbf.py    | 4 ++--
 aslprep/utils/cbf.py         | 6 +++---
 aslprep/workflows/asl/cbf.py | 4 ++--
 3 files changed, 7 insertions(+), 7 deletions(-)

diff --git a/aslprep/interfaces/cbf.py b/aslprep/interfaces/cbf.py
index 6ec4dbd9d..aced3b99d 100644
--- a/aslprep/interfaces/cbf.py
+++ b/aslprep/interfaces/cbf.py
@@ -681,7 +681,7 @@ class _ScoreAndScrubCBFInputSpec(BaseInterfaceInputSpec):
         usedefault=True,
         desc="Threshold for tissue probability maps.",
     )
-    wavelet_function = traits.Str(
+    cost_function = traits.Str(
         mandatory=False,
         default_value="huber",
         option=["bisquare", "andrews", "cauchy", "fair", "logistics", "ols", "talwar", "welsch"],
@@ -746,7 +746,7 @@ def _run_interface(self, runtime):
                 wm=wm_tpm_arr,
                 csf=csf_tpm_arr,
                 mask=mask_arr,
-                wavelet_function=self.inputs.wavelet_function,
+                cost_function=self.inputs.cost_function,
                 thresh=self.inputs.tpm_threshold,
             )
         else:
diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py
index 750b86073..a2789800f 100644
--- a/aslprep/utils/cbf.py
+++ b/aslprep/utils/cbf.py
@@ -403,7 +403,7 @@ def _robust_fit(
     return b
 
 
-def _scrub_cbf(cbf_ts, gm, wm, csf, mask, wavelet_function="huber", thresh=0.7):
+def _scrub_cbf(cbf_ts, gm, wm, csf, mask, cost_function="huber", thresh=0.7):
     """Apply SCRUB algorithm to CBF data.
 
     Parameters
@@ -477,7 +477,7 @@ def _scrub_cbf(cbf_ts, gm, wm, csf, mask, wavelet_function="huber", thresh=0.7):
     global_prior_full = c[0] * gm.flatten() + c[1] * wm.flatten()
     global_prior = global_prior_full[mask.flatten() == 1]
 
-    tune = _get_wfun_tuner(wfun=wavelet_function)
+    tune = _get_wfun_tuner(wfun=cost_function)
     masked_mean_cbf_scrub = _robust_fit(
         Y=masked_cbf_ts,
         mu=mu,
@@ -485,7 +485,7 @@ def _scrub_cbf(cbf_ts, gm, wm, csf, mask, wavelet_function="huber", thresh=0.7):
         modrobprior=modrobprior,
         lmd=0,
         localprior=0,
-        wfun=wavelet_function,
+        wfun=cost_function,
         tune=tune,
         flagstd=1,
         flagmodrobust=1,
diff --git a/aslprep/workflows/asl/cbf.py b/aslprep/workflows/asl/cbf.py
index 1d331eacf..e042c5874 100644
--- a/aslprep/workflows/asl/cbf.py
+++ b/aslprep/workflows/asl/cbf.py
@@ -365,7 +365,7 @@ def _getfiledir(file):
 
     if scorescrub:
         score_and_scrub_cbf = pe.Node(
-            ScoreAndScrubCBF(tpm_threshold=0.7, wavelet_function="huber"),
+            ScoreAndScrubCBF(tpm_threshold=0.7, cost_function="huber"),
             mem_gb=0.2,
             name="score_and_scrub_cbf",
             run_without_submitting=True,
@@ -766,7 +766,7 @@ def _getfiledir(file):
 """
 
         score_and_scrub_cbf = pe.Node(
-            ScoreAndScrubCBF(tpm_threshold=0.7, wavelet_function="huber"),
+            ScoreAndScrubCBF(tpm_threshold=0.7, cost_function="huber"),
             mem_gb=mem_gb,
             name="score_and_scrub_cbf",
             run_without_submitting=True,

From 410862500c9d6a9422aaccec0b71160c61f3b6fc Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Fri, 26 May 2023 15:11:29 -0400
Subject: [PATCH 25/27] Keep working on SCRUB.

---
 aslprep/utils/cbf.py | 35 +++++++++++++++++++++++++++++++----
 1 file changed, 31 insertions(+), 4 deletions(-)

diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py
index a2789800f..e01564fc5 100644
--- a/aslprep/utils/cbf.py
+++ b/aslprep/utils/cbf.py
@@ -406,6 +406,27 @@ def _robust_fit(
 def _scrub_cbf(cbf_ts, gm, wm, csf, mask, cost_function="huber", thresh=0.7):
     """Apply SCRUB algorithm to CBF data.
 
+    The SCRUB algorithm uses a Bayesian approach to compute the average CBF over a series of
+    CBF estimates (volumes), while taking into account a prior distribution of CBF values
+    (in the Bayesian sense).
+
+    SCRUB assumes that both the likelihood and prior term in the formula are Gaussian,
+    which means that the tuning parameter (lambda) is equal to the ratio of the variances of the
+    additive noise and the prior probability density.
+    A higher lambda value generally indicates less reliable CBF estimates,
+    driven by higher variance of the additive noise.
+    Basically, temporal variance of the CBF estimates is used as an indicator of unreliable
+    voxels/estimates.
+
+    SCRUB may account for other factors indicating unreliability by adjusting lambda,
+    such as voxels with negative CBF estimates.
+
+    The distribution of CBF is not Gaussian due to frequent presence of outliers.
+    SCRUB adds a data fidelity term to account for large outliers (rho?).
+
+    The prior is mu. The tuning parameter is lambda (consistent across voxels?).
+    The cost function is rho.
+
     Parameters
     ----------
     cbf_ts : 4D numpy.ndarray
@@ -414,8 +435,9 @@ def _scrub_cbf(cbf_ts, gm, wm, csf, mask, cost_function="huber", thresh=0.7):
         Gray matter, white matter, and CSF tissue probability maps.
     mask : 3D numpy.ndarray
         Binary brain mask.
-    wavelet_function : str
-        Wavelet function
+    cost_function : str
+        Wavelet function. This is the data fidelity term, which determines the cost of deviations
+        of CBF estimates from the solution.
     thresh : float
         Tissue probability threshold.
         Voxels in gm, wm, and csf with values greater than or equal to the threshold are
@@ -439,13 +461,17 @@ def _scrub_cbf(cbf_ts, gm, wm, csf, mask, cost_function="huber", thresh=0.7):
 
     mean_cbf = np.mean(cbf_ts, axis=3)
     masked_cbf_ts = cbf_ts[mask, :].T  # TxS array
-    masked_var_map = np.var(masked_cbf_ts, axis=0)  # S array
     n_voxels = masked_cbf_ts.shape[1]
     # "Since μ was designed based on whole brain, we designed the algorithm to only penalize values
     # whose normalized temporal variance is outside 99.99 percentile value (p99.99) of the chi
     # square distribution computed with (Number of brain voxels – 1) degrees of freedom."
-    p_thresh_0001 = chi2.ppf(0.0001, df=n_voxels - 1)
+    # (n - 1) * Sample variance (S squared) / Variance (sigma squared) follows chi-squared
+    # distribution with n - 1 degrees of freedom.
+    p_thresh_9999 = chi2.ppf(0.0001, df=n_voxels - 1)
     thresh1, thresh3 = _getchisquare(n_volumes)
+
+    # Construct voxel-wise prior map (i.e., mu_v)
+    masked_var_map = np.var(masked_cbf_ts, axis=0)  # S array of sigma squared
     mu1 = masked_var_map / (np.median(masked_var_map[gm_idx == 1]) * thresh3)
     mu = (
         ((mu1 > thresh1) & (mu1 < 10 * thresh1)) * (mu1 - thresh1)
@@ -456,6 +482,7 @@ def _scrub_cbf(cbf_ts, gm, wm, csf, mask, cost_function="huber", thresh=0.7):
     mu_map = np.zeros_like(mean_cbf)
     mu_map[mask] = mu
 
+    # Now I guess we work on lambda?
     modrobprior = mu / 10
     gm_idx2 = (
         [gm.flatten() > thresh] and [mu_map.flatten() == 0] and [wm.flatten() > csf.flatten()]

From d099a5ecf7f43701706d5eaf7096daa7824160cc Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Tue, 12 Dec 2023 16:57:29 -0500
Subject: [PATCH 26/27] Address style issues.

---
 aslprep/interfaces/cbf.py | 2 +-
 aslprep/utils/cbf.py      | 7 +++++--
 2 files changed, 6 insertions(+), 3 deletions(-)

diff --git a/aslprep/interfaces/cbf.py b/aslprep/interfaces/cbf.py
index 8e9d91ff5..1bbc23c61 100644
--- a/aslprep/interfaces/cbf.py
+++ b/aslprep/interfaces/cbf.py
@@ -736,7 +736,7 @@ def _run_interface(self, runtime):
             header=cbf_img.header,
         ).to_filename(self._results["mean_cbf_scrub"])
 
-        score_outlier_df = pd.DataFrame(columns=["score_outlier_index"], data=index_score)
+        score_outlier_df = pd.DataFrame(columns=["score_outlier_index"], data=score_outliers_idx)
         score_outlier_df.to_csv(self._results["score_outlier_index"], sep="\t", index=False)
 
         return runtime
diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py
index 878af76f3..adc472707 100644
--- a/aslprep/utils/cbf.py
+++ b/aslprep/utils/cbf.py
@@ -461,15 +461,18 @@ def _scrub_cbf(cbf_ts, gm, wm, csf, mask, cost_function="huber", thresh=0.7):
 
     mean_cbf = np.mean(cbf_ts, axis=3)
     masked_cbf_ts = cbf_ts[mask, :].T  # TxS array
-    n_voxels = masked_cbf_ts.shape[1]
     # "Since μ was designed based on whole brain, we designed the algorithm to only penalize values
     # whose normalized temporal variance is outside 99.99 percentile value (p99.99) of the chi
     # square distribution computed with (Number of brain voxels – 1) degrees of freedom."
     # (n - 1) * Sample variance (S squared) / Variance (sigma squared) follows chi-squared
     # distribution with n - 1 degrees of freedom.
-    p_thresh_9999 = chi2.ppf(0.0001, df=n_voxels - 1)
     thresh1, thresh3 = _getchisquare(n_volumes)
 
+    # XXX: It seems like the following should produce the same threshold as _getchisquare,
+    # but it doesn't. I haven't been able to figure out why.
+    # n_voxels = masked_cbf_ts.shape[1]
+    # p_thresh_9999 = chi2.ppf(0.0001, df=n_voxels - 1)
+
     # Construct voxel-wise prior map (i.e., mu_v)
     masked_var_map = np.var(masked_cbf_ts, axis=0)  # S array of sigma squared
     mu1 = masked_var_map / (np.median(masked_var_map[gm_idx == 1]) * thresh3)

From ac2356baf81e94d78cb4b2a7d3332e6ba4153343 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Tue, 12 Dec 2023 16:59:35 -0500
Subject: [PATCH 27/27] Remove unused import.

---
 aslprep/utils/cbf.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py
index adc472707..bf5772491 100644
--- a/aslprep/utils/cbf.py
+++ b/aslprep/utils/cbf.py
@@ -3,7 +3,7 @@
 
 import numpy as np
 from scipy.interpolate import interp1d
-from scipy.stats import chi2, median_abs_deviation
+from scipy.stats import median_abs_deviation
 
 from aslprep import config
 
@@ -470,6 +470,7 @@ def _scrub_cbf(cbf_ts, gm, wm, csf, mask, cost_function="huber", thresh=0.7):
 
     # XXX: It seems like the following should produce the same threshold as _getchisquare,
     # but it doesn't. I haven't been able to figure out why.
+    # from scipy.stats import chi2
     # n_voxels = masked_cbf_ts.shape[1]
     # p_thresh_9999 = chi2.ppf(0.0001, df=n_voxels - 1)