[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)