[ENH] Derive masks to support native-space data

aroma/aroma.py

@@ -22,6 +22,7 @@ def aroma_workflow(
     TR=None,
     overwrite=False,
     generate_plots=True,
+    csf=None,
 ):
     """Run the AROMA workflow.
 
@@ -151,53 +152,27 @@ def aroma_workflow(
 
     # Define/create mask. Either by making a copy of the specified mask, or by
     # creating a new one.
-    new_mask = op.join(out_dir, "mask.nii.gz")
-    if mask:
-        shutil.copyfile(mask, new_mask)
-    elif in_feat and op.isfile(op.join(in_feat, "example_func.nii.gz")):
-        # If a Feat directory is specified, and an example_func is present use
-        # example_func to create a mask
-        bet_command = "{0} {1} {2} -f 0.3 -n -m -R".format(
-            op.join(fsl_dir, "bet"),
-            op.join(in_feat, "example_func.nii.gz"),
-            op.join(out_dir, "bet"),
-        )
-        os.system(bet_command)
-        os.rename(op.join(out_dir, "bet_mask.nii.gz"), new_mask)
-        if op.isfile(op.join(out_dir, "bet.nii.gz")):
-            os.remove(op.join(out_dir, "bet.nii.gz"))
-    else:
-        if in_feat:
-            print(
-                " - No example_func was found in the Feat directory. "
-                "A mask will be created including all voxels with varying "
-                "intensity over time in the fMRI data. Please check!\n"
-            )
-        math_command = "{0} {1} -Tstd -bin {2}".format(
-            op.join(fsl_dir, "fslmaths"), in_file, new_mask
-        )
-        os.system(math_command)
+    masks = utils.derive_masks(in_file, csf=None)
 
     # Run ICA-AROMA
     print("Step 1) MELODIC")
-    utils.runICA(fsl_dir, in_file, out_dir, mel_dir, new_mask, dim, TR)
+    component_maps, mixing, mixing_FT = utils.runICA(
+        fsl_dir, in_file, out_dir, mel_dir, masks["brain"], dim, TR
+    )
 
     print("Step 2) Automatic classification of the components")
     print("  - registering the spatial maps to MNI")
-    mel_IC = op.join(out_dir, "melodic_IC_thr.nii.gz")
     mel_IC_MNI = op.join(out_dir, "melodic_IC_thr_MNI2mm.nii.gz")
-    utils.register2MNI(fsl_dir, mel_IC, mel_IC_MNI, affmat, warp)
+    utils.register2MNI(fsl_dir, component_maps, mel_IC_MNI, affmat, warp)
 
     print("  - extracting the CSF & Edge fraction features")
-    edge_fract, csf_fract = features.feature_spatial(mel_IC_MNI)
+    edge_fract, csf_fract = features.feature_spatial(mel_IC_MNI, masks)
 
     print("  - extracting the Maximum RP correlation feature")
-    mel_mix = op.join(out_dir, "melodic.ica", "melodic_mix")
-    max_RP_corr = features.feature_time_series(mel_mix, mc)
+    max_RP_corr = features.feature_time_series(mixing, mc)
 
     print("  - extracting the High-frequency content feature")
-    mel_FT_mix = op.join(out_dir, "melodic.ica", "melodic_FTmix")
-    HFC = features.feature_frequency(mel_FT_mix, TR)
+    HFC = features.feature_frequency(mixing_FT, TR)
 
     print("  - classification")
     motion_ICs = utils.classification(out_dir, max_RP_corr, edge_fract, HFC, csf_fract)
@@ -211,6 +186,6 @@ def aroma_workflow(
 
     if den_type != "no":
         print("Step 3) Data denoising")
-        utils.denoising(fsl_dir, in_file, out_dir, mel_mix, den_type, motion_ICs)
+        utils.denoising(fsl_dir, in_file, out_dir, mixing, den_type, motion_ICs)
 
     print("Finished")

aroma/cli/aroma.py

@@ -38,7 +38,11 @@ def _get_parser():
         dest="in_file",
         type=lambda x: is_valid_file(parser, x),
         required=False,
-        help="Input file name of fMRI data (.nii.gz)",
+        help=(
+            "Input file name of fMRI data (.nii.gz). "
+            "This file may be in standard (MNI152) or native space, with some restrictions. "
+            "The data should be smoothed prior to running AROMA."
+        ),
     )
     inputs.add_argument(
         "-f",
@@ -47,8 +51,9 @@ def _get_parser():
         required=False,
         type=lambda x: is_valid_path(parser, x),
         help=(
-            "Feat directory name (Feat should have been run without temporal "
-            "filtering and including registration to MNI152)"
+            "Path to FSL FEAT directory. "
+            "FEAT should have been run without temporal filtering, "
+            "but with registration to MNI152 space."
         ),
     )
 
@@ -113,6 +118,17 @@ def _get_parser():
     # Optional options
     optoptions = parser.add_argument_group("Optional arguments")
     optoptions.add_argument("-tr", dest="TR", help="TR in seconds", type=float)
+    optoptions.add_argument(
+        "--csf",
+        dest="csf",
+        type=lambda x: is_valid_file(parser, x),
+        default=None,
+        help=(
+            "Path to a cerebrospinal fluid (CSF) mask or tissue probability map. "
+            "If this file is not provided, then data are assumed to be in standard space, "
+            "and prepackaged masks will be used instead."
+        ),
+    )
     optoptions.add_argument(
         "-den",
         dest="den_type",

aroma/features.py

@@ -1,11 +1,9 @@
 """Functions to calculate ICA-AROMA features for component classification."""
-import os
-
 import nibabel as nib
 import numpy as np
 from nilearn import image, masking
 
-from .utils import cross_correlation, get_resource_path
+from .utils import cross_correlation
 
 
 def feature_time_series(mel_mix, mc):
@@ -151,7 +149,7 @@ def feature_frequency(mel_FT_mix, TR):
     return HFC
 
 
-def feature_spatial(mel_IC):
+def feature_spatial(mel_IC, masks):
     """Extract the spatial feature scores.
 
     For each IC it determines the fraction of the mixture modeled thresholded
@@ -163,6 +161,9 @@ def feature_spatial(mel_IC):
     mel_IC : str
         Full path of the nii.gz file containing mixture-modeled thresholded
         (p<0.5) Z-maps, registered to the MNI152 2mm template
+    masks : dict
+        Dictionary of masks. Keys are mask names and values are img_like
+        objects.
 
     Returns
     -------
@@ -177,11 +178,6 @@ def feature_spatial(mel_IC):
     mel_IC_img = nib.load(mel_IC)
     num_ICs = mel_IC_img.shape[3]
 
-    masks_dir = get_resource_path()
-    csf_mask = os.path.join(masks_dir, "mask_csf.nii.gz")
-    edge_mask = os.path.join(masks_dir, "mask_edge.nii.gz")
-    out_mask = os.path.join(masks_dir, "mask_out.nii.gz")
-
     # Loop over ICs
     edge_fract = np.zeros(num_ICs)
     csf_fract = np.zeros(num_ICs)
@@ -203,17 +199,17 @@ def feature_spatial(mel_IC):
 
         # Get sum of Z-values of the voxels located within the CSF
         # (calculate via the mean and number of non-zero voxels)
-        csf_data = masking.apply_mask(temp_IC, csf_mask)
+        csf_data = masking.apply_mask(temp_IC, masks["csf"])
         csf_sum = np.sum(csf_data)
 
         # Get sum of Z-values of the voxels located within the Edge
         # (calculate via the mean and number of non-zero voxels)
-        edge_data = masking.apply_mask(temp_IC, edge_mask)
+        edge_data = masking.apply_mask(temp_IC, masks["edge"])
         edge_sum = np.sum(edge_data)
 
         # Get sum of Z-values of the voxels located outside the brain
         # (calculate via the mean and number of non-zero voxels)
-        out_data = masking.apply_mask(temp_IC, out_mask)
+        out_data = masking.apply_mask(temp_IC, masks["out"])
         out_sum = np.sum(out_data)
 
         # Determine edge and CSF fraction

aroma/utils.py

@@ -2,10 +2,60 @@
 import os
 import os.path as op
 import shutil
+from tempfile import mkstemp
 
 import nibabel as nib
 import numpy as np
 from nilearn import image, masking
+from scipy import ndimage
+
+
+def derive_masks(in_file, csf=None):
+    """Estimate or load necessary masks based on inputs.
+
+    Parameters
+    ----------
+    in_file : str
+        4D EPI file.
+    csf : None or str, optional
+        CSF tissue probability map or binary mask.
+        If None, use precomputed, standard space masks packaged with AROMA.
+
+    Returns
+    -------
+    masks : dict
+        Dictionary with the different masks as img_like objects.
+    """
+    if csf is None:
+        print("No CSF TPM/mask provided. Using packaged masks.")
+        mask_dir = get_resource_path()
+        csf_img = nib.load(op.join(mask_dir, "mask_csf.nii.gz"))
+        out_img = nib.load(op.join(mask_dir, "mask_out.nii.gz"))
+        brain_img = image.math_img("1 - mask", mask=out_img)
+        edge_img = nib.load(op.join(mask_dir, "mask_edge.nii.gz"))
+    else:
+        brain_img = masking.compute_epi_mask(in_file)
+        out_img = image.math_img("1 - mask", mask=brain_img)
+        csf_img = nib.load(csf)
+        csf_data = csf_img.get_fdata()
+        if len(np.unique(csf_data)) == 2:
+            print("CSF mask provided. Inferring other masks.")
+        else:
+            print("CSF TPM provided. Inferring CSF and other masks.")
+            csf_img = image.math_img("csf >= 0.95", csf=csf_img)
+        gmwm_img = image.math_img("(brain - csf) > 0", brain=brain_img, csf=csf_img)
+        gmwm_data = gmwm_img.get_fdata()
+        eroded_data = ndimage.binary_erosion(gmwm_data, iterations=4)
+        edge_data = gmwm_data - eroded_data
+        edge_img = nib.Nifti1Image(edge_data, gmwm_img.affine, header=gmwm_img.header)
+
+    masks = {
+        "brain": brain_img,
+        "csf": csf_img,
+        "edge": edge_img,
+        "out": out_img,
+    }
+    return masks
 
 
 def runICA(fsl_dir, in_file, out_dir, mel_dir_in, mask, dim, TR):
@@ -30,17 +80,30 @@ def runICA(fsl_dir, in_file, out_dir, mel_dir_in, mask, dim, TR):
     TR : float
         TR (in seconds) of the fMRI data
 
+    Returns
+    -------
+    mel_IC_thr : str
+        Path to 4D nifti file containing thresholded component maps.
+    mel_IC_mix : str
+        Path to mixing matrix.
+    mel_FT_mix : str
+        Path to component power spectrum file.
+
     Output
     ------
     melodic.ica/: MELODIC directory
     melodic_IC_thr.nii.gz: Merged file containing the mixture modeling
                            thresholded Z-statistical maps located in
                            melodic.ica/stats/
     """
+    temp_file = mkstemp(suffix=".nii.gz")
+    mask.to_filename(temp_file)
+
     # Define the 'new' MELODIC directory and predefine some associated files
     mel_dir = op.join(out_dir, 'melodic.ica')
     mel_IC = op.join(mel_dir, 'melodic_IC.nii.gz')
     mel_IC_mix = op.join(mel_dir, 'melodic_mix')
+    mel_FT_mix = op.join(mel_dir, "melodic_FTmix")
     mel_IC_thr = op.join(out_dir, 'melodic_IC_thr.nii.gz')
 
     # When a MELODIC directory is specified,
@@ -98,7 +161,7 @@ def runICA(fsl_dir, in_file, out_dir, mel_dir_in, mask, dim, TR):
                                op.join(fsl_dir, 'melodic'),
                                in_file,
                                mel_dir,
-                               mask,
+                               temp_file,
                                dim,
                                TR
                            )
@@ -138,95 +201,8 @@ def runICA(fsl_dir, in_file, out_dir, mel_dir_in, mask, dim, TR):
     )
     zstat_4d_img.to_filename(mel_IC_thr)
 
-
-def register2MNI(fsl_dir, in_file, out_file, affmat, warp):
-    """Register an image (or time-series of images) to MNI152 T1 2mm.
-
-    If no affmat is defined, it only warps (i.e. it assumes that the data has
-    been registered to the structural scan associated with the warp-file
-    already). If no warp is defined either, it only resamples the data to 2mm
-    isotropic if needed (i.e. it assumes that the data has been registered to
-    a MNI152 template). In case only an affmat file is defined, it assumes that
-    the data has to be linearly registered to MNI152 (i.e. the user has a
-    reason not to use non-linear registration on the data).
-
-    Parameters
-    ----------
-    fsl_dir : str
-        Full path of the bin-directory of FSL
-    in_file : str
-        Full path to the data file (nii.gz) which has to be registerd to
-        MNI152 T1 2mm
-    out_file : str
-        Full path of the output file
-    affmat : str
-        Full path of the mat file describing the linear registration (if data
-        is still in native space)
-    warp : str
-        Full path of the warp file describing the non-linear registration (if
-        data has not been registered to MNI152 space yet)
-
-    Output
-    ------
-    melodic_IC_mm_MNI2mm.nii.gz : merged file containing the mixture modeling
-                                  thresholded Z-statistical maps registered to
-                                  MNI152 2mm
-    """
-    # Define the MNI152 T1 2mm template
-    fslnobin = fsl_dir.rsplit('/', 2)[0]
-    ref = op.join(fslnobin, 'data', 'standard', 'MNI152_T1_2mm_brain.nii.gz')
-
-    # If the no affmat- or warp-file has been specified, assume that the data
-    # is already in MNI152 space. In that case only check if resampling to
-    # 2mm is needed
-    if not affmat and not warp:
-        in_img = nib.load(in_file)
-        # Get 3D voxel size
-        pixdim1, pixdim2, pixdim3 = in_img.header.get_zooms()[:3]
-
-        # If voxel size is not 2mm isotropic, resample the data, otherwise
-        # copy the file
-        if (pixdim1 != 2) or (pixdim2 != 2) or (pixdim3 != 2):
-            os.system(' '.join([op.join(fsl_dir, 'flirt'),
-                                ' -ref ' + ref,
-                                ' -in ' + in_file,
-                                ' -out ' + out_file,
-                                ' -applyisoxfm 2 -interp trilinear']))
-        else:
-            os.copyfile(in_file, out_file)
-
-    # If only a warp-file has been specified, assume that the data has already
-    # been registered to the structural scan. In that case apply the warping
-    # without a affmat
-    elif not affmat and warp:
-        # Apply warp
-        os.system(' '.join([op.join(fsl_dir, 'applywarp'),
-                            '--ref=' + ref,
-                            '--in=' + in_file,
-                            '--out=' + out_file,
-                            '--warp=' + warp,
-                            '--interp=trilinear']))
-
-    # If only a affmat-file has been specified perform affine registration to
-    # MNI
-    elif affmat and not warp:
-        os.system(' '.join([op.join(fsl_dir, 'flirt'),
-                            '-ref ' + ref,
-                            '-in ' + in_file,
-                            '-out ' + out_file,
-                            '-applyxfm -init ' + affmat,
-                            '-interp trilinear']))
-
-    # If both a affmat- and warp-file have been defined, apply the warping
-    # accordingly
-    else:
-        os.system(' '.join([op.join(fsl_dir, 'applywarp'),
-                            '--ref=' + ref,
-                            '--in=' + in_file,
-                            '--out=' + out_file,
-                            '--warp=' + warp,
-                            '--premat=' + affmat,
-                            '--interp=trilinear']))
+    os.remove(temp_file)
+    return mel_IC_thr, mel_IC_mix, mel_FT_mix
 
 
 def cross_correlation(a, b):