options to create csf_mask from T2star map, and edge and out-of-the-b…

…rain mask from brain mask given by the user

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,13 +214,15 @@ 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,
         motion_params,
         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)

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,20 +441,22 @@ 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
     ----------
     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,