Merge pull request #1 from chrisfilo/enh

Fixes and cleanup -- Thanks Chris!

Author: ccraddock

__init__.py qap/__init__.py__init__.py renamed to qap/__init__.py

@@ -56,7 +56,6 @@ def load_mask(mask_file, ref_file=None):
     """
     mask_img    = nib.load(mask_file)
     mask_dat    = mask_img.get_data()
-    ref_img     = nib.load(ref_file)
 
     # Check that the specified mask is binary.
     mask_vals   = np.unique(mask_dat)
@@ -65,6 +64,7 @@ def load_mask(mask_file, ref_file=None):
         raise Exception("")
 
     if ref_file is not None:
+        ref_img     = nib.load(ref_file)
         # Verify that the mask and anatomical images have the same dimensions.
         if ref_img.shape != mask_img.shape:
             raise Exception("Error: Mask and anatomical image are different dimensions")

dvars.py qap/dvars.pydvars.py renamed to qap/dvars.py

@@ -3,7 +3,8 @@
 """
 
 from spatial_qc import summary_mask, snr, cnr, fber, efc, artifacts, fwhm, ghost_all
-from temporal_qc import mean_dvars_wrapper, mean_fd_wrapper, mean_outlier_timepoints, mean_quality_timepoints
+from temporal_qc import mean_dvars_wrapper, mean_outlier_timepoints, mean_quality_timepoints, summarize_fd
+from load import load_image, load_mask
 
 
 ###
@@ -94,7 +95,7 @@ def run_qc_temporal(func_path, mask_path,
 
     # Mean FD (Jenkinson)
     if verbose: print "...mean FD"
-    mean_fd     = mean_fd_wrapper(motion_matrix)
+    mean_fd     = summarize_fd(motion_matrix)
 
     # 3dTout
     if verbose: print "...3dTout"

load.py qap/load.pyload.py renamed to qap/load.py

@@ -159,12 +159,12 @@ def artifacts(anat_data, fg_mask_data, calculate_qi2=False):
 
         # divide by the standard deviation
         bgNoiseZ    = bgNoise / bgNoise.std()
-        bgChiParams = ss.chi.fit(bgNoiseZ)
+        bgChiParams = stats.chi.fit(bgNoiseZ)
         #print bgChiParams
 
         # now generate values that are consistent with the histogram
         yx          = range(0,H.size)/bgNoise.std()
-        rvs         = ss.chi.pdf(yx,bgChiParams[0],loc=bgChiParams[1],scale=bgChiParams[2])
+        rvs         = stats.chi.pdf(yx,bgChiParams[0],loc=bgChiParams[1],scale=bgChiParams[2])
 
         # now we can calculate the goodness of fit
         gof         = np.average(np.absolute(H[halfMaxRightTail:]-rvs[halfMaxRightTail:]))
@@ -251,8 +251,8 @@ def ghost_mask(epi_data, mask_data, direction, ref_file=None, out_file=None):
 
     # Save mask
     if ref_file is not None and out_file is not None:
-        ref = nib.load(ref_file)
-        out = nib.Nifti1Image(n2_mask_data, ref.get_affine(), ref.get_header()) 
+        ref = nb.load(ref_file)
+        out = nb.Nifti1Image(n2_mask_data, ref.get_affine(), ref.get_header()) 
         out.to_filename(out_file)
 
     return(n2_mask_data)

qc.py qap/qc.pyqc.py renamed to qap/qc.py

@@ -13,6 +13,7 @@
 
 # DVARS
 from dvars import mean_dvars_wrapper
+from utils import remove_oblique_warning
 
 # MeanFD
 ## borrowed from C-PAC
@@ -31,83 +32,51 @@ def fd_jenkinson(in_file):
     Method to calculate Framewise Displacement (FD) calculations
     (Jenkinson et al., 2002)
     Parameters; in_file : string
-    Returns; out_file : string
+    Returns; FDs : numpy.array
     NOTE: infile should have one 3dvolreg affine matrix in one row - NOT the motion parameters
     '''
     import numpy as np
     import os
     import sys
     import math
 
-    out_file = os.path.join(os.getcwd(), 'FD_J.1D')
-        
-    f = open(out_file, 'w')
-    #print in_file
     pm_ = np.genfromtxt(in_file)
 
     pm = np.zeros((pm_.shape[0],pm_.shape[1]+4))
     pm[:,:12]=pm_
     pm[:,12:]=[0.0, 0.0, 0.0, 1.0]
-       
-    flag = 0
-
-    #The default radius (as in FSL) of a sphere represents the brain
-    rmax = 80.0
 
-    #rigid body transformation matrix
-    T_rb_prev = np.matrix(np.eye(4))
+    # The default radius (as in FSL) of a sphere represents the brain
+    rmax = 50.0
 
-    for i in range(0, pm.shape[0]):
+    FDs = []
+    T_rb_prev = np.matrix(pm[0].reshape(4,4))
+    for i in range(1, pm.shape[0]):
         T_rb = np.matrix(pm[i].reshape(4,4)) # making use of the fact that the order of aff12 matrix is "row-by-row"
-        
-        if flag == 0:
-            flag = 1
-            # first timepoint
-            print >> f, 0
-        else:
-            M = np.dot(T_rb, T_rb_prev.I) - np.eye(4)
-            A = M[0:3, 0:3]
-            b = M[0:3, 3]
+        M = np.dot(T_rb, T_rb_prev.I) - np.eye(4)
+        A = M[0:3, 0:3]
+        b = M[0:3, 3]
 
-            FD_J = math.sqrt((rmax*rmax/5)*np.trace(np.dot(A.T, A)) + np.dot(b.T, b))
-            print >> f, '%.8f'%FD_J
-                
+        FD_J = math.sqrt((rmax*rmax/5)*np.trace(np.dot(A.T, A)) + np.dot(b.T, b))
+        FDs.append(FD_J)         
         T_rb_prev = T_rb
-    
-    f.close()
-    
-    return out_file
+        
+    return np.array(FDs)
 
-def summarize_fd(in_file, fd_out_file=None, threshold=0.2):
+def summarize_fd(in_file, threshold=0.2):
     # Threshold is in terms of mm, i think?
 
-    # Run in temporary working directory
-    tmpdir = mkdtemp()
-    curdir = os.getcwd()
-    os.chdir(tmpdir)
-    
     # Compute FD
-    out_file    = fd_jenkinson(in_file)
-    fd          = np.loadtxt(out_file)
+    fd    = fd_jenkinson(in_file)
 
     # Calculate Mean
     mean_fd     = fd.mean()
 
     # Calculate Outliers
     ## Number and Percent of frames (time points) where 
     ## movement (FD) exceeded threshold
-    num_fd      = np.float((fd>threshold).sum())
-    percent_fd  = (num_fd*100)/(len(fd)+1)    
-    
-    #data = np.loadtxt("/data/Projects/ABIDE_Initiative/CPAC/Output_2014-04-29_preproc/pipeline_abide_func_preproc/0050002_session_1/frame_wise_displacement/_scan_rest_1_rest/FD.1D")    
-    
-    # Clean-Up
-    if fd_out_file:
-        os.rename(out_file, fd_out_file)
-    else:
-        os.remove(out_file)
-    os.chdir(curdir)
-    os.rmdir(tmpdir)
+    num_fd      = (fd>threshold).sum()
+    percent_fd  = (num_fd*100.0)/(len(fd))    
 
     return (mean_fd, num_fd, percent_fd)
 
@@ -147,11 +116,13 @@ def outlier_timepoints(func_file, mask_file, out_fraction=True):
     # or -legendre to remove any trend
     cmd     = "3dToutcount %s" % str_opts
     out     = commands.getoutput(cmd)
-    
+    out     = remove_oblique_warning(out)
+
     # Extract time-series in output
     lines   = out.splitlines()
     ## remove general information
     lines   = [ l for l in lines if l[:2] != "++" ]
+    
     ## string => floats
     outliers= [ float(l.strip()) for l in lines ] # note: don't really need strip
 
@@ -164,32 +135,35 @@ def mean_outlier_timepoints(*args, **kwrds):
 
 
 # 3dTqual
-def quality_timepoints(func_file, automask=True):
+def quality_timepoints(func_file, mask=None):
     """
     Calculates a 'quality index' for each timepoint in the 4D functional dataset.
     Low values are good and indicate that the timepoint is not very different from the norm.
     """
     import subprocess
 
     opts    = []
-    if automask:
-        opts.append("-automask")
+    if mask:
+        if mask=="auto":
+            opts.append("-automask")
+        else:
+            opts.append("-mask %s"%mask)
     opts.append(func_file)
     str_opts= " ".join(opts)
 
     cmd     = "3dTqual %s" % str_opts
     p       = subprocess.Popen(cmd.split(" "), 
                                 stdout=subprocess.PIPE, 
                                 stderr=subprocess.PIPE)
-    out,err = p.communicate()
+    out, _ = p.communicate()
 
     #import code
     #code.interact(local=locals())
 
     # Extract time-series in output
     lines   = out.splitlines()
     ## remove general information
-    lines   = [ l for l in lines if l[:2] != "++" ]
+    lines   = [ l for l in lines if l[:2] not in "++" ]
     ## string => floats
     outliers= [ float(l.strip()) for l in lines ] # note: don't really need strip
 
@@ -199,3 +173,14 @@ def mean_quality_timepoints(*args, **kwrds):
     qualities       = quality_timepoints(*args, **kwrds)
     mean_qualities  = np.mean(qualities)
     return mean_qualities
+
+def median_tsnr(func_data):
+    """ Calculates median of temporal Signal to Noise ratio within a mask"""
+    
+    data_std = func_data.std(axis=0)
+    # exclude voxels with no variance
+    std_mask = (data_std !=0)
+    tsnr = func_data.mean(axis=0)[std_mask]/data_std[std_mask]
+    print tsnr.shape
+    
+    return np.median(tsnr)

spatial_qc.py qap/spatial_qc.pyspatial_qc.py renamed to qap/spatial_qc.py

@@ -86,3 +86,11 @@ def gen_file_map(sink_dir, resource_paths=None):
 
     return analysis_map
 
+def remove_oblique_warning(out_str):
+    new_lines = []
+    for i, line in enumerate(out_str.splitlines()):
+        if line.startswith("*+ WARNING:   If you are performing spatial transformations on an oblique dset,"):
+            break
+        new_lines.append(line)
+    new_lines += out_str.splitlines()[i+7:]
+    return "\n".join(new_lines)

temporal_qc.py qap/temporal_qc.pytemporal_qc.py renamed to qap/temporal_qc.py

@@ -0,0 +1,2 @@
+nibabel
+nitime

utils.py qap/utils.pyutils.py renamed to qap/utils.py

@@ -118,7 +118,7 @@ The spatial functional measures will make use of the mean functional image and i
 You should know the phase encoding direction to decide if you want to use `func_ghost_x` (RL/LR) or `func_ghost_y` (AP/PA).
 
 	from qap import ghost_all
-	func_ghost_x,func_ghost_y = ghost_all(mean_func_data, func_mask_data)
+	func_ghost_x, func_ghost_y, _ = ghost_all(mean_func_data, func_mask_data)
 
 
 ## Temporal Functional
@@ -129,6 +129,7 @@ You should know the phase encoding direction to decide if you want to use `func_
 * **Mean Fractional Displacement - Jenkinson [func_mean_fd]:** A measure of subject head motion, which compares the motion between the current and previous volumes. This is calculated by summing the absolute value of displacement changes in the x, y and z directions and rotational changes about those three axes. The rotational changes are given distance values based on the changes across the surface of a 50mm radius sphere [^5][^8]. _Uses functional time-series._
 * **Number of volumes with FD greater than 0.2mm [func_num_fd]:** _Uses functional time-series._
 * **Percent of volumes with FD greater than 0.2mm [func_perc_fd]:** _Uses functional time-series._
+* **Median temporal Signal to Noise ratio [tsnr]:** _Uses functional time-series and a mask._
 
 
 ### Standardized DVARS
@@ -147,17 +148,24 @@ The `out_fraction` option if `True` will return the mean _fraction_ of time-poin
 
 ### Median Distance Index
 
-The `auto_mask` option if `True` will automatically compute the brain mask from the data otherwise the whole dataset will be used.
+The `mask` allows you to specify a mask file. If set to "auto" mask will be calculated from data. If not set the whole dataset will be used.
 
 	from qap import mean_quality_timepoints
-	func_quality	= mean_quality_timepoints(func_file, automask=True)
+	func_quality	= mean_quality_timepoints(func_file, mask="auto")
 
 ### FD
 
 Here we describe computing `mean_fd`, `num_fd`, and `perc_fd`. This requires that you have the input coordinate transforms that is output by AFNI's `3dvolreg` during motion correction. The option `threshold` sets the threshold for determining the number and percent of volumes with FD greater than said threshold.
 
-	from qap import mean_fd_wrapper
+	from qap import summarize_fd
 	mean_fd, num_fd, perc_fd = summarize_fd(motion_matrix_file, threshold=0.2)
+	
+### tSNR
+
+tSNR is a ratio of signal (mean across time) and noise (standard deviation around the mean across time) averaged over all voxels withing the specified mask
+
+    from qap import median_tsnr
+	tsnr = median_tsnr(func_data, func_mask_data)
 
 
 ## Determining Outliers