Support ALFF for censored data (#1020)

docs/outputs.rst

@@ -244,8 +244,7 @@ Resting-state metric derivatives (ReHo and ALFF)
 
 .. important::
       ALFF will not be generated if bandpass filtering is disabled with the
-      ``--disable-bandpass-filtering`` parameter,
-      or if high-motion outlier censoring is enabled ``--fd-thresh`` is greater than zero.
+      ``--disable-bandpass-filtering`` parameter.
 
 *XCP-D* will also parcellate the ReHo and ALFF maps with each of the atlases used for the BOLD
 data.

docs/workflows.rst

@@ -387,10 +387,23 @@ ALFF
 ----
 :func:`~xcp_d.workflows.restingstate.init_alff_wf`
 
+Amplitude of low-frequency fluctuation (ALFF) is a measure that ostensibly localizes
+spontaneous neural activity in resting-state BOLD data.
+It is calculated by the following:
+
+1. The ``filtered, interpolated, denoised BOLD`` is passed along to the ALFF workflow.
+2. Voxel-wise BOLD time series are normalized (mean-centered and scaled to unit standard deviation)
+   over time.
+3. The power spectrum and associated frequencies are estimated from the BOLD data.
+   -  If censoring+interpolation was not performed, then this uses :func:`scipy.signal.periodogram`.
+   -  If censoring+interpolation was performed, then this uses :func:`scipy.signal.lombscargle`.
+4. The square root of the power spectrum is calculated.
+5. The power spectrum values corresponding to the frequency range retained by the
+   temporal filtering step are extracted from the full power spectrum.
+6. The mean of the within-band power spectrum is calculated and multiplied by 2.
+
 ALFF will only be calculated if the bandpass filter is enabled
-(i.e., if the ``--disable-bandpass-filter`` flag is not used)
-and censoring is disabled
-(i.e., if ``--fd-thresh`` is set to a value less than or equal to zero).
+(i.e., if the ``--disable-bandpass-filter`` flag is not used).
 
 Smoothed ALFF derivatives will also be generated if the ``--smoothing`` flag is used.
 

xcp_d/data/boilerplate.bib

@@ -724,3 +724,53 @@ @article{najdenovska2018vivo
   url={https://doi.org/10.1038/sdata.2018.270},
   doi={10.1038/sdata.2018.270}
 }
+
+@article{taylorlomb,
+  title={Lomb-Scargle your way to RSFC parameter estimation in AFNI-FATCAT},
+  author={Taylor, Paul A and Chen, Gang and Glen, Daniel R and Reynolds, Richard C and Cox, Robert W},
+  year={2018},
+  organization={International Society for Magnetic Resonance in Medicine},
+  journal={International Society for Magnetic Resonance in Medicine},
+  note={Description of Lomb-Scargle applied to scrubbed fMRI data in AFNI.}
+}
+
+@article{lomb1976least,
+  title={Least-squares frequency analysis of unequally spaced data},
+  author={Lomb, Nicholas R},
+  journal={Astrophysics and space science},
+  volume={39},
+  pages={447--462},
+  year={1976},
+  publisher={Springer}
+}
+
+@article{scargle1982studies,
+  title={Studies in astronomical time series analysis. II-Statistical aspects of spectral analysis of unevenly spaced data},
+  author={Scargle, Jeffrey D},
+  journal={Astrophysical Journal, Part 1, vol. 263, Dec. 15, 1982, p. 835-853.},
+  volume={263},
+  pages={835--853},
+  year={1982}
+}
+
+@article{townsend2010fast,
+  title={Fast calculation of the Lomb--Scargle periodogram using graphics processing units},
+  author={Townsend, RHD},
+  journal={The Astrophysical Journal Supplement Series},
+  volume={191},
+  number={2},
+  pages={247},
+  year={2010},
+  publisher={IOP Publishing}
+}
+
+@article{yu2007altered,
+  title={Altered baseline brain activity in children with ADHD revealed by resting-state functional MRI},
+  author={Yu-Feng, Zang and Yong, He and Chao-Zhe, Zhu and Qing-Jiu, Cao and Man-Qiu, Sui and Meng, Liang and Li-Xia, Tian and Tian-Zi, Jiang and Yu-Feng, Wang},
+  journal={Brain and Development},
+  volume={29},
+  number={2},
+  pages={83--91},
+  year={2007},
+  publisher={Elsevier}
+}

xcp_d/interfaces/restingstate.py

@@ -8,6 +8,7 @@
 import os
 import shutil
 
+import pandas as pd
 from nipype import logging
 from nipype.interfaces.afni.preprocess import Despike, DespikeInputSpec
 from nipype.interfaces.afni.utils import ReHoInputSpec, ReHoOutputSpec
@@ -16,6 +17,7 @@
     File,
     SimpleInterface,
     TraitedSpec,
+    Undefined,
     traits,
     traits_extension,
 )
@@ -88,44 +90,70 @@ class _ComputeALFFInputSpec(BaseInterfaceInputSpec):
     TR = traits.Float(mandatory=True, desc="repetition time")
     low_pass = traits.Float(
         mandatory=True,
-        default_value=0.10,
         desc="low_pass filter in Hz",
     )
     high_pass = traits.Float(
         mandatory=True,
-        default_value=0.01,
         desc="high_pass filter in Hz",
     )
     mask = File(
         exists=True,
         mandatory=False,
         desc=" brain mask for nifti file",
     )
+    temporal_mask = traits.Either(
+        File(exists=True),
+        Undefined,
+        mandatory=False,
+        desc="Temporal mask.",
+    )
 
 
 class _ComputeALFFOutputSpec(TraitedSpec):
     alff = File(exists=True, mandatory=True, desc=" alff")
 
 
 class ComputeALFF(SimpleInterface):
-    """Compute ALFF."""
+    """Compute amplitude of low-frequency fluctuation (ALFF).
+
+    Notes
+    -----
+    The ALFF implementation is based on :footcite:t:`yu2007altered`,
+    although the ALFF values are not scaled by the mean ALFF value across the brain.
+
+    If censoring is applied (i.e., ``fd_thresh > 0``), then the power spectrum will be estimated
+    using a Lomb-Scargle periodogram
+    :footcite:p:`lomb1976least,scargle1982studies,townsend2010fast,taylorlomb`.
+
+    References
+    ----------
+    .. footbibliography::
+    """
 
     input_spec = _ComputeALFFInputSpec
     output_spec = _ComputeALFFOutputSpec
 
     def _run_interface(self, runtime):
         # Get the nifti/cifti into matrix form
         data_matrix = read_ndata(datafile=self.inputs.in_file, maskfile=self.inputs.mask)
+
+        sample_mask = None
+        temporal_mask = self.inputs.temporal_mask
+        if isinstance(temporal_mask, str) and os.path.isfile(temporal_mask):
+            censoring_df = pd.read_table(temporal_mask)
+            # Invert the temporal mask to make retained volumes 1s and dropped volumes 0s.
+            sample_mask = ~censoring_df["framewise_displacement"].values.astype(bool)
+
         # compute the ALFF
         alff_mat = compute_alff(
             data_matrix=data_matrix,
             low_pass=self.inputs.low_pass,
             high_pass=self.inputs.high_pass,
             TR=self.inputs.TR,
+            sample_mask=sample_mask,
         )
 
         # Write out the data
-
         if self.inputs.in_file.endswith(".dtseries.nii"):
             suffix = "_alff.dscalar.nii"
         elif self.inputs.in_file.endswith(".nii.gz"):