[ENH] Derive masks to support native-space data

aroma/aroma.py

@@ -5,8 +5,8 @@
 import os.path as op
 import shutil
 
-import nibabel as nib
 import pandas as pd
+from nilearn._utils import load_niimg
 
 from aroma import utils, features, _version
 
@@ -29,6 +29,8 @@ def aroma_workflow(
     generate_plots=True,
     debug=False,
     quiet=False,
+    csf=None,
+    brain=None,
     mc_source="auto",
 ):
     """Run the AROMA workflow.
@@ -169,7 +171,7 @@ def aroma_workflow(
     fsl_dir = op.join(os.environ["FSLDIR"], "bin", "")
     # Get TR of the fMRI data, if not specified
     if not TR:
-        in_img = nib.load(in_file)
+        in_img = load_niimg(in_file)
         TR = in_img.header.get_zooms()[3]
 
     # Check TR
@@ -189,57 +191,37 @@ 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:
-            LGR.warning(
-                " - 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)
+    brain_img, csf_img, out_img, edge_img = utils.load_masks(
+        in_file, csf=csf, brain=brain
+    )
 
     # Run ICA-AROMA
     LGR.info("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, brain_img, dim, TR
+    )
 
     LGR.info("Step 2) Automatic classification of the components")
     LGR.info("  - 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)
 
     LGR.info("  - extracting the CSF & Edge fraction features")
     features_df = pd.DataFrame()
     features_df["edge_fract"], features_df["csf_fract"] = features.feature_spatial(
-        mel_IC_MNI
+        z_maps=mel_IC_MNI,
+        csf_mask=csf_img,
+        brain_mask=brain_img,
+        edge_mask=edge_img,
+        out_mask=out_img,
     )
 
     LGR.info("  - extracting the Maximum RP correlation feature")
-    mel_mix = op.join(out_dir, "melodic.ica", "melodic_mix")
     mc = utils.load_motpars(mc, source=mc_source)
-    features_df["max_RP_corr"] = features.feature_time_series(mel_mix, mc)
+    features_df["max_RP_corr"] = features.feature_time_series(mixing, mc)
 
     LGR.info("  - extracting the High-frequency content feature")
-    mel_FT_mix = op.join(out_dir, "melodic.ica", "melodic_FTmix")
-    features_df["HFC"] = features.feature_frequency(mel_FT_mix, TR)
+    features_df["HFC"] = features.feature_frequency(mixing_FT, TR)
 
     LGR.info("  - classification")
     motion_ICs = utils.classification(features_df, out_dir)
@@ -254,6 +236,6 @@ def aroma_workflow(
     if den_type != "no":
         LGR.info("Step 3) Data denoising")
         utils.denoising(fsl_dir, in_file, out_dir,
-                        mel_mix, den_type, motion_ICs)
+                        mixing, den_type, motion_ICs)
 
     LGR.info("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",
@@ -48,8 +52,9 @@ def _get_parser():
         default=None,
         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."
         ),
     )
 
@@ -130,6 +135,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,10 +1,9 @@
 """Functions to calculate ICA-AROMA features for component classification."""
 import logging
-import os
 
-import nibabel as nib
 import numpy as np
 from nilearn import image, masking
+from nilearn._utils import load_niimg
 
 from . import utils
 
@@ -157,7 +156,7 @@ def feature_frequency(mel_FT_mix, TR):
     return HFC
 
 
-def feature_spatial(mel_IC):
+def feature_spatial(*, z_maps, csf_mask, edge_mask, brain_mask, out_mask):
     """Extract the spatial feature scores.
 
     For each IC it determines the fraction of the mixture modeled thresholded
@@ -169,6 +168,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
     -------
@@ -180,20 +182,15 @@ def feature_spatial(mel_IC):
         mel_IC file
     """
     # Get the number of ICs
-    mel_IC_img = nib.load(mel_IC)
+    mel_IC_img = load_niimg(z_maps)
     num_ICs = mel_IC_img.shape[3]
 
-    masks_dir = utils.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)
     for i in range(num_ICs):
         # Extract IC from the merged melodic_IC_thr2MNI2mm file
-        temp_IC = image.index_img(mel_IC, i)
+        temp_IC = image.index_img(mel_IC_img, i)
 
         # Change to absolute Z-values
         temp_IC = image.math_img("np.abs(img)", img=temp_IC)

aroma/tests/test_features.py

@@ -1,8 +1,7 @@
+"""Tests for the aroma.features module."""
 import numpy as np
-import os.path as op
-import pandas as pd
 
-from aroma import features
+from aroma import features, utils
 
 
 def test_feature_time_series(mel_mix, mc, max_correls):
@@ -35,8 +34,19 @@ def test_feature_spatial(mel_IC, edgeFract, csfFract):
 
     np.random.seed(1)
 
+    # Get masks
+    brain_img, csf_img, out_img, edge_img = utils.load_masks(
+        mel_IC, csf=None, brain=None
+    )
+
     # Run feature_spatial
-    edge_fract, csf_fract = features.feature_spatial(mel_IC)
+    edge_fract, csf_fract = features.feature_spatial(
+        z_maps=mel_IC,
+        csf_mask=csf_img,
+        brain_mask=brain_img,
+        edge_mask=edge_img,
+        out_mask=out_img,
+    )
 
     # Read features csv
     edgeFract = np.load(edgeFract)