Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[ENH] Derive masks to support native-space data #24

Open
wants to merge 7 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from 4 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
54 changes: 16 additions & 38 deletions aroma/aroma.py
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -28,7 +28,9 @@ def aroma_workflow(
overwrite=False,
generate_plots=True,
debug=False,
quiet=False
quiet=False,
csf=None,
brain=None,
):
"""Run the AROMA workflow.

Expand Down Expand Up @@ -168,7 +170,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
Expand All @@ -188,56 +190,32 @@ 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
mel_IC_MNI, csf_img, brain_img, edge_img, out_img
)

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)
Expand All @@ -252,6 +230,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")
22 changes: 19 additions & 3 deletions aroma/cli/aroma.py
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand All @@ -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."
),
)

Expand Down Expand Up @@ -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",
Expand Down
19 changes: 8 additions & 11 deletions aroma/features.py
Original file line number Diff line number Diff line change
@@ -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 nilearn._utils import load_niimg

from .utils import cross_correlation, get_resource_path
from .utils import cross_correlation

LGR = logging.getLogger(__name__)

Expand Down Expand Up @@ -154,7 +153,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
Expand All @@ -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
-------
Expand All @@ -177,20 +179,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 = 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)
Expand Down
Loading