Source code for CPAC.easy_thresh.easy_thresh

import os
import re
import subprocess

import nibabel as nib
from CPAC.pipeline import nipype_pipeline_engine as pe
import nipype.interfaces.fsl as fsl
import nipype.interfaces.utility as util

[docs]def easy_thresh(wf_name): """ Workflow for carrying out cluster-based thresholding and colour activation overlaying Parameters ---------- wf_name : string Workflow name Returns ------- easy_thresh : object Easy thresh workflow object Notes ----- `Source <>`_ Workflow Inputs:: inputspec.z_stats : string (nifti file) z_score stats output for t or f contrast from flameo inputspec.merge_mask : string (nifti file) mask generated from 4D Merged derivative file inputspec.z_threshold : float Z Statistic threshold value for cluster thresholding. It is used to determine what level of activation would be statistically significant. Increasing this will result in higher estimates of required effect. inputspec.p_threshold : float Probability threshold for cluster thresholding. inputspec.paramerters : string (tuple) tuple containing which MNI and FSLDIR path information Workflow Outputs:: outputspec.cluster_threshold : string (nifti files) the thresholded Z statistic image for each t contrast outputspec.cluster_index : string (nifti files) image of clusters for each t contrast; the values in the clusters are the index numbers as used in the cluster list. outputspec.overlay_threshold : string (nifti files) 3D color rendered stats overlay image for t contrast After reloading this image, use the Statistics Color Rendering GUI to reload the color look-up-table outputspec.overlay_rendered_image : string (nifti files) 2D color rendered stats overlay picture for each t contrast outputspec.cluster_localmax_txt : string (text files) local maxima text file, defines the coordinates of maximum value in the cluster Order of commands: - Estimate smoothness of the image:: smoothest --mask= merge_mask.nii.gz --zstat=.../flameo/stats/zstat1.nii.gz arguments --mask : brain mask volume --zstat : filename of zstat/zfstat image - Create mask. For details see `fslmaths <>`_:: fslmaths ../flameo/stats/zstat1.nii.gz -mas merge_mask.nii.gz zstat1_mask.nii.gz arguments -mas : use (following image>0) to mask current image - Copy Geometry image dimensions, voxel dimensions, voxel dimensions units string, image orientation/origin or qform/sform info) from one image to another:: fslcpgeom MNI152_T1_2mm_brain.nii.gz zstat1_mask.nii.gz - Cluster based thresholding. For details see `FEAT <>`_:: cluster --dlh = 0.0023683100 --in = zstat1_mask.nii.gz --oindex = zstat1_cluster_index.nii.gz --olmax = zstat1_cluster_localmax.txt --othresh = zstat1_cluster_threshold.nii.gz --pthresh = 0.0500000000 --thresh = 2.3000000000 --volume = 197071 arguments --in : filename of input volume --dlh : smoothness estimate = sqrt(det(Lambda)) --oindex : filename for output of cluster index --othresh : filename for output of thresholded image --olmax : filename for output of local maxima text file --volume : number of voxels in the mask --pthresh : p-threshold for clusters --thresh : threshold for input volume Z statistic image is thresholded to show which voxels or clusters of voxels are activated at a particular significance level. A Z statistic threshold is used to define contiguous clusters. Then each cluster's estimated significance level (from GRF-theory) is compared with the cluster probability threshold. Significant clusters are then used to mask the original Z statistic image. - Get the maximum intensity value of the output thresholded image. This used is while rendering the Z statistic image:: fslstats zstat1_cluster_threshold.nii.gz -R arguments -R : output <min intensity> <max intensity> - Rendering. For details see `FEAT <>`_:: overlay 1 0 MNI152_T1_2mm_brain.nii.gz -a zstat1_cluster_threshold.nii.gz 2.30 15.67 zstat1_cluster_threshold_overlay.nii.gz slicer zstat1_cluster_threshold_overlay.nii.gz -L -A 750 zstat1_cluster_threshold_overlay.png The Z statistic range selected for rendering is automatically calculated by default, to run from red (minimum Z statistic after thresholding) to yellow (maximum Z statistic, here maximum intensity). High Level Workflow Graph: .. exec:: from CPAC.easy_thresh import easy_thresh wf = easy_thresh('easy_thresh_wf') wf.write_graph( graph2use='orig', dotfilename='./images/generated/' ) .. image:: ../../images/generated/easy_thresh.png :width: 800 Detailed Workflow Graph: .. image:: ../../images/generated/generated.easy_thresh_detailed.png :width: 800 Examples -------- >>> from CPAC.easy_thresh import easy_thresh >>> preproc = easy_thresh("new_workflow") >>> preproc.inputs.inputspec.z_stats= 'flameo/stats/zstat1.nii.gz' >>> preproc.inputs.inputspec.merge_mask = 'merge_mask/alff_Z_fn2standard_merged_mask.nii.gz' >>> preproc.inputs.inputspec.z_threshold = 2.3 >>> preproc.inputs.inputspec.p_threshold = 0.05 >>> preproc.inputs.inputspec.parameters = ('/usr/local/fsl/', 'MNI152') >>> # doctest: +SKIP """ easy_thresh = pe.Workflow(name=wf_name) inputnode = pe.Node(util.IdentityInterface(fields=['z_stats', 'merge_mask', 'z_threshold', 'p_threshold', 'parameters']), name='inputspec') outputnode = pe.Node(util.IdentityInterface(fields=['cluster_threshold', 'cluster_index', 'cluster_localmax_txt', 'overlay_threshold', 'rendered_image']), name='outputspec') # fsl easythresh # estimate image smoothness smooth_estimate = pe.MapNode(interface=fsl.SmoothEstimate(), name='smooth_estimate', iterfield=['zstat_file']) # run clustering after fixing stats header for talspace zstat_mask = pe.MapNode(interface=fsl.MultiImageMaths(), name='zstat_mask', iterfield=['in_file']) # operations to perform #-mas use (following image>0) to mask current image zstat_mask.inputs.op_string = '-mas %s' # fslcpgeom #copy certain parts of the header information (image dimensions, #voxel dimensions, voxel dimensions units string, image orientation/origin #or qform/sform info) from one image to another geo_imports = ['import subprocess'] copy_geometry = pe.MapNode(util.Function(input_names=['infile_a', 'infile_b'], output_names=['out_file'], function=copy_geom, imports=geo_imports), name='copy_geometry', iterfield=['infile_a', 'infile_b']) ##cluster-based thresholding #After carrying out the initial statistical test, the resulting #Z statistic image is then normally thresholded to show which voxels or #clusters of voxels are activated at a particular significance level. #A Z statistic threshold is used to define contiguous clusters. #Then each cluster's estimated significance level (from GRF-theory) is #compared with the cluster probability threshold. Significant clusters #are then used to mask the original Z statistic image for later production #of colour blobs.This method of thresholding is an alternative to #Voxel-based correction, and is normally more sensitive to activation. # cluster = pe.MapNode(interface=fsl.Cluster(), # name='cluster', # iterfield=['in_file', 'volume', 'dlh']) # #output of cluster index (in size order) # cluster.inputs.out_index_file = True # #thresholded image # cluster.inputs.out_threshold_file = True # #local maxima text file # #defines the cluster cordinates # cluster.inputs.out_localmax_txt_file = True cluster_imports = ['import os', 'import re', 'import subprocess'] cluster = pe.MapNode(util.Function(input_names=['in_file', 'volume', 'dlh', 'threshold', 'pthreshold', 'parameters'], output_names=['index_file', 'threshold_file', 'localmax_txt_file'], function=call_cluster, imports=cluster_imports), name='cluster', iterfield=['in_file', 'volume', 'dlh']) # max and minimum intensity values image_stats = pe.MapNode(interface=fsl.ImageStats(), name='image_stats', iterfield=['in_file']) image_stats.inputs.op_string = '-R' # create tuple of z_threshold and max intensity value of threshold file create_tuple = pe.MapNode(util.Function(input_names=['infile_a', 'infile_b'], output_names=['out_file'], function=get_tuple), name='create_tuple', iterfield=['infile_b']) # colour activation overlaying overlay = pe.MapNode(interface=fsl.Overlay(), name='overlay', iterfield=['stat_image', 'stat_thresh']) overlay.inputs.transparency = True overlay.inputs.auto_thresh_bg = True overlay.inputs.out_type = 'float' # colour rendering slicer = pe.MapNode(interface=fsl.Slicer(), name='slicer', iterfield=['in_file']) # set max picture width slicer.inputs.image_width = 750 # set output all axial slices into one picture slicer.inputs.all_axial = True # function mapnode to get the standard fsl brain image based on parameters # as FSLDIR,MNI and voxel size get_bg_imports = ['import os', 'import nibabel as nib'] get_backgroundimage = pe.MapNode(util.Function(input_names=['in_file', 'file_parameters'], output_names=['out_file'], function=get_standard_background_img, imports=get_bg_imports), name='get_bckgrndimg1', iterfield=['in_file']) # function node to get the standard fsl brain image # outputs single file get_backgroundimage2 = pe.Node(util.Function(input_names=['in_file', 'file_parameters'], output_names=['out_file'], function=get_standard_background_img, imports=get_bg_imports), name='get_backgrndimg2') # connections easy_thresh.connect(inputnode, 'z_stats', smooth_estimate, 'zstat_file') easy_thresh.connect(inputnode, 'merge_mask', smooth_estimate, 'mask_file') easy_thresh.connect(inputnode, 'z_stats', zstat_mask, 'in_file') easy_thresh.connect(inputnode, 'merge_mask', zstat_mask, 'operand_files') easy_thresh.connect(zstat_mask, 'out_file', get_backgroundimage, 'in_file') easy_thresh.connect(inputnode, 'parameters', get_backgroundimage, 'file_parameters') easy_thresh.connect(get_backgroundimage, 'out_file', copy_geometry, 'infile_a') easy_thresh.connect(zstat_mask, 'out_file', copy_geometry, 'infile_b') easy_thresh.connect(copy_geometry, 'out_file', cluster, 'in_file') easy_thresh.connect(inputnode, 'z_threshold', cluster, 'threshold') easy_thresh.connect(inputnode, 'p_threshold', cluster, 'pthreshold') easy_thresh.connect(smooth_estimate, 'volume', cluster, 'volume') easy_thresh.connect(smooth_estimate, 'dlh', cluster, 'dlh') easy_thresh.connect(inputnode, 'parameters', cluster, 'parameters') easy_thresh.connect(cluster, 'threshold_file', image_stats, 'in_file') easy_thresh.connect(image_stats, 'out_stat', create_tuple, 'infile_b') easy_thresh.connect(inputnode, 'z_threshold', create_tuple, 'infile_a') easy_thresh.connect(cluster, 'threshold_file', overlay, 'stat_image') easy_thresh.connect(create_tuple, 'out_file', overlay, 'stat_thresh') easy_thresh.connect(inputnode, 'merge_mask', get_backgroundimage2, 'in_file' ) easy_thresh.connect(inputnode, 'parameters', get_backgroundimage2, 'file_parameters') easy_thresh.connect(get_backgroundimage2, 'out_file', overlay, 'background_image') easy_thresh.connect(overlay, 'out_file', slicer, 'in_file') easy_thresh.connect(cluster, 'threshold_file', outputnode, 'cluster_threshold') easy_thresh.connect(cluster, 'index_file', outputnode, 'cluster_index') easy_thresh.connect(cluster, 'localmax_txt_file', outputnode, 'cluster_localmax_txt') easy_thresh.connect(overlay, 'out_file', outputnode, 'overlay_threshold') easy_thresh.connect(slicer, 'out_file', outputnode, 'rendered_image') return easy_thresh
def call_cluster(in_file, volume, dlh, threshold, pthreshold, parameters): out_name = re.match(r'z(\w)*stat(\d)+', os.path.basename(in_file)) filename, ext = os.path.splitext(os.path.basename(in_file)) ext= os.path.splitext(filename)[1] + ext filename = os.path.splitext(filename)[0] if out_name: out_name= else: out_name = filename FSLDIR = parameters[0] index_file = os.path.join(os.getcwd(), 'cluster_mask_' + out_name + ext) threshold_file = os.path.join(os.getcwd(), 'thresh_' + out_name + ext) localmax_txt_file = os.path.join(os.getcwd(), 'cluster_' + out_name + '.txt') cmd_path = os.path.join(FSLDIR, 'bin/cluster') f = open(localmax_txt_file,'wb') cmd = subprocess.Popen([ cmd_path, '--dlh=' + str(dlh), '--in=' + in_file, '--oindex=' + index_file, '--othresh=' + threshold_file, '--pthresh=' + str(pthreshold), '--thresh=' + str(threshold), '--volume=' + str(volume)], stdout= f) stdout_value, stderr_value = cmd.communicate() f.close() return index_file, threshold_file, localmax_txt_file
[docs]def copy_geom(infile_a, infile_b): """ Method to call fsl fslcpgeom command to copy certain parts of the header information (image dimensions, voxel dimensions, voxel dimensions units string, image orientation/origin or qform/sform info) from one image to another Parameters ----------- infile_a : nifti file input volume from which the geometry is copied from infile_b : nifti file input volume from which the geometry is copied to Returns ------- out_file : nifti file Input volume infile_b with modified geometry information in the header. Raises ------ Exception If fslcpgeom fails """ try: out_file = infile_b cmd = ['fslcpgeom', infile_a, out_file] retcode = subprocess.check_output(cmd) return out_file except Exception: raise Exception("Error while using fslcpgeom to copy geometry")
[docs]def get_standard_background_img(in_file, file_parameters): """ Method to get the standard brain image from FSL standard data directory Parameters ---------- in_file : nifti file Merged 4D Zmap volume file_parameters : tuple Value FSLDIR and MNI from config file Returns ------- standard_path : string Standard FSL Image file path Raises ------ Exception If nibabel cannot load the input nifti volume """ try: img = nib.load(in_file) hdr = img.header group_mm = int(hdr.get_zooms()[2]) FSLDIR, MNI = file_parameters standard_path = \ os.path.join(FSLDIR, 'data/standard/', '{0}_T1_{1}mm_brain.nii.gz'.format(MNI, group_mm)) return os.path.abspath(standard_path) except Exception: raise Exception("Error while loading background image")
[docs]def get_tuple(infile_a, infile_b): """ Simple method to return tuple of z_threhsold maximum intensity values of Zstatistic image for input to the overlay. Parameters ---------- z_theshold : float z threshold value intensity_stat : tuple of float values minimum and maximum intensity values Returns ------- img_min_max : tuple (float) tuple of zthreshold and maximum intensity value of z statistic image """ out_file = (infile_a, infile_b[1]) return out_file