Added sliding window and dvars(mask) approaches to outfind.py for finding outliers, modified dvars() and added a new function dvars_voxel() in metrics.py and modified mad_time_detector() in detectors.py

Author: jagruti8

findoutlie/detectors.py

@@ -45,7 +45,7 @@ def mad_voxel_detector(img,threshold=3.5):
     outlier_tf = np.abs(img-med)>(threshold*mad)
     return outlier_tf
 
-def mad_time_detector(measures, threshold=3.5):
+def mad_time_detector(measures, lower_bound, threshold=3.5):
     """ Detect outliers in 'measures' using median absolute deviation.
     Returns 1D vector of same length as 'measures', where True means the corresponsding 
     value in 'measures' is an outlier.
@@ -71,7 +71,10 @@ def mad_time_detector(measures, threshold=3.5):
     # Calculate median absoulte deviation of measures
     mad = np.median(np.abs(measures-med))
     # Calculate the outliers
-    outlier_tf = measures>med+threshold*mad
+    if lower_bound:
+        outlier_tf = np.abs(measures-med)>threshold*mad
+    else:
+        outlier_tf = measures>med+threshold*mad
     return outlier_tf
 
 def iqr_detector(measures, iqr_proportion=1.5):

findoutlie/metrics.py

@@ -5,6 +5,26 @@
 # +++your code here+++
 import numpy as np
 
+def dvars_voxel(voxels):
+    """ Calculate dvars metric on 2D array with voxels in rows and time-points(volumes) in columns
+
+    The dvars calculation between two volumes is defined as the square root of 
+    (the mean of the (voxel differences square)).
+
+    Parameters
+    ----------
+    voxels : 2D array
+
+    Returns
+    -------
+    dvals : 1D array
+        One-dimensional array with n-1 elements, where n is the number of 
+        volumes in 'img'.
+    """
+    vol_diff = voxels[..., 1:] - voxels[..., :-1]
+    dvar_val = np.sqrt(np.mean(vol_diff ** 2, axis=0))
+    return dvar_val
+
 
 def dvars(img):
     """ Calculate dvars metric on Nibabel image `img`
@@ -33,14 +53,8 @@ def dvars(img):
     data = img.get_fdata()
 
     voxel_by_time = np.reshape(data, (-1, data.shape[-1]))
+    dvar_val = dvars_voxel(voxel_by_time)
 
-    vol_diff = voxel_by_time[..., 1:] - voxel_by_time[..., :-1]  # 2D array
-    # print(vol_diff.shape())
-    # print(vol_diff)
-    # vol_diff_1D=vol_diff.flatten()
-    dvar_val = np.sqrt(np.mean(vol_diff ** 2, axis=0))
-    # print(dvar_val.shape())
-    # print(dvar_val)
     return dvar_val
 
     raise NotImplementedError("Code up this function")

findoutlie/outfind.py

@@ -9,7 +9,7 @@
 
 from skimage.filters import threshold_otsu
 
-from .metrics import dvars
+from .metrics import dvars,dvars_voxel
 from .detectors import iqr_detector,mad_voxel_detector,mad_time_detector
 
 def segment_brain(img):
@@ -27,14 +27,12 @@ def segment_brain(img):
     mean_img = np.mean(img, axis=-1)
     # calculate the threshold for segmenting brain from background
     threshold = threshold_otsu(mean_img)
-    mask = np.expand_dims(mean_img > threshold, axis=1)
-    mask_2D = np.tile(mask, (1, img.shape[-1]))
-    thresholded_img = np.where(mask_2D, img, np.nan)
+    mask = mean_img > threshold
     # filter only brain voxels
-    brain_voxels = thresholded_img[~np.isnan(thresholded_img).all(axis=1)]
+    brain_voxels = img[mask]
     return brain_voxels
 
-def detect_outliers_mean_absolute_deviation_mask(fname):
+def detect_outliers_mad_median_absolute_deviation_mask(fname):
     """ Detect outliers given image file path 'filename'
      
     Parameters
@@ -47,7 +45,7 @@ def detect_outliers_mean_absolute_deviation_mask(fname):
     outliers : array
         Indices of outlier volumes.
     """
-    # A mask is used to first segment the brain regions from the background, then mean absolute deviation is used to detect outliers
+    # A mask is used to first segment the brain regions from the background, then median absolute deviation is used to detect outliers
     img = nib.load(fname)
     img_data = img.get_fdata()
     # reshape from 4D to 2D
@@ -59,9 +57,71 @@ def detect_outliers_mean_absolute_deviation_mask(fname):
     # calculate the number of outlying voxels for each time point
     voxel_outliers_per_time = np.nansum(outliers_voxel,axis=0)
     # find the outliers in the time-series    
-    outliers_time = mad_time_detector(voxel_outliers_per_time)
+    outliers_time = mad_time_detector(voxel_outliers_per_time, lower_bound=False)
     # Return indices of True values from Boolean array. 
-    return np.nonzero(outliers_time)[0]
+    return np.nonzero(outliers_time)
+
+def detect_outliers_mad_dvars_mask(fname):
+    """ Detect outliers given image file path 'filename'
+    
+    Parameters
+    ----------
+    fname : str or Path
+        Filename of 4D image, as string or Path object
+    
+    Returns
+    -------
+    outliers : array
+        Indices of outlier volumes.
+    """
+    # A mask is used to first segment the brain regions from the background, dvars is calculated and then median absolute deviation is used to detect outliers 
+    img = nib.load(fname)
+    img_data = img.get_fdata()
+    # reshape from 4D to 2D
+    img_data_2D = np.reshape(img_data, (-1,img_data.shape[-1]))
+    # segment brain from background
+    brain_voxels = segment_brain(img_data_2D)
+    # calculate dvars
+    dvs = dvars_voxel(brain_voxels)
+    # detect outliers
+    is_outlier = mad_time_detector(dvs, lower_bound=True)
+    return np.nonzero(is_outlier)
+
+def detect_outliers_mad_sliding_window_mask(fname):
+    """ Detect outliers given image file path 'filename'
+    
+    Parameters
+    ----------
+    fname : str or Path
+        Filename of 4D image, as string or Path object
+    
+    Returns 
+    -------
+    outliers : array
+        Indices of outlier volumes.
+    """
+    # A mask is used to first segment the brain regions from the background, sliding window approach is used to detect outliers in each window using median absolute deviation
+    img = nib.load(fname)
+    img_data = img.get_fdata()
+    # reshape from 4D to 2D
+    img_data_2D = np.reshape(img_data, (-1,img_data.shape[-1]))
+    # segment brain from background
+    brain_voxels = segment_brain(img_data_2D)
+    # apply sliding window
+    overlap = 10
+    window_length = 20
+    outliers = []
+    for i in range(0, brain_voxels.shape[-1],overlap):
+       if i+window_length>=brain_voxels.shape[-1]:
+           elements = np.mean(brain_voxels[:,i:],axis=0)
+           outliers_1 = np.nonzero(mad_time_detector(elements, lower_bound=True))[0] + i
+           outliers.extend(outliers_1)
+           break
+       else:
+           elements = np.mean(brain_voxels[:,i:i+window_length],axis=0)
+           outliers_1 = np.nonzero(mad_time_detector(elements, lower_bound=True))[0] + i
+           outliers.extend(outliers_1)
+    return np.unique(outliers)
 
 
 def detect_outliers(fname):
@@ -82,7 +142,7 @@ def detect_outliers(fname):
     dvs = dvars(img)
     is_outlier = iqr_detector(dvs, iqr_proportion=2)
     # Return indices of True values from Boolean array.
-    return np.nonzero(is_outlier)[0]
+    return np.nonzero(is_outlier)
 
 
 def find_outliers(data_directory):
@@ -102,7 +162,13 @@ def find_outliers(data_directory):
     image_fnames = Path(data_directory).glob("**/sub-*.nii.gz")
     outlier_dict = {}
     for fname in image_fnames:
-        outliers = detect_outliers_mean_absolute_deviation_mask(fname)
+        # detect outliers using mad over voxels and mad over time points
+        outliers_mad = detect_outliers_mad_median_absolute_deviation_mask(fname)
+        # detect outliers using dvars and mad over time points
+        outliers_dvars = detect_outliers_mad_dvars_mask(fname)
+        # detect outliers using sliding window and mad over time points
+        outliers_sliding_window = np.array(detect_outliers_mad_sliding_window_mask(fname))
+        outliers = np.intersect1d(outliers_sliding_window,np.union1d(outliers_mad, outliers_dvars))
         #outliers = detect_outliers(fname)
         outlier_dict[fname] = outliers
     return outlier_dict