[ENH] Derive masks to support native-space data

aroma/aroma.py

@@ -28,7 +28,8 @@ def aroma_workflow(
  overwrite=False,
  generate_plots=True,
  debug=False,
- quiet=False
+ quiet=False,
+ csf=None,
 ):
  """Run the AROMA workflow.
 
@@ -188,56 +189,28 @@ 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)
+ masks = utils.derive_masks(in_file, csf=csf)
 
  # 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, masks["brain"], 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
+ mel_IC_MNI, masks
  )
 
  LGR.info(" - extracting the Maximum RP correlation feature")
- mel_mix = op.join(out_dir, "melodic.ica", "melodic_mix")
- 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)
@@ -252,6 +225,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",
@@ -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,12 +1,11 @@
 """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 .utils import cross_correlation, get_resource_path
+from .utils import cross_correlation
 
 LGR = logging.getLogger(__name__)
 
@@ -154,7 +153,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
@@ -166,6 +165,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,11 +182,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)
@@ -206,17 +203,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