confounds.py

# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
"""Workflows for calculating confounds for ASL data."""
from nipype.algorithms import confounds as nac
from nipype.interfaces import utility as niu
from nipype.pipeline import engine as pe
from templateflow.api import get as get_template

from aslprep.config import DEFAULT_MEMORY_MIN_GB
from aslprep.interfaces import ASLSummary, DerivativesDataSink, GatherConfounds
from aslprep.interfaces.ants import ApplyTransforms
from aslprep.niworkflows.engine.workflows import LiterateWorkflow as Workflow
from aslprep.niworkflows.interfaces.utils import AddTSVHeader


def init_asl_confounds_wf(
    mem_gb,
    name="asl_confounds_wf",
):
    """Build a workflow to generate and write out confounding signals.

    This workflow calculates confounds for a asl series, and aggregates them
    into a :abbr:`TSV (tab-separated value)` file, for use as nuisance
    regressors in a :abbr:`GLM (general linear model)`.
    The following confounds are calculated, with column headings in parentheses:

    #. DVARS - original and standardized variants (``dvars``, ``std_dvars``)
    #. Framewise displacement, based on head-motion parameters
       (``framewise_displacement``)
    #. Estimated head-motion parameters, in mm and rad
       (``trans_x``, ``trans_y``, ``trans_z``, ``rot_x``, ``rot_y``, ``rot_z``)

    Workflow Graph
        .. workflow::
            :graph2use: orig
            :simple_form: yes

            from aslprep.workflows.asl.confounds import init_asl_confounds_wf

            wf = init_asl_confounds_wf(
                mem_gb=1,
            )

    Parameters
    ----------
    mem_gb : :obj:`float`
        Size of asl file in GB - please note that this size
        should be calculated after resamplings that may extend
        the FoV
    metadata : :obj:`dict`
        BIDS metadata for asl file
    name : :obj:`str`
        Name of workflow (default: ``asl_confounds_wf``)


    Inputs
    ------
    asl
        asl image, after the prescribed corrections (HMC and SDC)
        when available.
    asl_mask
        asl series mask
    movpar_file
        SPM-formatted motion parameters file
    rmsd_file
        Framewise displacement as measured by ``fsl_motion_outliers``.
    skip_vols
        number of non steady state volumes
    t1w_mask
        Mask of the skull-stripped template image
    t1w_tpms
        List of tissue probability maps in T1w space
    anat_to_aslref_xfm
        Affine matrix that maps the T1w space into alignment with
        the native asl space

    Outputs
    -------
    confounds_file
        TSV of all aggregated confounds
    confounds_metadata
        Confounds metadata dictionary.
    """
    workflow = Workflow(name=name)
    workflow.__desc__ = """\
Several confounding timeseries were calculated, including both framewise displacement
(FD) and DVARS. FD and DVARS are calculated using the implementations in in *Nipype*
(following the definition by [@power_fd_dvars]) for each ASL run.  ASLPrep summarizes
in-scanner motion as the mean framewise displacement and relative root-mean square displacement.
"""
    inputnode = pe.Node(
        niu.IdentityInterface(
            fields=[
                "asl",
                "asl_mask",
                "movpar_file",
                "rmsd_file",
                "skip_vols",
                "t1w_mask",
                "t1w_tpms",
                "anat_to_aslref_xfm",
            ]
        ),
        name="inputnode",
    )
    outputnode = pe.Node(
        niu.IdentityInterface(fields=["confounds_file", "confounds_metadata"]), name="outputnode"
    )

    # DVARS
    dvars = pe.Node(
        nac.ComputeDVARS(save_nstd=True, save_std=True, remove_zerovariance=True),
        name="dvars",
        mem_gb=mem_gb,
    )

    # Frame displacement
    fdisp = pe.Node(nac.FramewiseDisplacement(parameter_source="SPM"), name="fdisp", mem_gb=mem_gb)

    # Global and segment regressors
    # signals_class_labels = ["csf", "white_matter", "global_signal"]

    # Arrange confounds
    add_dvars_header = pe.Node(
        AddTSVHeader(columns=["dvars"]),
        name="add_dvars_header",
        mem_gb=0.01,
        run_without_submitting=True,
    )
    add_std_dvars_header = pe.Node(
        AddTSVHeader(columns=["std_dvars"]),
        name="add_std_dvars_header",
        mem_gb=0.01,
        run_without_submitting=True,
    )
    add_motion_headers = pe.Node(
        AddTSVHeader(columns=["trans_x", "trans_y", "trans_z", "rot_x", "rot_y", "rot_z"]),
        name="add_motion_headers",
        mem_gb=0.01,
        run_without_submitting=True,
    )
    add_rmsd_header = pe.Node(
        AddTSVHeader(columns=["rmsd"]),
        name="add_rmsd_header",
        mem_gb=0.01,
        run_without_submitting=True,
    )
    concat = pe.Node(GatherConfounds(), name="concat", mem_gb=0.01, run_without_submitting=True)

    # Expand model to include derivatives and quadratics
    # fmt:off
    workflow.connect([
        # Connect inputnode to each non-anatomical confound node
        (inputnode, dvars, [
            ("asl", "in_file"),
            ("asl_mask", "in_mask"),
        ]),
        (inputnode, fdisp, [("movpar_file", "in_file")]),
        # Collate computed confounds together
        (inputnode, add_motion_headers, [("movpar_file", "in_file")]),
        (inputnode, add_rmsd_header, [("rmsd_file", "in_file")]),
        (dvars, add_dvars_header, [("out_nstd", "in_file")]),
        (dvars, add_std_dvars_header, [("out_std", "in_file")]),
        (fdisp, concat, [("out_file", "fd")]),
        (add_motion_headers, concat, [("out_file", "motion")]),
        (add_rmsd_header, concat, [("out_file", "rmsd")]),
        (add_dvars_header, concat, [("out_file", "dvars")]),
        (add_std_dvars_header, concat, [("out_file", "std_dvars")]),
        # Set outputs
        (concat, outputnode, [("confounds_file", "confounds_file")]),
    ])
    # fmt:on

    return workflow


def init_carpetplot_wf(mem_gb, metadata, name="carpetplot_wf"):
    """Build a workflow to generate carpet plots.

    Resamples the MNI parcellation (ad-hoc parcellation derived from the
    Harvard-Oxford template and others).

    Parameters
    ----------
    mem_gb : :obj:`float`
        Size of ASL file in GB - please note that this size
        should be calculated after resamplings that may extend
        the FoV
    metadata : :obj:`dict`
        BIDS metadata for ASL file
    name : :obj:`str`
        Name of workflow (default: ``carpetplot_wf``)

    Inputs
    ------
    asl
        asl image, after the prescribed corrections (HMC and SDC)
        when available.
    asl_mask
        ASL series mask
    confounds_file
        TSV of all aggregated confounds
    anat_to_aslref_xfm
        Affine matrix that maps the T1w space into alignment with
        the native ASL space
    template_to_anat_xfm
        ANTs-compatible affine-and-warp transform file
    """
    workflow = Workflow(name=name)

    inputnode = pe.Node(
        niu.IdentityInterface(
            fields=[
                "asl",
                "asl_mask",
                "confounds_file",
                "anat_to_aslref_xfm",
                "template_to_anat_xfm",
            ]
        ),
        name="inputnode",
    )

    # List transforms
    mrg_xfms = pe.Node(niu.Merge(2), name="mrg_xfms")

    # Warp segmentation into EPI space
    resample_parc = pe.Node(
        ApplyTransforms(
            float=True,
            input_image=str(
                get_template(
                    "MNI152NLin2009cAsym",
                    resolution=1,
                    desc="carpet",
                    suffix="dseg",
                    extension=[".nii", ".nii.gz"],
                )
            ),
            dimension=3,
            default_value=0,
            interpolation="MultiLabel",
        ),
        name="resample_parc",
    )

    # Carpetplot and confounds plot
    conf_plot = pe.Node(
        ASLSummary(
            tr=metadata.get("RepetitionTime", metadata["RepetitionTimePreparation"]),
            confounds_list=[("std_dvars", None, "DVARS"), ("framewise_displacement", "mm", "FD")],
        ),
        name="conf_plot",
        mem_gb=mem_gb,
    )
    ds_report_asl_conf = pe.Node(
        DerivativesDataSink(desc="carpetplot", datatype="figures", keep_dtype=True),
        name="ds_report_asl_conf",
        run_without_submitting=True,
        mem_gb=DEFAULT_MEMORY_MIN_GB,
    )

    # fmt:off
    workflow.connect([
        (inputnode, mrg_xfms, [
            ("anat_to_aslref_xfm", "in1"),
            ("template_to_anat_xfm", "in2"),
        ]),
        (inputnode, resample_parc, [("asl_mask", "reference_image")]),
        (mrg_xfms, resample_parc, [("out", "transforms")]),
        # Carpetplot
        (inputnode, conf_plot, [
            ("asl", "in_func"),
            ("asl_mask", "in_mask"),
            ("confounds_file", "confounds_file"),
        ]),
        (resample_parc, conf_plot, [("output_image", "in_segm")]),
        (conf_plot, ds_report_asl_conf, [("out_file", "in_file")]),
    ])
    # fmt:on

    return workflow