From b6ffaf05ba1175ee594da78f10299939a5bde00b Mon Sep 17 00:00:00 2001 From: Cesar Caballero Gaudes Date: Fri, 25 Nov 2022 16:28:08 +0100 Subject: [PATCH] options to create csf_mask from T2star map, and edge and out-of-the-brain mask from brain mask given by the user --- aroma/aroma.py | 42 +++++++++++++++++++++++++++++++++++++++++- aroma/utils.py | 24 +++++++++++++----------- 2 files changed, 54 insertions(+), 12 deletions(-) diff --git a/aroma/aroma.py b/aroma/aroma.py index 5c36d6e..85c2008 100644 --- a/aroma/aroma.py +++ b/aroma/aroma.py @@ -27,7 +27,13 @@ def aroma_workflow( debug=False, quiet=False, mc_source="auto", - f_hp=0.01 + f_hp=0.01, + T2star_map=None, + brain_mask=None, + T2star_thresh=0.08, + erode_csf_voxels=2, + erode_edge_voxels=2, + dilate_out_voxels=0 ): """Run the AROMA workflow. @@ -56,6 +62,21 @@ def aroma_workflow( What format is the mc file in? f_hp : float, optional High-pass cutoff frequency in spectrum computations. + T2star_map : str, optional + path to nifti file with T2star map. + brain_mask: str, optional + path to nifti file with brain mask + T2star_thresh: float, optional + mask will only consider voxels with T2-star values + above T2star_thresh. Default value of 0.080 sec assumes 3T + acquisition. + erode_voxels: : int, optional + number of voxels to erode (radius of ball) + erode_edge_voxels: int, optional + Number of voxels to erode in order to compute the edge mask. + dilate_out_voxels: int, optional + Number of voxels to erode in order to compute the out-of-brain mask + """ if not op.isfile(in_file): raise FileNotFoundError(f"Input file does not exist: {in_file}") @@ -167,6 +188,23 @@ def aroma_workflow( f"number of rows in motion parameters ({motion_params.shape[0]})." ) + # Compute spatial metrics (CSF & Edge fraction) + + # create csf mask + if (T2star_map is None) and (brain_mask is None) : + masks_dir = utils.get_resource_path() + csf_mask = os.path.join(masks_dir, "mask_csf.nii.gz") + else: + csf_mask = utils.csf_mask_from_T2star(T2star_map,brain_mask,out_dir,T2star_thresh,erode_csf_voxels) + + # create edge and out-of-the-brain masks + if brain_mask is None: + masks_dir = utils.get_resource_path() + edge_mask = os.path.join(masks_dir, "mask_edge.nii.gz") + out_mask = os.path.join(masks_dir, "mask_out.nii.gz") + else: + edge_mask,out_mask = utils.edge_out_mask(brain_mask,out_dir,erode_edge_voxels,dilate_out_voxels) + LGR.info(" - extracting the CSF & Edge fraction features") metric_metadata = {} features_df = pd.DataFrame() @@ -176,6 +214,7 @@ def aroma_workflow( metric_metadata ) = features.feature_spatial(component_maps, csf_mask, edge_mask, out_mask, metric_metadata) + # Compute temporal metrics (correlation with RP) LGR.info(" - extracting the Maximum RP correlation feature") features_df["max_RP_corr"], metric_metadata = features.feature_time_series( mixing, @@ -183,6 +222,7 @@ def aroma_workflow( metric_metadata, ) + # Compute frequency metrics (half-power frequency) LGR.info(" - extracting the High-frequency content feature") # Should probably check that the frequencies match up with MELODIC's outputs mel_FT_mix, FT_freqs = utils.get_spectrum(mixing, TR) diff --git a/aroma/utils.py b/aroma/utils.py index 83e11f5..fc4b5c2 100644 --- a/aroma/utils.py +++ b/aroma/utils.py @@ -3,7 +3,6 @@ import logging import os.path as op import shutil -import tempfile import nibabel as nib import numpy as np @@ -442,7 +441,7 @@ def get_spectrum(data: np.array, tr: float): idx = np.argsort(freqs) return power_spectrum[idx, :], freqs[idx] -def csf_mask_from_T2star(T2star_map:str, brain_mask:str, T2star_thresh:float= 0.080, erode_voxels:int=2): +def csf_mask_from_T2star(T2star_map:str, brain_mask:str, out_dir:str, T2star_thresh:float= 0.080, erode_voxels:int=2): """"Create mask of ventricle CSF based on T2-star map. Parameters @@ -450,12 +449,14 @@ def csf_mask_from_T2star(T2star_map:str, brain_mask:str, T2star_thresh:float= 0. T2star_map: str path to nifti file with T2star map. brain_mask: str - path to nifti file with brain mask + path to nifti file with brain mask. + out_dir: str + Output directory. T2star_thresh: :obj:`float` mask will only consider voxels with T2-star values above T2star_thresh. Default value of 0.080 sec assumes 3T acquisition. - erode_voxels: :obj:`int`` + erode_voxels: :obj:`int` number of voxels to erode (radius of ball) Returns @@ -471,7 +472,6 @@ def csf_mask_from_T2star(T2star_map:str, brain_mask:str, T2star_thresh:float= 0. # TODO: pass number of iterations as input argument # save to nifti file - temp_dir = tempfile.mkdtemp() path_csf_mask = op.join(temp_dir,"mask_csf.nii.gz") # Apply threshold to identify CSF voxels @@ -482,7 +482,7 @@ def csf_mask_from_T2star(T2star_map:str, brain_mask:str, T2star_thresh:float= 0. return path_csf_mask -def edge_out_mask(brain_mask:str, erode_edge_voxels:int=2, dilate_out_voxels:int=0): +def edge_out_mask(brain_mask:str, out_dir:str, erode_edge_voxels:int=2, dilate_out_voxels:int=0): """"Create mask of edge voxels and out-of-the-brain voxels. The edge mask will be the subtraction of the brain mask minus this mask eroded by erode_voxels, whereas the out_mask will @@ -492,10 +492,13 @@ def edge_out_mask(brain_mask:str, erode_edge_voxels:int=2, dilate_out_voxels:int Parameters ---------- - brain_img: path to nifti file with a brain to compute the within-brain mask + brain_mask: str, + path to nifti file with a brain to compute the within-brain mask. + out_dir: str, + Output directory. erode_edge_voxels: :obj:`int` Number of voxels to erode in order to compute the edge mask. - dilate_out_voxels: :obj:`int`` + dilate_out_voxels: :obj:`int` Number of voxels to erode in order to compute the out-of-brain mask Returns @@ -516,9 +519,8 @@ def edge_out_mask(brain_mask:str, erode_edge_voxels:int=2, dilate_out_voxels:int out_mask_bool = binary_dilation(brain_mask_bool, structure=ball(dilate_out_voxels)) # save to nifti file - temp_dir = tempfile.mkdtemp() - path_edge_mask = op.join(temp_dir,"mask_edge.nii.gz") - path_out_mask = op.join(temp_dir,"mask_out.nii.gz") + path_edge_mask = op.join(out_dir,"mask_edge.nii.gz") + path_out_mask = op.join(out_dir,"mask_out.nii.gz") # Compute masks edge_mask = nib.Nifti1Image(brain_edge_bool,