diff --git a/brainiak/utils/fmrisim.py b/brainiak/utils/fmrisim.py index c93cbcee7..1e9d07833 100644 --- a/brainiak/utils/fmrisim.py +++ b/brainiak/utils/fmrisim.py @@ -17,37 +17,58 @@ Simulate fMRI data for a single subject. This code provides a set of functions necessary to produce realistic -simulations of neural data. +simulations of fMRI data. There are two main steps: characterizing the +signal and generating the noise model, which are then combined to simulate +brain data. Tools are included to support the creation of different types +of signal, such as region specific differences in univariate +activity. To create the noise model the parameters can either be set +manually or can be estimated from real fMRI data with reasonable accuracy ( +works best when fMRI data has not been preprocessed) -Steps: +Functions: generate_signal -Specify the volume (or volumes) which represent the signal in the neural data. +Create a volume with activity, of a specified shape and either multivariate +or univariate pattern, in a specific region to represent the signal in the +neural data. generate_stimfunction -Create a function to describe the stimulus onsets/durations +Create a timecourse of the signal activation. This can be specified using event +onsets and durations from a timing file. This is the time course before +convolution and therefore can be at any temporal precision. -export_stimfunction: +export_3_column: Generate a three column timing file that can be used with software like FSL +to represent event event onsets and duration -double_gamma_hrf -Convolve the stimulus function with the HRF to model when events are expected. +export_epoch_file: +Generate an epoch file from the time course which can be used as an input to +brainiak functions + +convolve_hrf +Convolve the signal timecourse with the HRF to model the expected evoked +activity apply_signal -Combine the volume and the HRF +Combine the signal volume with the HRF, thus giving the signal the temporal +properties of the HRF (such as smoothing and lag) calc_noise -Estimate the noise properties of a given volume +Estimate the noise properties of a given fMRI volume. Prominently, estimate +the smoothing and SFNR of the data generate_noise -Create the noise for this run. This creates temporal, task and white noise. -Various parameters can be tuned depending on need +Create the noise for this run. This creates temporal, spatial task and white +noise. Various parameters can be tuned depending on need mask_brain -Mask the volume to look like a volume. Based on MNI standard space +Create a mask volume that has similar contrast as an fMRI image. Defaults to +use an MNI grey matter atlas but any image can be supplied to create an +estimate. plot_brain -Show the brain as it unfolds over time with a given opacity. +Display the brain, timepoint by timepoint, with above threshold voxels +highlighted against the outline of the brain. Authors: @@ -62,13 +83,15 @@ import numpy as np from pkg_resources import resource_stream from scipy import stats +from scipy import signal import scipy.ndimage as ndimage __all__ = [ "generate_signal", "generate_stimfunction", - "export_stimfunction", - "double_gamma_hrf", + "export_3_column", + "export_epoch_file", + "convolve_hrf", "apply_signal", "calc_noise", "generate_noise", @@ -85,18 +108,19 @@ def _generate_feature(feature_type, thickness=1): """Generate features corresponding to signal - Generate signal in specific regions of the brain with for a single - volume. This will then be convolved with the HRF across time + Generate a single feature, that can be inserted into the signal volume. + A feature is a region of activation with a specific shape such as cube + or ring Parameters ---------- feature_type : str - What feature_type of signal is being inserted? Options are cube, - loop, cavity, sphere. + What shape signal is being inserted? Options are 'cube', + 'loop' (aka ring), 'cavity' (aka hollow sphere), 'sphere'. feature_size : int - How big is the signal? + How big is the signal in diameter? signal_magnitude : float Set the signal size, a value of 1 means the signal is one standard @@ -108,8 +132,8 @@ def _generate_feature(feature_type, Returns ---------- - 3 dimensional array - The volume representing the signal to be outputed + signal : 3 dimensional array + The volume representing the signal """ @@ -121,20 +145,12 @@ def _generate_feature(feature_type, if feature_type == 'cube': # Preset the size of the signal - signal = np.ones(np.power(feature_size, 3)) - - # Reorganize the signal into a 3d matrix - signal = signal.reshape([feature_size, - feature_size, - feature_size]) + signal = np.ones((feature_size, feature_size, feature_size)) elif feature_type == 'loop': # First make a cube of zeros - signal = np.zeros(np.power(feature_size, - 3)).reshape([feature_size, - feature_size, - feature_size]) + signal = np.zeros((feature_size, feature_size, feature_size)) # Make a mesh grid of the space seq = np.linspace(0, feature_size - 1, @@ -159,6 +175,11 @@ def _generate_feature(feature_type, # Subtract the two disks to get a loop loop = outer != inner + # Check if the loop is a disk + if np.all(inner is False): + logger.warning('Loop feature reduces to a disk because the loop ' + 'is too thick') + # If there is complete overlap then make the signal just the # outer one if np.all(loop is False): @@ -199,6 +220,11 @@ def _generate_feature(feature_type, # Subtract the two disks to get a loop signal = outer != inner + # Check if the cavity is a sphere + if np.all(inner is False): + logger.warning('Cavity feature reduces to a sphere because ' + 'the cavity is too thick') + # If there is complete overlap then make the signal just the # outer one if np.all(signal is False): @@ -212,7 +238,7 @@ def _generate_feature(feature_type, def _insert_idxs(feature_centre, feature_size, dimensions): - """Returns the indexes of where to put the signal into Volume_Static + """Returns the indices of where to put the signal into the signal volume Parameters ---------- @@ -221,7 +247,7 @@ def _insert_idxs(feature_centre, feature_size, dimensions): List of coordinates for the centre location of the signal feature_size : list, int - How big is the signal. + How big is the signal's diameter. dimensions : 3 length array, int What are the dimensions of the volume you wish to create @@ -281,7 +307,7 @@ def generate_signal(dimensions, ): """Generate volume containing signal - Generate signal in specific regions of the brain with for a single + Generate signal, of a specific shape in specific regions, for a single volume. This will then be convolved with the HRF across time Parameters @@ -294,19 +320,20 @@ def generate_signal(dimensions, What are the feature_coordinates of the signal being created. Be aware of clipping: features far from the centre of the brain will be clipped. If you wish to have multiple features - then list these as an features x 3 array. To create a signal of - a specific shape then supply all the individual - feature_coordinates and set the feature_size to 1. + then list these as a features x 3 array. To create a feature of + a unique shape then supply all the individual + feature_coordinates of the shape and set the feature_size to 1. feature_size : list, int - How big is the signal. If m=1 then only one value is accepted, - if m>1 then either one value must be supplied or m values + How big is the signal. If feature_coordinates=1 then only one value is + accepted, if feature_coordinates>1 then either one value must be + supplied or m values feature_type : list, string What feature_type of signal is being inserted? Options are cube, - loop, cavity, sphere. If features = 1 then - only one value is accepted, if features > 1 then either one value - must be supplied or m values + loop, cavity, sphere. If feature_coordinates = 1 then + only one value is accepted, if feature_coordinates > 1 then either + one value must be supplied or m values signal_magnitude : list, float What is the (average) magnitude of the signal being generated? A @@ -314,18 +341,19 @@ def generate_signal(dimensions, noise signal_constant : list, bool - Is the signal constant or is it a random pattern (with the same - average magnitude) + Is the signal constant across the feature (for univariate activity) + or is it a random pattern of a given magnitude across the feature (for + multivariate activity) Returns ---------- - volume_static : 3 dimensional array, float + volume_signal : 3 dimensional array, float Creates a single volume containing the signal """ # Preset the volume - volume_static = np.zeros(dimensions) + volume_signal = np.zeros(dimensions) feature_quantity = round(feature_coordinates.shape[0]) @@ -342,7 +370,7 @@ def generate_signal(dimensions, signal_magnitude = signal_magnitude * feature_quantity # Iterate through the signals and insert in the data - for signal_counter in list(range(0, feature_quantity)): + for signal_counter in range(feature_quantity): # What is the centre of this signal if len(feature_size) > 1: @@ -369,10 +397,10 @@ def generate_signal(dimensions, dimensions) # Insert the signal into the Volume - volume_static[x_idx[0]:x_idx[1], y_idx[0]:y_idx[1], z_idx[0]:z_idx[ + volume_signal[x_idx[0]:x_idx[1], y_idx[0]:y_idx[1], z_idx[0]:z_idx[ 1]] = signal - return volume_static + return volume_signal def generate_stimfunction(onsets, @@ -382,43 +410,47 @@ def generate_stimfunction(onsets, timing_file=None, temporal_resolution=1000.0, ): - """Return the function for the onset of events + """Return the function for the timecourse events When do stimuli onset, how long for and to what extent should you resolve the fMRI time course. There are two ways to create this, either by supplying onset, duration and weight information or by supplying a - timing file (in the three column format) + timing file (in the three column format used by FSL). Parameters ---------- - onsets : list, int - What are the timestamps for when an event you want to generate - onsets? + onsets : list, int + What are the timestamps (in s) for when an event you want to + generate onsets? - event_durations : list, int - What are the durations of the events you want to generate? If - there is only one value then this will be assigned to all onsets + event_durations : list, int + What are the durations (in s) of the events you want to + generate? If there is only one value then this will be assigned + to all onsets - total_time : int - How long is the experiment in total. + total_time : int + How long (in s) is the experiment in total. - weights : list, float - How large is the box car for each of the onsets. If there is - only one value then this will be assigned to all onsets + weights : list, float + What is the weight for each event (how high is the box car)? If + there is only one value then this will be assigned to all onsets - timing_file : string - The filename (with path) to a three column timing file (FSL) to - make the events. Still requires tr_duration and total_time + timing_file : string + The filename (with path) to a three column timing file (FSL) to + make the events. Still requires total_time to work - temporal_resolution : float - How many elements per second are you modeling for the stim function + temporal_resolution : float + How many elements per second are you modeling for the + timecourse. This is useful when you want to model the HRF at an + arbitrarily high resolution (and then downsample to your TR later). Returns ---------- - Iterable[float] - The time course of stimulus related activation + stim_function : 1 by timepoint array, float + The time course of stimulus evoked activation. This has a temporal + resolution of temporal resolution / 1.0 elements per second """ @@ -449,9 +481,14 @@ def generate_stimfunction(onsets, if len(weights) == 1: weights = weights * len(onsets) + # Check files + if np.max(onsets) > total_time: + print('Onsets outside of range of total time. Aborting') + exit() + # Generate the time course as empty, each element is a millisecond by # default - stimfunction = [0] * int(round(total_time * temporal_resolution)) + stimfunction = np.zeros((int(round(total_time * temporal_resolution)), 1)) # Cycle through the onsets for onset_counter in list(range(len(onsets))): @@ -460,63 +497,65 @@ def generate_stimfunction(onsets, # Adjust for the resolution offset_idx = int(np.floor((onsets[onset_counter] + event_durations[ - onset_counter])) * temporal_resolution) + onset_counter]) * temporal_resolution)) - # For the appropriate number of indexes and duration, make this value 1 - idx_n = round(event_durations[onset_counter] * temporal_resolution) - stimfunction[onset_idx:offset_idx] = [weights[onset_counter]] * idx_n + # Store the weights + stimfunction[onset_idx:offset_idx, 0] = [weights[onset_counter]] # Shorten the data if it's too long - if len(stimfunction) > total_time * temporal_resolution: - stimfunction = stimfunction[0:int(total_time * temporal_resolution)] + if stimfunction.shape[0] > total_time * temporal_resolution: + stimfunction = stimfunction[0:int(total_time * temporal_resolution), 0] return stimfunction -def export_stimfunction(stimfunction, - filename, - temporal_resolution=1000.0 - ): - """ Output a tab separated timing file +def export_3_column(stimfunction, + filename, + temporal_resolution=1000.0 + ): + """ Output a tab separated three column timing file This produces a three column tab separated text file, with the three - columns representing onset time, event duration and weight, respectively - - Useful if you want to run the simulated data through FEAT analyses + columns representing onset time (s), event duration (s) and weight, + respectively. Useful if you want to run the simulated data through FEAT + analyses. In a way, this is the reverse of generate_stimfunction Parameters ---------- - stimfunction : list - The stimulus function describing the time course of events - filename : str - The name of the three column text file to be output + stimfunction : timepoint by 1 array + The stimulus function describing the time course of events. For + instance output from generate_stimfunction. + + filename : str + The name of the three column text file to be output - temporal_resolution : float - How many elements per second are you modeling for the stim function + temporal_resolution : float + How many elements per second are you modeling with the + stimfunction? """ # Iterate through the stim function stim_counter = 0 event_counter = 0 - while stim_counter < len(stimfunction): + while stim_counter < stimfunction.shape[0]: # Is it an event? - if stimfunction[stim_counter] != 0: + if stimfunction[stim_counter, 0] != 0: # When did the event start? event_onset = str(stim_counter / temporal_resolution) # The weight of the stimulus - weight = str(stimfunction[stim_counter]) + weight = str(stimfunction[stim_counter, 0]) # Reset event_duration = 0 # Is the event still ongoing? - while stimfunction[stim_counter] != 0 & stim_counter <= len( - stimfunction): + while stimfunction[stim_counter, 0] != 0 & stim_counter <= \ + stimfunction.shape[0]: # Add one millisecond to each duration event_duration = event_duration + 1 @@ -539,56 +578,155 @@ def export_stimfunction(stimfunction, stim_counter = stim_counter + 1 -def double_gamma_hrf(stimfunction, - tr_duration, - response_delay=6, - undershoot_delay=12, - response_dispersion=0.9, - undershoot_dispersion=0.9, - response_scale=1, - undershoot_scale=0.035, - scale_function=1, - temporal_resolution=1000.0, - ): - """Return a double gamma HRF +def export_epoch_file(stimfunction, + filename, + tr_duration, + temporal_resolution=1000.0 + ): + """ Output an epoch file, necessary for some inputs into brainiak + + This takes in the time course of stimulus events and outputs the epoch + file used in Brainiak. The epoch file is a way to structure the timing + information in fMRI that allows you to flexibly input different stimulus + sequences. This is a list with each entry a 3d matrix corresponding to a + participant. The dimensions of the 3d matrix are condition by epoch by time Parameters ---------- - stimfunction : list, bool - What is the time course of events to be modelled in this - experiment - tr_duration : float + stimfunction : list of timepoint by condition arrays + The stimulus function describing the time course of events. Each + list entry is from a different participant, each row is a different + timepoint (with the given temporal precision), each column is a + different condition. export_epoch_file is looking for differences in + the value of stimfunction to identify the start and end of an + epoch. If epochs in stimfunction are coded with the same weight and + there is no time between blocks then export_epoch_file won't be able to + label them as different epochs + + filename : str + The name of the three column text file to be output + + tr_duration : float + How long is each TR in seconds + + temporal_resolution : float + How many elements per second are you modeling with the + stimfunction? - response_delay : float - How many seconds until the peak of the HRF + """ + + # Cycle through the participants, different entries in the list + epoch_file = [0] * len(stimfunction) + for participant_counter in range(len(stimfunction)): + + # What is the time course for the participant (binarized) + stimfunction_ppt = np.abs(stimfunction[participant_counter]) > 0 + + # Cycle through conditions + conditions = stimfunction_ppt.shape[1] + for condition_counter in range(conditions): + + # Down sample the stim function + stride = tr_duration * temporal_resolution + stimfunction_temp = stimfunction_ppt[:, condition_counter] + stimfunction_temp = stimfunction_temp[::int(stride)] + + if condition_counter == 0: + # Calculates the number of event onsets (max of all + # conditions). This uses changes in value to reflect + # different epochs. This might be false in some cases (the + # weight is supposed to unfold over an epoch or there is no + # break between identically weighted epochs). In such cases + # this will not work + weight_change = (np.diff(stimfunction_temp, 1, 0) != 0) + epochs = int(np.max(np.sum(weight_change, 0)) / 2) + + # Get other information + trs = stimfunction_temp.shape[0] + + # Make a timing file for this participant + epoch_file[participant_counter] = np.zeros((conditions, + epochs, trs)) + + epoch_counter = 0 + tr_counter = 0 + while tr_counter < stimfunction_temp.shape[0]: + + # Is it an event? + if stimfunction_temp[tr_counter] == 1: + + # Add a one for this TR + epoch_file[participant_counter][condition_counter, + epoch_counter, + tr_counter] = 1 + + # Find the next non event value + end_idx = np.where(stimfunction_temp[tr_counter:] == 0)[ + 0][0] + tr_idxs = list(range(tr_counter, tr_counter + end_idx)) + + # Add ones to all the trs within this event time frame + epoch_file[participant_counter][condition_counter, + epoch_counter, + tr_idxs] = 1 + + # Start from this index + tr_counter += end_idx + + # Increment + epoch_counter += 1 + + # Increment the counter + tr_counter += 1 + + # Save the file + np.save(filename, epoch_file) + + +def _double_gamma_hrf(response_delay=6, + undershoot_delay=12, + response_dispersion=0.9, + undershoot_dispersion=0.9, + response_scale=1, + undershoot_scale=0.035, + temporal_resolution=1000.0, + ): + """Create the double gamma HRF with the timecourse evoked activity. + Default values are based on Glover, 1999 and Walvaert, Durnez, + Moerkerke, Verdoolaege and Rosseel, 2011 + + Parameters + ---------- - undershoot_delay : float - How many seconds until the trough of the HRF + response_delay : float + How many seconds until the peak of the HRF - response_dispersion : float - How wide is the rising peak dispersion + undershoot_delay : float + How many seconds until the trough of the HRF - undershoot_dispersion : float - How wide is the undershoot dispersion + response_dispersion : float + How wide is the rising peak dispersion - response_scale : float - How big is the response relative to the peak + undershoot_dispersion : float + How wide is the undershoot dispersion - undershoot_scale :float - How big is the undershoot relative to the trough + response_scale : float + How big is the response relative to the peak - scale_function : bool - Do you want to scale the function to a range of 1 + undershoot_scale :float + How big is the undershoot relative to the trough - temporal_resolution : float - How many elements per second are you modeling for the stim function + scale_function : bool + Do you want to scale the function to a range of 1 + + temporal_resolution : float + How many elements per second are you modeling for the stimfunction Returns ---------- - one dimensional array - The time course of the HRF convolved with the stimulus function - + hrf : multi dimensional array + A double gamma HRF to be used for convolution. """ @@ -624,57 +762,153 @@ def double_gamma_hrf(stimfunction, # For this time point find the value of the HRF hrf[hrf_counter] = response_model - undershoot_model - # Convolve the hrf that was created with the boxcar input - signal_function = np.convolve(stimfunction, hrf) + return hrf + + +def convolve_hrf(stimfunction, + tr_duration, + hrf_type='double_gamma', + scale_function=True, + temporal_resolution=1000.0, + ): + """ Convolve the specified hrf with the timecourse + + Parameters + ---------- + + stimfunction : timepoint by timecourse array + What is the time course of events to be modelled in this + experiment. This can specify one or more timecourses of events. + The events can be weighted or binary + + tr_duration : float + How long (in s) between each volume onset + + hrf_type : str or list + Takes in a string describing the hrf that ought to be created. + Can instead take in a vector describing the HRF as it was + specified by any function. The default is 'double_gamma' in which + an initial rise and an undershoot are modelled. + + scale_function : bool + Do you want to scale the function to a range of 1 + + temporal_resolution : float + How many elements per second are you modeling for the stimfunction + Returns + ---------- + + signal_function : timepoint by timecourse array + The time course of the HRF convolved with the stimulus function. + This can have multiple time courses specified as different + columns in this array. + + """ + + # Generate the hrf to use in the convolution + if hrf_type == 'double_gamma': + hrf = _double_gamma_hrf(temporal_resolution=temporal_resolution) + elif isinstance(hrf_type, list): + hrf = hrf_type + + # How many timecourses are there + list_num = stimfunction.shape[1] + + # Create signal functions for each list in the stimfunction + for list_counter in range(list_num): + + # Take the stim function + stimfunction_temp = stimfunction[:, list_counter] + + signal_function_temp = np.convolve(stimfunction_temp, hrf) - # Decimate the signal function so that it only has one element per TR - decimate_interval = int(tr_duration * temporal_resolution) - signal_function = signal_function[0::decimate_interval] + # Decimate the signal function so that it only has one element per TR + decimate_interval = int(tr_duration * temporal_resolution) + signal_function_temp = signal_function_temp[0::decimate_interval] - # Cut off the HRF - signal_function = signal_function[0:int((len(stimfunction) / - tr_duration) / - temporal_resolution)] + # Cut off the HRF + last_timepoint = stimfunction_temp.shape[0] / tr_duration + last_timepoint /= temporal_resolution + signal_function_temp = signal_function_temp[0:int(last_timepoint)] - # Scale the function so that the peak response is 1 - if scale_function == 1: - signal_function = signal_function / np.max(signal_function) + # Scale the function so that the peak response is 1 + if scale_function: + signal_function_temp = signal_function_temp / np.max( + signal_function_temp) + + # Add this function to the stack + if list_counter == 0: + signal_function = np.zeros((len(signal_function_temp), list_num)) + + signal_function[:, list_counter] = signal_function_temp return signal_function def apply_signal(signal_function, - volume_static, + volume_signal, ): - """Apply the convolution and stimfunction + """Combine the signal volume with its timecourse Apply the convolution of the HRF and stimulus time course to the volume. Parameters ---------- - signal_function : list, float - The one dimensional timecourse of the signal over time. - Found by convolving the HRF with the stimulus time course. - volume_static : multi dimensional array, float + signal_function : timepoint by timecourse array, float + The timecourse of the signal over time. If there is only one column + then the same timecourse is applied to all non-zero voxels in + volume_signal. If there is more than one column then each column is + paired with a non-zero voxel in the volume_signal (a 3d numpy array + generated in generate_signal). + + volume_signal : multi dimensional array, float + The volume containing the signal to be convolved with the same + dimensions as the output volume. The elements in volume_signal + indicate how strong each signal in signal_function are modulated by + in the output volume Returns ---------- - multidimensional array, float - Generates the spatial noise volume for these parameters """ + signal : multidimensional array, float + The convolved signal volume with the same 3d as volume signal and + the same 4th dimension as signal_function + + """ + + # How many timecourses are there within the signal_function + timepoints = signal_function.shape[0] + timecourses = signal_function.shape[1] # Preset volume - signal = np.ndarray([volume_static.shape[0], volume_static.shape[ - 1], volume_static.shape[2], len(signal_function)]) + signal = np.zeros([volume_signal.shape[0], volume_signal.shape[ + 1], volume_signal.shape[2], timepoints]) + + # Find all the non-zero voxels in the brain + idxs = np.where(volume_signal != 0) + if timecourses == 1: + # If there is only one time course supplied then duplicate it for + # every voxel + signal_function = np.matlib.repmat(signal_function, 1, len(idxs[0])) + + elif len(idxs[0]) != timecourses: + raise IndexError('The number of non-zero voxels in the volume and ' + 'the number of timecourses does not match. Aborting') + + # For each coordinate with a non zero voxel, fill in the timecourse for + # that voxel + for idx_counter in range(len(idxs[0])): + x = idxs[0][idx_counter] + y = idxs[1][idx_counter] + z = idxs[2][idx_counter] + + # Pull out the function for this voxel + signal_function_temp = signal_function[:, idx_counter] - # Iterate through the brain, multiplying the volume by the HRF - for tr_counter in list(range(0, len(signal_function))): - signal[0:volume_static.shape[0], - 0:volume_static.shape[1], - 0:volume_static.shape[2], - tr_counter] = signal_function[tr_counter] * volume_static + # Multiply the voxel value by the function timecourse + signal[x, y, z, :] = volume_signal[x, y, z] * signal_function_temp return signal @@ -684,25 +918,27 @@ def _calc_fwhm(volume, voxel_size=[1.0, 1.0, 1.0], ): """ Calculate the FWHM of a volume - Takes in a 3d volume and mask and outputs the FWHM (mm) of this - volume for the non-masked voxels + Estimates the FWHM (mm) of a volume's non-masked voxels Parameters ---------- + volume : 3 dimensional array - Functional data to have the FWHM measured. + Functional data to have the FWHM measured. mask : 3 dimensional array - A mask of the voxels to have the FWHM measured from + A binary mask of the brain voxels in volume voxel_size : length 3 list, float - Millimeters per voxel for x, y and z. + Millimeters per voxel for x, y and z. Returns ------- - float, list - Returns the FWHM of each TR in mm""" + fwhm : float, list + Returns the FWHM of each TR in mm + + """ # What are the dimensions of the volume dimensions = volume.shape @@ -788,118 +1024,165 @@ def _calc_fwhm(volume, def _calc_sfnr(volume, mask, ): - """ Calculate the SFNR of a volume - This takes the middle of the volume and averages the signal within the - brain and compares to the temporal standard deviation in the voxels - outside the brain. + """ Calculate the the SFNR of a volume + Calculates the Signal to Fluctuation Noise Ratio, the mean divided + by the detrended standard deviation of each brain voxel. Based on + Friedman and Glover, 2006 Parameters ---------- + volume : 4d array, float - Take a volume time series to extract the middle slice from the middle TR + Take a volume time series - mask : 4d array, float - A mask the same size as the volume, specifying the mask (values=0 -> 1) + mask : 3d array, binary + A binary mask the same size as the volume Returns ------- - float
 - The sd of the temporal variability of brain voxels.

 + snr : float
 + The SFNR of the volume - float
 - The sfnr of the volume (mean brain activity divided by
 temporal - variability in the average non brain voxels)

 + """ - float
 - What is the max activity measured here, a point of reference for masking + # Make a matrix of brain voxels by time + brain_voxels = volume[mask > 0] - """ + # Take the means of each voxel over time + mean_voxels = np.nanmean(brain_voxels, 1) + + # Detrend (second order polynomial) the voxels over time and then + # calculate the standard deviation. + order = 2 + seq = np.linspace(1, brain_voxels.shape[1], brain_voxels.shape[1]) + detrend_poly = np.polyfit(seq, brain_voxels.transpose(), order) + + # Detrend for each voxel + detrend_voxels = np.zeros(brain_voxels.shape) + for voxel in range(brain_voxels.shape[0]): + trend = detrend_poly[0, voxel] * seq ** 2 + detrend_poly[1, voxel] * \ + seq + detrend_poly[2, voxel] + detrend_voxels[voxel, :] = brain_voxels[voxel, :] - trend - # What are the midpoints to be extracted - mid_x_idx = int(np.ceil(volume.shape[0] / 2)) - mid_tr_idx = int(np.ceil(volume.shape[3] / 2)) + std_voxels = np.nanstd(detrend_voxels, 1) - # Pull out the slices - slice_volume = volume[mid_x_idx, :, :, mid_tr_idx] - slice_mask = mask[mid_x_idx, :, :, mid_tr_idx] + # Calculate the sfnr of all voxels across the brain + sfnr_voxels = mean_voxels / std_voxels - # Divide by the mask (if the middle slice was low in grey matter mass - # then you would have a lower mean signal by default) - mean_signal = (slice_volume[slice_mask > 0]).mean() + # Return the average sfnr + return np.mean(sfnr_voxels) - # What are the brain and non-brain voxels - brain_voxels = volume[mask[:, :, :, 0] > 0] - mask_voxels = volume[mask[:, :, :, 0] == 0] - # What is the noise in the masked voxels - mask_noise = np.std(mask_voxels, 1).mean() +def _calc_snr(volume, + mask, + tr=None, + ): + """ Calculate the the SNR of a volume + Calculates the Signal to Noise Ratio, the mean of brain voxels + divided by the standard deviation across non-brain voxels. Specify a TR + value to calculate the mean and standard deviation for that TR. To + calculate the standard deviation this subtracts any baseline structure + in the non-brain voxels, hence getting at deviations due to the system + noise and not something like high baseline values in non-brain parts of + the body. + + Parameters + ---------- - # Assume the mask noise is a combination of random + drift. Thus average - # all masked voxels and find the variation over time, this is the - # variance due to drift. - drift_noise = np.mean(mask_voxels, 0).std() + volume : 4d array, float + Take a volume time series - # Subtract the drift variance from total mask noise to find the system - # noise - background_noise = np.sqrt(mask_noise ** 2 - drift_noise ** 2) + mask : 3d array, binary + A binary mask the same size as the volume - # Find the noise to brain voxels - temporal_noise = np.std(brain_voxels, 1).mean() + tr : int + Integer specifying TR to calculate the SNR for - # Convert temporal noise into percent signal change - temporal_noise = temporal_noise / mean_signal * 100 + Returns + ------- - # Calculate sfnr - sfnr = mean_signal / background_noise + snr : float
 + The SNR of the volume - # Convert from memmap - sfnr = float(sfnr) + """ - # What is the max activation of this volume - max_activity = volume.max() + # If no TR is specified then take the middle one + if tr is None: + tr = int(np.ceil(volume.shape[3] / 2)) - return temporal_noise, sfnr, max_activity + # Make a matrix of brain and non_brain voxels by time + brain_voxels = volume[mask > 0] + nonbrain_voxels = volume[:, :, :, tr][mask == 0] + + # Find the mean of the non_brain voxels (deals with structure that may + # exist outside of the mask) + nonbrain_voxels_mean = np.mean(volume[mask == 0], 1) + + # Take the means of each voxel over time + mean_voxels = np.nanmean(brain_voxels) + std_voxels = np.nanstd(nonbrain_voxels - nonbrain_voxels_mean) + + # Return the snr + return mean_voxels / std_voxels def _calc_temporal_noise(volume, mask, auto_reg_order=1, ): - """ Calculate the mix of autoregressive and drift noise. + """ Calculate the the temporal noise of a volume + This calculates the variability of the volume over time and the + proportion of variance over time that is due to autoregression and how + much is due to scanner drift. Parameters ---------- - volume : 4d masarray, float - Input data to be calculate the autoregression - mask : 4d array, float - What voxels of the input are within the brain + volume : 4d array, float + Take a volume time series to extract the middle slice from the + middle TR + + mask : 3d array, binary + A binary mask the same size as the volume auto_reg_order : int - What order of the autoregression do you want to pull out + What order of the autoregression do you want to pull out + Returns ------- - float - A sigma of the autoregression in the data - float - Sigma of the drift in the data + + sfnr : float
 + The SFNR of the volume (mean brain activity divided by
 temporal + variability in the averaged non brain voxels)

 + + auto_reg_sigma : float + A sigma of the autoregression in the data + + drift_sigma : float + Sigma of the drift in the data """ - # Calculate the time course - timecourse = np.mean(volume[mask[:, :, :, 0] > 0], 0) + # Calculate sfnr and convert from memmap + sfnr = _calc_sfnr(volume, + mask, + ) + + # Calculate the time course of voxels within the brain + timecourse = np.mean(volume[mask > 0], 0) + demeaned_timecourse = timecourse-timecourse.mean() # Pull out the AR values (depends on order) - auto_reg_sigma = ar.AR_est_YW(timecourse, auto_reg_order)[1] - auto_reg_sigma = np.sqrt(auto_reg_sigma) + auto_reg_sigma = ar.AR_est_YW(demeaned_timecourse, auto_reg_order) + auto_reg_sigma = np.sqrt(auto_reg_sigma[1]) # What is the size of the change in the time course drift_sigma = timecourse.std().tolist() - return auto_reg_sigma, drift_sigma + return sfnr, auto_reg_sigma, drift_sigma def calc_noise(volume, @@ -914,113 +1197,165 @@ def calc_noise(volume, Parameters ---------- + volume : 4d numpy array, float Take in a functional volume (either the file name or the numpy array) to be used to estimate the noise properties of this - mask : 4d numpy array, float - The mask to be used, the same size as the volume + mask : 3d numpy array, binary + A binary mask of the brain, the same size as the volume Returns ------- - dict - Return a dictionary of the calculated noise parameters of the provided - dataset + noise_dict : dict + Return a dictionary of the calculated noise parameters of the provided + dataset """ - # Preset - - # Create the mask + # Create the mask if not supplied and set the mask size if mask is None: - mask = np.ones(volume.shape) + mask = np.ones(volume.shape[:-1]) - # Update noise dict + # Update noise dict if it is not yet created if noise_dict is None: noise_dict = {'voxel_size': [1.0, 1.0, 1.0]} elif 'voxel_size' not in noise_dict: noise_dict['voxel_size'] = [1.0, 1.0, 1.0] + # What is the max activation of the mean of this voxel (allows you to + # convert between the mask and the mean of the brain volume) + noise_dict['max_activity'] = np.nanmax(np.mean(volume, 3)) + # Since you are deriving the 'true' values then you want your noise to # be set to that level - # Calculate the temporal_noise noise and SFNR of the volume - noise_dict['temporal_noise'], noise_dict['sfnr'], noise_dict[ - 'max_activity'] = _calc_sfnr(volume, mask) + # Calculate the temporal variability of the volume + sfnr, auto_reg, drift = _calc_temporal_noise(volume, mask) + noise_dict['sfnr'] = sfnr # Calculate the fwhm on a subset of volumes - if volume.shape[3] > 100: # Take only 100 shuffled TRs - trs = np.arange(volume.shape[3]) - np.random.shuffle(trs) - trs = trs[0:100] + trs = np.random.choice(volume.shape[3], size=100, replace=False) else: trs = list(range(0, volume.shape[3])) # Go through the trs and pull out the fwhm fwhm = [0] * len(trs) - for tr in list(range(0, len(trs))): + for tr in range(len(trs)): fwhm[tr] = _calc_fwhm(volume[:, :, :, trs[tr]], - mask[:, :, :, trs[tr]], + mask, noise_dict['voxel_size'], ) # Keep only the mean noise_dict['fwhm'] = np.mean(fwhm) - - # Calculate the autoregressive and drift noise - auto_reg_sigma, drift_sigma = _calc_temporal_noise(volume, mask) - - # Calibrate for how sigma is originally calculated - auto_reg_sigma = auto_reg_sigma / noise_dict['temporal_noise'] + noise_dict['snr'] = _calc_snr(volume, + mask, + ) # Total temporal noise, since these values only make sense relatively - total_temporal_noise = auto_reg_sigma + drift_sigma + total_temporal_noise = auto_reg + drift # What proportion of noise is accounted for by these variables? - noise_dict['auto_reg_sigma'] = auto_reg_sigma / total_temporal_noise - noise_dict['drift_sigma'] = drift_sigma / total_temporal_noise + noise_dict['auto_reg_sigma'] = auto_reg / total_temporal_noise + noise_dict['drift_sigma'] = drift / total_temporal_noise # Return the noise dictionary return noise_dict def _generate_noise_system(dimensions_tr, + spatial_sd, + temporal_sd, + spatial_noise_type='exponential', + temporal_noise_type='exponential', ): """Generate the scanner noise - Generate the noise that is typical of a scanner. This is comprised - of two types of noise, Rician and Gaussian + Generate system noise, either rician or exponential, for the scanner. + Low SNR scans tend to have rician noise whereas high SNR scans (>30) are + better modelled by exponential noise. Generates a distribution with a SD + of 1. Parameters ---------- + dimensions_tr : n length array, int What are the dimensions of the volume you wish to insert noise into. This can be a volume of any size + noise_type : str + String specifying the noise type. Rician is appropriate when the SNR is + low but is insufficiently skewed to appropriately model high SNR data. + Returns ---------- - system_noise : multidimensional array, float - Create a volume with system noise + system_noise : multidimensional array, float + Create a volume with system noise + + """ + def generate_noise_volume(dimensions, + noise_type, + ): + + if noise_type == 'rician': + # Generate the Rician noise (has an SD of 1) + noise = stats.rice.rvs(b=0, loc=0, scale=1.527, size=dimensions) + elif noise_type == 'exponential': + # Make an exponential distribution (has an SD of 1) + noise = stats.expon.rvs(0, scale=1, size=dimensions) + elif noise_type == 'gaussian': + noise = np.random.randn(np.prod(dimensions)).reshape(dimensions) + + # Return the noise + return noise + + # Get just the xyz coordinates + dimensions = np.asarray([dimensions_tr[0], + dimensions_tr[1], + dimensions_tr[2], + 1]) + + # Generate noise + spatial_noise = generate_noise_volume(dimensions, spatial_noise_type) + temporal_noise = generate_noise_volume(dimensions_tr, temporal_noise_type) + + # Since you are combining spatial and temporal noise, you need to + # subtract the variance of the two to get the spatial sd + if spatial_sd > temporal_sd: + spatial_sd = np.sqrt(spatial_sd ** 2 - temporal_sd ** 2) + else: + # If this is below zero then all the noise will be temporal + spatial_sd = 0 + + # # Mean centre, while preserving the SD + # spatial_noise = spatial_noise - spatial_noise.mean() + + # Make the system noise have a specific spatial variability + spatial_noise *= spatial_sd + + # # Mean centre, while preserving the SD + # temporal_noise = temporal_noise - temporal_noise.mean() - """ + # Set the size of the noise + temporal_noise *= temporal_sd - # Generate the Rician noise - noise_rician = stats.rice.rvs(b=0, loc=0, scale=1.527, size=dimensions_tr) - # - # # Apply the gaussian noise - # noise_gaussian = np.random.normal(0, 1, size=dimensions_tr) - # - # # Combine these two noise types - # noise_system = noise_rician + noise_gaussian - # - # # Normalize - # noise_system = stats.zscore(noise_system) + # The mean in time of system noise needs to be zero, so subtract the + # means of the temporal noise in time and spatial noise + temporal_noise_mean = np.mean(temporal_noise, 3).reshape(dimensions[0], + dimensions[1], + dimensions[2], + 1) + temporal_noise = temporal_noise - (temporal_noise_mean - spatial_noise) - return noise_rician + # Save the size of the noise + system_noise = spatial_noise + temporal_noise + + return system_noise def _generate_noise_temporal_task(stimfunction_tr, @@ -1028,8 +1363,8 @@ def _generate_noise_temporal_task(stimfunction_tr, ): """Generate the signal dependent noise - This noise depends on things like the signal or the timing of the - experiment. + Create noise specific to the signal, for instance there is variability + in how the signal manifests on each event Parameters ---------- @@ -1043,13 +1378,14 @@ def _generate_noise_temporal_task(stimfunction_tr, Returns ---------- - one dimensional array, float - Generates the temporal task noise timecourse + noise_task : one dimensional array, float + Generates the temporal task noise timecourse """ # Make the noise to be added + stimfunction_tr = stimfunction_tr != 0 if motion_noise == 'gaussian': noise = stimfunction_tr * np.random.normal(0, 1, size=len( stimfunction_tr)) @@ -1067,44 +1403,43 @@ def _generate_noise_temporal_task(stimfunction_tr, def _generate_noise_temporal_drift(trs, tr_duration, + period=300, ): """Generate the drift noise - According to AFNI (https://afni.nimh.nih.gov/pub/dist/doc/ - program_help/3dDeconvolve.html) the appropriate order of the - polynomial to fit for temporal drift is calculated as follows + Create a sinewave, of a given period and random phase, to represent the + drift of the signal over time Parameters ---------- trs : int - How many TRs are there + How many volumes (aka TRs) are there tr_duration : float - How long is each TR + How long in seconds is each volume acqusition + period : int + How many seconds is the period of oscillation of the drift Returns ---------- - one dimensional array, float - Generates the autoregression noise timecourse + noise_drift : one dimensional array, float + The drift timecourse of activity """ - # What time points are sampled by a TR? - timepoints = list(range(0, trs * tr_duration))[::tr_duration] + # Calculate the cycles of the drift for a given function. + cycles = trs * tr_duration / period - # Calculate the coefficients of the drift for a given function - degree = round(trs * tr_duration / 150) + 1 - if degree > 50: - degree = 50 # Max out in order to avoid precision errors - coefficients = np.random.normal(0, 1, size=degree) + # Create a sine wave with a given number of cycles and random phase + timepoints = np.linspace(0, trs - 1, trs) + phaseshift = np.pi * 2 * np.random.random() + phase = (timepoints / (trs - 1) * cycles * 2 * np.pi) + phaseshift + noise_drift = np.sin(phase) - # What are the values of this drift - noise_drift = np.polyval(coefficients, timepoints) - - # Normalize + # Normalize so the sigma is 1 noise_drift = stats.zscore(noise_drift) # Return noise @@ -1113,10 +1448,12 @@ def _generate_noise_temporal_drift(trs, def _generate_noise_temporal_autoregression(timepoints, auto_reg_order=1, - auto_reg_rho=[1], + auto_reg_rho=[0.5], ): """Generate the autoregression noise + Make a slowly drifting timecourse with the given autoregression + parameters. The output should have an autoregression coefficient of 1 Parameters ---------- @@ -1130,21 +1467,26 @@ def _generate_noise_temporal_autoregression(timepoints, auto_reg_rho : float What is the scaling factor on the predictiveness of the previous - time point + time point. This value is below 1 to avoid brownian motion (and + growing variance). Values near or greater than one may produce drift or + other unwanted trends. Returns ---------- - one dimensional array, float + noise_autoregression : one dimensional array, float Generates the autoregression noise timecourse - """ + + """ if len(auto_reg_rho) == 1: auto_reg_rho = auto_reg_rho * auto_reg_order # Duplicate this so that # there is one # for each value + # Generate a random variable at each time point that is a decayed value + # of the previous time points noise_autoregression = [] - for tr_counter in list(range(0, len(timepoints))): + for tr_counter in range(len(timepoints)): if tr_counter == 0: noise_autoregression.append(np.random.normal(0, 1)) @@ -1156,13 +1498,14 @@ def _generate_noise_temporal_autoregression(timepoints, if tr_counter - pCounter >= 0: past_trs = noise_autoregression[int(tr_counter - pCounter)] past_reg = auto_reg_rho[pCounter - 1] - random = np.random.normal(0, 1) - temp.append(past_trs * past_reg + random) + temp.append(past_trs * past_reg) - noise_autoregression.append(np.mean(temp)) + random = np.random.normal(0, 1) + noise_autoregression.append(np.sum(temp) + random) - # Normalize - noise_autoregression = stats.zscore(noise_autoregression) + # N.B. You don't want to normalize. Although that may make the sigma of + # this timecourse 1, it will change the autoregression coefficient to be + # much lower. return noise_autoregression @@ -1172,12 +1515,15 @@ def _generate_noise_temporal_phys(timepoints, heart_freq=1.17, ): """Generate the physiological noise. + Create noise representing the heart rate and respiration of the data. + Default values based on Walvaert, Durnez, Moerkerke, Verdoolaege and + Rosseel, 2011 Parameters ---------- timepoints : 1 Dimensional array - What time points are sampled by a TR + What time points, in seconds, are sampled by a TR resp_freq : float What is the frequency of respiration @@ -1187,19 +1533,23 @@ def _generate_noise_temporal_phys(timepoints, Returns ---------- - one dimensional array, float + noise_phys : one dimensional array, float Generates the physiological temporal noise timecourse + """ noise_phys = [] # Preset + resp_phase = (np.random.rand(1) * 2 * np.pi)[0] + heart_phase = (np.random.rand(1) * 2 * np.pi)[0] for tr_counter in timepoints: - # Calculate the radians for each variable at this given TR - resp_radians = resp_freq * tr_counter * 2 * np.pi - heart_radians = heart_freq * tr_counter * 2 * np.pi + + # Calculate the radians for each variable at this + # given TR + resp_radians = resp_freq * tr_counter * 2 * np.pi + resp_phase + heart_radians = heart_freq * tr_counter * 2 * np.pi + heart_phase # Combine the two types of noise and append - noise_phys.append(np.cos(resp_radians) + np.sin(heart_radians) + - np.random.normal(0, 1)) + noise_phys.append(np.cos(resp_radians) + np.sin(heart_radians)) # Normalize noise_phys = stats.zscore(noise_phys) @@ -1208,6 +1558,7 @@ def _generate_noise_temporal_phys(timepoints, def _generate_noise_spatial(dimensions, + template=None, mask=None, fwhm=4.0, ): @@ -1223,11 +1574,16 @@ def _generate_noise_spatial(dimensions, Parameters ---------- + dimensions : 3 length array, int What is the shape of the volume to be generated - mask : 3d array - The mask describing the boundaries of the brain + template : 3d array, float + A continuous (0 -> 1) volume describing the likelihood a voxel is in + the brain. This can be used to contrast the brain and non brain. + + mask : 3 dimensional array, binary + The masked brain, thresholded to distinguish brain and non-brain fwhm : float What is the full width half max of the gaussian fields being created. @@ -1242,15 +1598,16 @@ def _generate_noise_spatial(dimensions, this will manifest differently on every draw. Secondly, because of the masking and dimensions of the generated volume, this does not behave simply- wrapping effects matter (the outputs are - closer to this value if you have no mask). + closer to the input value if you have no mask). Use _calc_fwhm on this volume alone if you have concerns about the accuracy of the fwhm. Returns ---------- - 3d array, float + noise_spatial : 3d array, float Generates the spatial noise volume for these parameters + """ if len(dimensions) == 4: @@ -1261,6 +1618,7 @@ def logfunc(x, a, b, c): Parameters ---------- + x : float x value of log function @@ -1315,10 +1673,10 @@ def Pk2(idxs): noise_spatial = np.fft.ifftn(noise * amplitude) # Mask or not, then z score - if mask is not None: + if mask is not None and template is not None: # Mask the output - noise_spatial = noise_spatial.real * mask + noise_spatial = noise_spatial.real * template # Z score the specific to the brain noise_spatial[mask > 0] = stats.zscore(noise_spatial[mask > 0]) @@ -1331,6 +1689,7 @@ def Pk2(idxs): def _generate_noise_temporal(stimfunction_tr, tr_duration, dimensions, + template, mask, fwhm, motion_sigma, @@ -1339,10 +1698,9 @@ def _generate_noise_temporal(stimfunction_tr, physiological_sigma, ): """Generate the temporal noise - - To increase or decrease the amount of total noise change the - temporal_noise noise_dict entry. To change the relative mixing of the - noise components, change the sigma's specified below. + Generate the time course of the average brain voxel. To change the + relative mixing of the noise components, change the sigma's specified + below. Parameters ---------- @@ -1354,15 +1712,27 @@ def _generate_noise_temporal(stimfunction_tr, tr_duration : int How long is a TR, in seconds - motion_sigma : float + dimensions : 3 length array, int + What is the shape of the volume to be generated + + template : 3d array, float + A continuous (0 -> 1) volume describing the likelihood a voxel is in + the brain. This can be used to contrast the brain and non brain. + + mask : 3 dimensional array, binary + The masked brain, thresholded to distinguish brain and non-brain + + fwhm : float + What is the full width half max of the gaussian fields being created + to model spatial noise. - How much noise is left over after pre-processing has been - done. This is noise specifically on the task events + motion_sigma : float + This is noise that only occurs for the task events, potentially + representing something like noise due to motion drift_sigma : float - What is the sigma on the distribution that coefficients are - randomly sampled from + What is the sigma on the size of the sine wave auto_reg_sigma : float, list How large is the sigma on the autocorrelation. Higher means more @@ -1376,18 +1746,18 @@ def _generate_noise_temporal(stimfunction_tr, Returns ---------- - one dimensional array, float - Generates the temporal noise timecourse for these parameters + noise_temporal : one dimensional array, float + Generates the temporal noise timecourse for these parameters - """ + """ # Set up common parameters # How many TRs are there trs = len(stimfunction_tr) # What time points are sampled by a TR? - timepoints = list(range(0, trs * tr_duration))[::tr_duration] + timepoints = list(np.linspace(0, (trs - 1) * tr_duration, trs)) noise_drift = _generate_noise_temporal_drift(trs, tr_duration, @@ -1400,18 +1770,19 @@ def _generate_noise_temporal(stimfunction_tr, ) # Generate the volumes that will differ depending on the type of noise - # that it will be used for + # that it will be used for. For drift you want the volume to not have + # the shape of the brain, for the other types of noise you want them to + # have brain shapes volume_drift = np.ones(dimensions) - # volume_drift = _generate_noise_spatial(dimensions=dimensions, - # fwhm=fwhm, - # ) volume_phys = _generate_noise_spatial(dimensions=dimensions, + template=template, mask=mask, fwhm=fwhm, ) volume_autoreg = _generate_noise_spatial(dimensions=dimensions, + template=template, mask=mask, fwhm=fwhm, ) @@ -1435,6 +1806,7 @@ def _generate_noise_temporal(stimfunction_tr, noise_task = _generate_noise_temporal_task(stimfunction_tr, ) volume_task = _generate_noise_spatial(dimensions=dimensions, + template=template, mask=mask, fwhm=fwhm, ) @@ -1444,33 +1816,41 @@ def _generate_noise_temporal(stimfunction_tr, # Finally, z score each voxel so things mix nicely noise_temporal = stats.zscore(noise_temporal, 3) + # If it is a nan it is because you just divided by zero (since some + # voxels are zeros in the template) + noise_temporal[np.isnan(noise_temporal)] = 0 + return noise_temporal def mask_brain(volume, - mask_name=None, - mask_threshold=0.1, + template_name=None, + mask_threshold=None, mask_self=0, ): """ Mask the simulated volume - This takes in a volume and will output the masked volume. if a one - dimensional array is supplied then the output will be a volume of the - dimensions specified in the array. The mask can be created from the - volume by averaging it. All masks will be bounded to the range of 0 to 1. + This creates a mask specifying the likelihood (kind of) a voxel is + part of the brain. All values are bounded to the range of 0 to 1. An + appropriate threshold to isolate brain voxels is >0.2 Parameters ---------- volume : multidimensional array - Either numpy array of the volume that has been simulated or a tuple - describing the dimensions of the mask to be created + Either numpy array of a volume or a tuple describing the dimensions + of the mask to be created - mask_name : str - What is the path to the mask to be loaded? If empty then it defaults - to an MNI152 grey matter mask. + template_name : str + What is the path to the template to be loaded? If empty then it + defaults to an MNI152 grey matter mask. This is ignored if mask_self is + True. mask_threshold : float - What is the threshold (0 -> 1) for including a voxel in the mask? + What is the threshold (0 -> 1) for including a voxel in the mask? If + None then the program will try and identify the last wide peak in a + histogram of the template (assumed to be the brain voxels) and takes + the minima before that peak as the threshold. Won't work when the + data is not bimodal. mask_self : bool If set to true then it makes a mask from the volume supplied (by @@ -1478,8 +1858,14 @@ def mask_brain(volume, Returns ---------- - mask : multidimensional array, float - The masked brain + + mask : 3 dimensional array, binary + The masked brain, thresholded to distinguish brain and non-brain + + template : 3 dimensional array, float + A continuous (0 -> 1) volume describing the likelihood a voxel is in + the brain. This can be used to contrast the brain and non brain. + """ # If the volume supplied is a 1d array then output a volume of the @@ -1488,24 +1874,23 @@ def mask_brain(volume, volume = np.ones(volume) # Load in the mask - if mask_name is None: + if mask_self is True: + mask_raw = volume + elif template_name is None: mask_raw = np.load(resource_stream(__name__, "grey_matter_mask.npy")) else: - mask_raw = np.load(mask_name) + mask_raw = np.load(template_name) - # Is the mask based on the volume - if mask_self is True: - mask_raw = np.zeros([volume.shape[0], volume.shape[1], volume.shape[ - 2], 1]) - - if len(volume.shape) == 4: - mask_raw[:, :, :, 0] = np.mean(volume, 3) - else: - mask_raw[:, :, :, 0] = np.array(volume) - - mask_max = volume.max() + # Make the masks 3d + if len(mask_raw.shape) == 3: + mask_raw = np.array(mask_raw) + elif len(mask_raw.shape) == 4 and mask_raw.shape[3] == 1: + mask_raw = np.array(mask_raw[:, :, :, 0]) else: - mask_max = 1 + mask_raw = np.mean(mask_raw, 3) + + # Find the max value (so you can calulate these as proportions) + mask_max = mask_raw.max() # Make sure the mask values range from 0 to 1 (make out of max of volume # so that this is invertible later) @@ -1524,48 +1909,79 @@ def mask_brain(volume, zoom_factor = (brain_dim[0] / mask_dim[0], brain_dim[1] / mask_dim[1], brain_dim[2] / mask_dim[2], - 1) + ) # Scale the mask according to the input brain - # You might get a warning but ignore it - mask = ndimage.zoom(mask_raw, zoom_factor, order=2) + # You might get a warning if the zoom_factor is not an integer but you + # can safely ignore that. + template = ndimage.zoom(mask_raw, zoom_factor, order=2) + template[template < 0] = 0 - # Any proportion that is below threshold (presumably the exterior of the - # brain), make zero - mask[mask < mask_threshold] = 0 + # If the mask threshold is not supplied then guess it is a minima + # between the two peaks of the bimodal distribution of voxel activity + if mask_threshold is None: - # create the mask in 4d - mask = np.ones(volume.shape) * mask + # How many bins on either side of a peak will be compared + order = 5 - return mask + # Make the histogram + template_vector = template.reshape(brain_dim[0] * brain_dim[1] * + brain_dim[2]) + template_hist = np.histogram(template_vector, 100) + + # Zero pad the values + binval = np.concatenate([np.zeros((order,)), template_hist[0]]) + bins = np.concatenate([np.zeros((order,)), template_hist[1]]) + + # Identify the first two peaks + peaks = signal.argrelmax(binval, order=order)[0][0:2] + + # What is the minima between peaks + minima = binval[peaks[0]:peaks[1]].min() + + # What is the index of the last idx with this min value (since if + # zero, there may be many) + minima_idx = (np.where(binval[peaks[0]:peaks[1]] == minima) + peaks[ + 0])[-1] + + # Convert the minima into a threshold + mask_threshold = bins[minima_idx][0] + + # Mask the template based on the threshold + mask = np.zeros(template.shape) + mask[template > mask_threshold] = 1 + + return mask, template def _noise_dict_update(noise_dict): """ - Update the noise dictionary parameters, in case any were missing + Update the noise dictionary parameters with default values, in case any + were missing Parameters ---------- - noise_dict : dict - A dictionary specifying the types of noise in this experiment. The noise - types interact in important ways. First, all noise types ending with - sigma (e.g. motion sigma) are mixed together in - _generate_temporal_noise. These values describe the proportion of mixing - of these elements. However critically, temporal_noise is the parameter - that describes how much noise these components contribute to the brain. + noise_dict : dict + A dictionary specifying the types of noise in this experiment. The + noise types interact in important ways. First, all noise types + ending with sigma (e.g. motion sigma) are mixed together in + _generate_temporal_noise. These values describe the proportion of + mixing of these elements. However critically, SFNR is the + parameter that describes how much noise these components contribute + to the brain. Returns ------- - Updated dictionary + + noise_dict : dict + Updated dictionary """ # Check what noise is in the dictionary and add if necessary. Numbers # determine relative proportion of noise - if 'temporal_noise' not in noise_dict: - noise_dict['temporal_noise'] = 5 if 'motion_sigma' not in noise_dict: noise_dict['motion_sigma'] = 0 if 'drift_sigma' not in noise_dict: @@ -1576,6 +1992,8 @@ def _noise_dict_update(noise_dict): noise_dict['physiological_sigma'] = 0.1 if 'sfnr' not in noise_dict: noise_dict['sfnr'] = 30 + if 'snr' not in noise_dict: + noise_dict['snr'] = 30 if 'max_activity' not in noise_dict: noise_dict['max_activity'] = 1000 if 'voxel_size' not in noise_dict: @@ -1589,6 +2007,7 @@ def _noise_dict_update(noise_dict): def generate_noise(dimensions, stimfunction_tr, tr_duration, + template, mask=None, noise_dict=None, ): @@ -1596,10 +2015,11 @@ def generate_noise(dimensions, Default noise parameters will create a noise volume with a standard deviation of 0.1 (where the signal defaults to a value of 1). This has built into estimates of how different types of noise mix. All noise - values can be set by the user + values can be set by the user or estimated with calc_noise. Parameters ---------- + dimensions : nd array What is the shape of the volume to be generated @@ -1609,18 +2029,22 @@ def generate_noise(dimensions, tr_duration : float What is the duration, in seconds, of each TR? - mask : 4d array, float - The mask of the brain volume, using + template : 3d array, float + A continuous (0 -> 1) volume describing the likelihood a voxel is in + the brain. This can be used to contrast the brain and non brain. + + mask : 3d array, binary + The mask of the brain volume, distinguishing brain from non-brain noise_dict : dictionary, float This is a dictionary which describes the noise parameters of the - data. If there are no other variables provided then it will default + data. If there are no other variables provided then it will use default values Returns ---------- - multidimensional array, float + noise : multidimensional array, float Generates the noise volume for these parameters """ @@ -1640,13 +2064,14 @@ def generate_noise(dimensions, # Get the mask of the brain and set it to be 3d if mask is None: - mask = np.ones(dimensions_tr) + mask = np.ones(dimensions) # Generate the noise noise_temporal = _generate_noise_temporal(stimfunction_tr=stimfunction_tr, tr_duration=tr_duration, dimensions=dimensions, - mask=mask[:, :, :, 0], + template=template, + mask=mask, fwhm=noise_dict[ 'fwhm'], motion_sigma=noise_dict[ @@ -1659,40 +2084,37 @@ def generate_noise(dimensions, 'physiological_sigma'], ) - # Create the base (this inverts the process to make the mask) - base = mask * noise_dict['max_activity'] - - # Set the amount of background based on the SFNR value + # Create the base (this inverts the process to make the template) + base = template * noise_dict['max_activity'] - # What are the midpoints to be extracted - mid_x_idx = int(np.ceil(base.shape[0] / 2)) - mid_tr_idx = int(np.ceil(base.shape[3] / 2)) + # What is the mean signal of the non masked voxels in this template? + mean_signal = (base[mask > 0]).mean() - # Pull out the slices - slice_volume = base[mid_x_idx, :, :, mid_tr_idx] - slice_mask = mask[mid_x_idx, :, :, mid_tr_idx] + # Convert SFNR into the size of the standard deviation of temporal + # variability + temporal_sd = (mean_signal / noise_dict['sfnr']) - # What is the mean signal of the non masked voxels in this slice? - mean_signal = (slice_volume[slice_mask > 0]).mean() - - # Set up the machine noise - noise_system = _generate_noise_system(dimensions_tr=dimensions_tr) + # Calculate the sd that is necessary to be combined with itself in order + # to generate the temporal_sd + temporal_sd_element = np.sqrt(temporal_sd ** 2 / 2) # What is the standard deviation of the background activity - # (N.B. You need to subtract the mean system noise from this number - # since system_sigma * noise_system.mean() will later be added to the base) - system_sigma = mean_signal / (noise_dict['sfnr'] - noise_system.mean()) + spatial_sd = mean_signal / noise_dict['snr'] - # Increase the size of the system noise based on the SFNR - noise_system *= system_sigma + # Set up the machine noise + noise_system = _generate_noise_system(dimensions_tr=dimensions_tr, + spatial_sd=spatial_sd, + temporal_sd=temporal_sd_element, + ) - # Convert temporal noise (in percent) to real numbers - abs_change = noise_dict['temporal_noise'] * mean_signal / 100 + # Reshape the base (to be the same size as the volume to be created) + base = base.reshape(dimensions[0], dimensions[1], dimensions[2], 1) + base = np.ones(dimensions_tr) * base # Sum up the noise of the brain - noise = base + (noise_temporal * abs_change) + noise_system + noise = base + (noise_temporal * temporal_sd_element) + noise_system - # Reject negative values (only happens outside of the brain + # Reject negative values (only happens outside of the brain) noise[noise < 0] = 0 return noise @@ -1704,6 +2126,8 @@ def plot_brain(fig, percentile=99, ): """ Display the brain that has been generated with a given threshold + Will display the voxels above the given percentile and then a shadow of + all voxels in the mask Parameters ---------- diff --git a/brainiak/utils/utils.py b/brainiak/utils/utils.py index 7fea61920..78c030d7b 100644 --- a/brainiak/utils/utils.py +++ b/brainiak/utils/utils.py @@ -15,7 +15,7 @@ import re import warnings import os.path -from .fmrisim import generate_stimfunction, double_gamma_hrf +from .fmrisim import generate_stimfunction, _double_gamma_hrf, convolve_hrf """ Some utility functions that can be used by different algorithms @@ -432,8 +432,8 @@ def gen_design(stimtime_files, scan_duration, TR, style='FSL', response_delay = hrf_para['response_delay'] undershoot_delay = hrf_para['undershoot_delay'] - response_dispersion = hrf_para['response_dispersion'] - undershoot_dispersion = hrf_para['undershoot_dispersion'] + response_disp = hrf_para['response_dispersion'] + undershoot_disp = hrf_para['undershoot_dispersion'] undershoot_scale = hrf_para['undershoot_scale'] # generate design matrix for i_s in range(n_S): @@ -444,13 +444,15 @@ def gen_design(stimtime_files, scan_duration, TR, style='FSL', total_time=scan_duration[i_s], weights=design_info[i_s][i_c]['weight'], temporal_resolution=1.0/temp_res) - design[i_s][:, i_c] = double_gamma_hrf( - stimfunction, TR, response_delay=response_delay, - undershoot_delay=undershoot_delay, - response_dispersion=response_dispersion, - undershoot_dispersion=undershoot_dispersion, - undershoot_scale=undershoot_scale, scale_function=0, - temporal_resolution=1.0/temp_res) * temp_res + hrf = _double_gamma_hrf(response_delay=response_delay, + undershoot_delay=undershoot_delay, + response_dispersion=response_disp, + undershoot_dispersion=undershoot_disp, + undershoot_scale=undershoot_scale, + temporal_resolution=1.0/temp_res) + design[i_s][:, i_c] = convolve_hrf( + stimfunction, TR, hrf_type=hrf, scale_function=0, + temporal_resolution=1.0 / temp_res).transpose() * temp_res # We multiply the resulting design matrix with # the temporal resolution to normalize it. # We do not use the internal normalization diff --git a/examples/fcma/generate_fcma_data.py b/examples/fcma/generate_fcma_data.py index 996478957..93e74e1b5 100644 --- a/examples/fcma/generate_fcma_data.py +++ b/examples/fcma/generate_fcma_data.py @@ -104,7 +104,7 @@ for cond in list(range(conds)): # Generate a volume representing the location and quality of the signal - volume_static = sim.generate_signal(dimensions=dimensions, + volume_signal = sim.generate_signal(dimensions=dimensions, feature_coordinates=coordinates[cond], feature_type=feature_type, feature_size=feature_size, @@ -126,7 +126,7 @@ # Multiply the HRF timecourse with the signal signal_cond = sim.apply_signal(signal_function=signal_function, - volume_static=volume_static, + volume_signal=volume_signal, ) # Concatenate all the signal and function files @@ -138,11 +138,14 @@ signal += signal_cond # Generate the mask of the signal -mask = sim.mask_brain(signal) +mask, template = sim.mask_brain(signal) # Mask the signal to the shape of a brain (does not attenuate signal according # to grey matter likelihood) -signal *= mask > 0 +signal *= mask.reshape(dimensions[0], dimensions[1], dimensions[2], 1) + +# Downsample the stimulus function to generate it in TR time +stimfunction_tr = stimfunction[::int(tr_duration * 1000)] # Iterate through the participants and store participants epochs = [] @@ -156,8 +159,9 @@ # Create the noise volumes (using the default parameters noise = sim.generate_noise(dimensions=dimensions, - stimfunction=stimfunction, + stimfunction_tr=stimfunction_tr, tr_duration=tr_duration, + template=template, mask=mask, ) @@ -173,7 +177,5 @@ np.save(directory + 'epoch_labels.npy', epochs) # Store the mask -mask[mask > 0] = 1 -mask = mask[:, :, :, 0] brain_nifti = nibabel.Nifti1Image(mask, affine_matrix) nibabel.save(brain_nifti, directory + 'mask.nii') diff --git a/examples/utils/fmrisim_example.py b/examples/utils/fmrisim_example.py index f436228c9..1c4bb194c 100644 --- a/examples/utils/fmrisim_example.py +++ b/examples/utils/fmrisim_example.py @@ -50,14 +50,14 @@ savename = 'examples/utils/example.nii' # Generate a volume representing the location and quality of the signal -volume_static_A = sim.generate_signal(dimensions=dimensions, +volume_signal_A = sim.generate_signal(dimensions=dimensions, feature_coordinates=coordinates_A, feature_type=feature_type, feature_size=feature_size, signal_magnitude=signal_magnitude, ) -volume_static_B = sim.generate_signal(dimensions=dimensions, +volume_signal_B = sim.generate_signal(dimensions=dimensions, feature_coordinates=coordinates_B, feature_type=feature_type, feature_size=feature_size, @@ -67,7 +67,7 @@ # Visualize the signal that was generated for condition A fig = plt.figure() sim.plot_brain(fig, - volume_static_A) + volume_signal_A) plt.show() # Create the time course for the signal to be generated @@ -84,23 +84,23 @@ ) # Convolve the HRF with the stimulus sequence -signal_function_A = sim.double_gamma_hrf(stimfunction=stimfunction_A, - tr_duration=tr_duration, - temporal_resolution=temporal_res, - ) +signal_function_A = sim.convolve_hrf(stimfunction=stimfunction_A, + tr_duration=tr_duration, + temporal_resolution=temporal_res, + ) -signal_function_B = sim.double_gamma_hrf(stimfunction=stimfunction_B, - tr_duration=tr_duration, - temporal_resolution=temporal_res, - ) +signal_function_B = sim.convolve_hrf(stimfunction=stimfunction_B, + tr_duration=tr_duration, + temporal_resolution=temporal_res, + ) # Multiply the HRF timecourse with the signal signal_A = sim.apply_signal(signal_function=signal_function_A, - volume_static=volume_static_A, + volume_signal=volume_signal_A, ) signal_B = sim.apply_signal(signal_function=signal_function_B, - volume_static=volume_static_B, + volume_signal=volume_signal_B, ) # Combine the signals from the two conditions @@ -111,11 +111,11 @@ stimfunction_tr = stimfunction[::int(tr_duration * temporal_res)] # Generate the mask of the signal -mask = sim.mask_brain(signal) +mask, template = sim.mask_brain(signal, mask_threshold=0.2) # Mask the signal to the shape of a brain (attenuates signal according to grey # matter likelihood) -signal *= mask +signal *= mask.reshape(dimensions[0], dimensions[1], dimensions[2], 1) # Generate original noise dict for comparison later orig_noise_dict = sim._noise_dict_update({}) @@ -125,6 +125,7 @@ stimfunction_tr=stimfunction_tr, tr_duration=tr_duration, mask=mask, + template=template, noise_dict=orig_noise_dict, ) @@ -142,7 +143,7 @@ # Get the axis to be plotted ax = sim.plot_brain(fig, brain[:, :, :, tr_counter], - mask=mask[:, :, :, tr_counter], + mask=mask, percentile=99.9) # Wait for an input @@ -167,9 +168,9 @@ stimfunction_tr = stimfunction[::int(tr_duration * temporal_res)] # Calculate the mask -mask = sim.mask_brain(volume=volume, - mask_self=True, - ) +mask, template = sim.mask_brain(volume=volume, + mask_self=True, + ) # Calculate the noise parameters noise_dict = sim.calc_noise(volume=volume, @@ -180,6 +181,7 @@ noise = sim.generate_noise(dimensions=dimensions, tr_duration=tr_duration, stimfunction_tr=stimfunction_tr, + template=template, mask=mask, noise_dict=noise_dict, ) diff --git a/examples/utils/fmrisim_multivariate_example.ipynb b/examples/utils/fmrisim_multivariate_example.ipynb new file mode 100644 index 000000000..02bece51e --- /dev/null +++ b/examples/utils/fmrisim_multivariate_example.ipynb @@ -0,0 +1,5525 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Copyright 2016 Intel Corporation\n", + "\n", + "Licensed under the Apache License, Version 2.0 (the \"License\");\n", + "you may not use this file except in compliance with the License.\n", + "You may obtain a copy of the License at\n", + "http://www.apache.org/licenses/LICENSE-2.0\n", + "Unless required by applicable law or agreed to in writing, software\n", + "distributed under the License is distributed on an \"AS IS\" BASIS,\n", + "WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n", + "See the License for the specific language governing permissions and\n", + "limitations under the License.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# fMRI Simulator example script for multivariate analyses\n", + "\n", + "Example script to demonstrate fmrisim functionality. This generates\n", + "data for a two condition, event related design in which each condition\n", + "evokes different activity within the same voxels. It then runs simple \n", + "univariate and multivariate analyses on the data\n", + "\n", + "Authors: Cameron Ellis (Yale) 2017\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### **1.\tSet parameters.** \n", + "\n", + "It is necessary to set various parameters that describe how the signal and the noise will be generated." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*1.1 Import necessary Python packages*" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "%matplotlib notebook\n", + "\n", + "from pathlib import Path\n", + "from brainiak.utils import fmrisim\n", + "import nibabel\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import scipy.spatial.distance as sp_distance\n", + "import sklearn.manifold as manifold\n", + "import scipy.stats as stats\n", + "import sklearn.model_selection\n", + "import sklearn.svm" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*1.2 Load participant data*\n", + "\n", + "Any 4 dimensional fMRI data that is readible by nibabel can be used as input to this pipeline. For this example, data is taken from the open access repository DataSpace: http://arks.princeton.edu/ark:/88435/dsp01dn39x4181. This file is unzipped and placed in the home directory with the name Corr_MVPA " + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "home = str(Path.home())\n", + "nii = nibabel.load(home + '/Corr_MVPA/Participant_01_rest_run01.nii')\n", + "volume = nii.get_data()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*1.2\tSpecify participant dimensions and resolution*\n", + "\n", + "It is possible to manually specify all parameters necessary for fmrisim. However, it is also possible to instead provide an fMRI dataset as input and extract the necessary parameters from that dataset. Such an example is described below and will be followed throughout. Here the size of the volume and the resolution of the voxels within it are determined." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "dim = volume.shape\n", + "dimsize = nii.header.get_zooms()\n", + "tr = dimsize[3]\n", + "if tr > 100: # If high then these values are likely in ms\n", + " tr /= 1000" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*1.3 Generate an activity template and a mask*\n", + "\n", + "Functions in fmrisim require a continuous map that describes the appropriate average MR value for each voxel in the brain and a mask which specifies voxels in the brain versus voxels outside of the brain. One way to generate both of these volumes is the mask_brain function. At a minimum, this takes as an input the fMRI volume to be simulated. To create the template this volume is averaged over time and bounded to a range from 0 to 1. In other words, voxels with a high value in the template have high activity over time. To create a mask, the template is thresholded. This threshold can be set manually or instead an appropriate value can be determined by looking for the minima between the two first peaks in the histogram of voxel values.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "mask, template = fmrisim.mask_brain(volume=volume, \n", + " mask_self=True,\n", + " )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*1.4 Determine noise parameters*\n", + "\n", + "A critical step in the fmrisim toolbox is determining the noise parameters of the volume to be created. Many noise parameters are available for specification and if any are not set then they will default to reasonable values. As before, it is instead possible to provide raw fMRI data that will be used to estimate these noise parameters. The goal of the noise estimation is to calculate general descriptive statistics about the noise in the brain that are thought to be important. The simulations are thought to be useful for understanding how signals will survive analyses when embedded in realistic neural noise. \n", + "\n", + "Now the disclaimers: the values here are only an estimate and will depend on noise properties combining in the ways specified. In addition, because of the non-linearity and stochasticity of this simulation, this estimation is not fully invertible: if you generate a dataset with a set of noise parameters it will have similar but not the same noise parameters as a result. Moreover, complex interactions between brain regions that likely better describe brain noise are not modelled here: this toolbox pays no attention to regions of the brain or their interactions. Finally, for best results use raw fMRI because if the data has been preprocessed then assumptions this algorithm makes are likely to be erroneous. For instance, if the brain has been masked then this will eliminate variance in non-brain voxels which will mean that calculations of noise dependent on those voxels as a reference will fail.\n", + "\n", + "This toolbox separates noise in two: spatial noise and temporal noise. To estimate spatial noise both the smoothness and the amount of non-brain noise of the data must be quantified. For smoothness, the Full Width Half Max (FWHM) of the volume is averaged for the X, Y and Z dimension and then averaged across a sample of time points. To calculate the Signal to Noise Ratio the mean activity in brain voxels for the middle time point is divided by the standard deviation in activity across non-brain voxels for that time point. For temporal noise the drift, temporal autocorrelation, and functional variability is estimated. The drift is estimated by averaging all non-brain voxels and looking at the variance in this average across time. This time course is also used to estimate the temporal autoregression by taking the first slope coefficient of an autoregression estimation function from the Nitime package . The Signal to Fluctuation Noise Ratio is calculated by dividing the average activity of voxels in the brain with that voxel’s noise (Friedman & Glover, 2006). That noise is calculated by taking the standard deviation of that voxel over time after it has been detrended with a second order polynomial. Other types of noise can be generated, such as physiological noise, but are not estimated by this function.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "noise_dict = {'voxel_size': [dimsize[0], dimsize[1], dimsize[2]]}\n", + "noise_dict = fmrisim.calc_noise(volume=volume,\n", + " mask=mask,\n", + " noise_dict=noise_dict,\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Noise parameters of the data were estimated as follows:\n", + "SNR: 69.5116700015\n", + "SFNR: 70.7171164885\n", + "FWHM: 5.65860178162\n" + ] + } + ], + "source": [ + "print('Noise parameters of the data were estimated as follows:')\n", + "print('SNR: ' + str(noise_dict['snr']))\n", + "print('SFNR: ' + str(noise_dict['sfnr']))\n", + "print('FWHM: ' + str(noise_dict['fwhm']))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### **2. Generate signal**\n", + "\n", + "fmrisim can be used to generate signal in a number of different ways depending on the type of effect being simulated. Several tools are supplied to help with different types of signal that may be required; however, custom scripts may be necessary for unique effects. Below an experiment will be simulated in which two conditions, A and B, evoke different patterns of activity in the same set of voxels in the brain. This pattern does not manifest as a univariate change in voxel activity across voxels but instead each condition evokes a consistent pattern across voxels. These conditions are randomly intermixed trial by trial. This code could be easily changed to instead compare univariate changes evoked by stimuli in different brain regions. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*2.1 Establish effect size*\n", + "\n", + "When specifying the signal we must determine the amount of activity change each voxel undergoes. A useful metric for this is the SFNR value determined from noise calculations because it can be used to estimate the variability in the average voxel. For a univariate effect, to estimate activity with a Cohen’s d of 1, the size of the change must be equivalent to one standard deviation. For multivariate effects the effect size depends on multiple factors including the number of voxels and conditions. Different measures for effect size could also be calculated, such as percent signal change." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "effect_size = 1\n", + "temporal_sd = (template[mask > 0]).mean() * noise_dict['max_activity'] / noise_dict['sfnr']\n", + "effect_signal_change = effect_size * temporal_sd" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*2.2 Characterize signal for voxels*\n", + "\n", + "Specify the pattern of activity across a given number of voxels that characterizes each condition. This pattern can simply be random, as is done here, or can be structured, like the position of voxels in high dimensional representation space." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "feature_size = 3\n", + "voxels = feature_size ** 3\n", + "pattern_A = np.random.randn(voxels).reshape((voxels, 1)) * effect_signal_change\n", + "pattern_B = np.random.randn(voxels).reshape((voxels, 1)) * effect_signal_change" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "window.mpl = {};\n", + "\n", + "\n", + "mpl.get_websocket_type = function() {\n", + " if (typeof(WebSocket) !== 'undefined') {\n", + " return WebSocket;\n", + " } else if (typeof(MozWebSocket) !== 'undefined') {\n", + " return MozWebSocket;\n", + " } else {\n", + " alert('Your browser does not have WebSocket support.' +\n", + " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", + " 'Firefox 4 and 5 are also supported but you ' +\n", + " 'have to enable WebSockets in about:config.');\n", + " };\n", + "}\n", + "\n", + "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", + " this.id = figure_id;\n", + "\n", + " this.ws = websocket;\n", + "\n", + " this.supports_binary = (this.ws.binaryType != undefined);\n", + "\n", + " if (!this.supports_binary) {\n", + " var warnings = document.getElementById(\"mpl-warnings\");\n", + " if (warnings) {\n", + " warnings.style.display = 'block';\n", + " warnings.textContent = (\n", + " \"This browser does not support binary websocket messages. \" +\n", + " \"Performance may be slow.\");\n", + " }\n", + " }\n", + "\n", + " this.imageObj = new Image();\n", + "\n", + " this.context = undefined;\n", + " this.message = undefined;\n", + " this.canvas = undefined;\n", + " this.rubberband_canvas = undefined;\n", + " this.rubberband_context = undefined;\n", + " this.format_dropdown = undefined;\n", + "\n", + " this.image_mode = 'full';\n", + "\n", + " this.root = $('
');\n", + " this._root_extra_style(this.root)\n", + " this.root.attr('style', 'display: inline-block');\n", + "\n", + " $(parent_element).append(this.root);\n", + "\n", + " this._init_header(this);\n", + " this._init_canvas(this);\n", + " this._init_toolbar(this);\n", + "\n", + " var fig = this;\n", + "\n", + " this.waiting = false;\n", + "\n", + " this.ws.onopen = function () {\n", + " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", + " fig.send_message(\"send_image_mode\", {});\n", + " if (mpl.ratio != 1) {\n", + " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", + " }\n", + " fig.send_message(\"refresh\", {});\n", + " }\n", + "\n", + " this.imageObj.onload = function() {\n", + " if (fig.image_mode == 'full') {\n", + " // Full images could contain transparency (where diff images\n", + " // almost always do), so we need to clear the canvas so that\n", + " // there is no ghosting.\n", + " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", + " }\n", + " fig.context.drawImage(fig.imageObj, 0, 0);\n", + " };\n", + "\n", + " this.imageObj.onunload = function() {\n", + " this.ws.close();\n", + " }\n", + "\n", + " this.ws.onmessage = this._make_on_message_function(this);\n", + "\n", + " this.ondownload = ondownload;\n", + "}\n", + "\n", + "mpl.figure.prototype._init_header = function() {\n", + " var titlebar = $(\n", + " '
');\n", + " var titletext = $(\n", + " '
');\n", + " titlebar.append(titletext)\n", + " this.root.append(titlebar);\n", + " this.header = titletext[0];\n", + "}\n", + "\n", + "\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._init_canvas = function() {\n", + " var fig = this;\n", + "\n", + " var canvas_div = $('
');\n", + "\n", + " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", + "\n", + " function canvas_keyboard_event(event) {\n", + " return fig.key_event(event, event['data']);\n", + " }\n", + "\n", + " canvas_div.keydown('key_press', canvas_keyboard_event);\n", + " canvas_div.keyup('key_release', canvas_keyboard_event);\n", + " this.canvas_div = canvas_div\n", + " this._canvas_extra_style(canvas_div)\n", + " this.root.append(canvas_div);\n", + "\n", + " var canvas = $('');\n", + " canvas.addClass('mpl-canvas');\n", + " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", + "\n", + " this.canvas = canvas[0];\n", + " this.context = canvas[0].getContext(\"2d\");\n", + "\n", + " var backingStore = this.context.backingStorePixelRatio ||\n", + "\tthis.context.webkitBackingStorePixelRatio ||\n", + "\tthis.context.mozBackingStorePixelRatio ||\n", + "\tthis.context.msBackingStorePixelRatio ||\n", + "\tthis.context.oBackingStorePixelRatio ||\n", + "\tthis.context.backingStorePixelRatio || 1;\n", + "\n", + " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", + "\n", + " var rubberband = $('');\n", + " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", + "\n", + " var pass_mouse_events = true;\n", + "\n", + " canvas_div.resizable({\n", + " start: function(event, ui) {\n", + " pass_mouse_events = false;\n", + " },\n", + " resize: function(event, ui) {\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " stop: function(event, ui) {\n", + " pass_mouse_events = true;\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " });\n", + "\n", + " function mouse_event_fn(event) {\n", + " if (pass_mouse_events)\n", + " return fig.mouse_event(event, event['data']);\n", + " }\n", + "\n", + " rubberband.mousedown('button_press', mouse_event_fn);\n", + " rubberband.mouseup('button_release', mouse_event_fn);\n", + " // Throttle sequential mouse events to 1 every 20ms.\n", + " rubberband.mousemove('motion_notify', mouse_event_fn);\n", + "\n", + " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", + " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", + "\n", + " canvas_div.on(\"wheel\", function (event) {\n", + " event = event.originalEvent;\n", + " event['data'] = 'scroll'\n", + " if (event.deltaY < 0) {\n", + " event.step = 1;\n", + " } else {\n", + " event.step = -1;\n", + " }\n", + " mouse_event_fn(event);\n", + " });\n", + "\n", + " canvas_div.append(canvas);\n", + " canvas_div.append(rubberband);\n", + "\n", + " this.rubberband = rubberband;\n", + " this.rubberband_canvas = rubberband[0];\n", + " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", + " this.rubberband_context.strokeStyle = \"#000000\";\n", + "\n", + " this._resize_canvas = function(width, height) {\n", + " // Keep the size of the canvas, canvas container, and rubber band\n", + " // canvas in synch.\n", + " canvas_div.css('width', width)\n", + " canvas_div.css('height', height)\n", + "\n", + " canvas.attr('width', width * mpl.ratio);\n", + " canvas.attr('height', height * mpl.ratio);\n", + " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", + "\n", + " rubberband.attr('width', width);\n", + " rubberband.attr('height', height);\n", + " }\n", + "\n", + " // Set the figure to an initial 600x600px, this will subsequently be updated\n", + " // upon first draw.\n", + " this._resize_canvas(600, 600);\n", + "\n", + " // Disable right mouse context menu.\n", + " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", + " return false;\n", + " });\n", + "\n", + " function set_focus () {\n", + " canvas.focus();\n", + " canvas_div.focus();\n", + " }\n", + "\n", + " window.setTimeout(set_focus, 100);\n", + "}\n", + "\n", + "mpl.figure.prototype._init_toolbar = function() {\n", + " var fig = this;\n", + "\n", + " var nav_element = $('
')\n", + " nav_element.attr('style', 'width: 100%');\n", + " this.root.append(nav_element);\n", + "\n", + " // Define a callback function for later on.\n", + " function toolbar_event(event) {\n", + " return fig.toolbar_button_onclick(event['data']);\n", + " }\n", + " function toolbar_mouse_event(event) {\n", + " return fig.toolbar_button_onmouseover(event['data']);\n", + " }\n", + "\n", + " for(var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " // put a spacer in here.\n", + " continue;\n", + " }\n", + " var button = $('');\n", + " button.click(method_name, toolbar_event);\n", + " button.mouseover(tooltip, toolbar_mouse_event);\n", + " nav_element.append(button);\n", + " }\n", + "\n", + " // Add the status bar.\n", + " var status_bar = $('');\n", + " nav_element.append(status_bar);\n", + " this.message = status_bar[0];\n", + "\n", + " // Add the close button to the window.\n", + " var buttongrp = $('
');\n", + " var button = $('');\n", + " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", + " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", + " buttongrp.append(button);\n", + " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", + " titlebar.prepend(buttongrp);\n", + "}\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(el){\n", + " var fig = this\n", + " el.on(\"remove\", function(){\n", + "\tfig.close_ws(fig, {});\n", + " });\n", + "}\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(el){\n", + " // this is important to make the div 'focusable\n", + " el.attr('tabindex', 0)\n", + " // reach out to IPython and tell the keyboard manager to turn it's self\n", + " // off when our div gets focus\n", + "\n", + " // location in version 3\n", + " if (IPython.notebook.keyboard_manager) {\n", + " IPython.notebook.keyboard_manager.register_events(el);\n", + " }\n", + " else {\n", + " // location in version 2\n", + " IPython.keyboard_manager.register_events(el);\n", + " }\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._key_event_extra = function(event, name) {\n", + " var manager = IPython.notebook.keyboard_manager;\n", + " if (!manager)\n", + " manager = IPython.keyboard_manager;\n", + "\n", + " // Check for shift+enter\n", + " if (event.shiftKey && event.which == 13) {\n", + " this.canvas_div.blur();\n", + " // select the cell after this one\n", + " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", + " IPython.notebook.select(index + 1);\n", + " }\n", + "}\n", + "\n", + "mpl.figure.prototype.handle_save = function(fig, msg) {\n", + " fig.ondownload(fig, null);\n", + "}\n", + "\n", + "\n", + "mpl.find_output_cell = function(html_output) {\n", + " // Return the cell and output element which can be found *uniquely* in the notebook.\n", + " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", + " // IPython event is triggered only after the cells have been serialised, which for\n", + " // our purposes (turning an active figure into a static one), is too late.\n", + " var cells = IPython.notebook.get_cells();\n", + " var ncells = cells.length;\n", + " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", + " data = data.data;\n", + " }\n", + " if (data['text/html'] == html_output) {\n", + " return [cell, data, j];\n", + " }\n", + " }\n", + " }\n", + " }\n", + "}\n", + "\n", + "// Register the function which deals with the matplotlib target/channel.\n", + "// The kernel may be null if the page has been refreshed.\n", + "if (IPython.notebook.kernel != null) {\n", + " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", + "}\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Display signal\n", + "plt.figure()\n", + "\n", + "response = stats.zscore(signal_func[0:100,0])\n", + "plt.title('Example event time course and voxel response')\n", + "Event_A = plt.plot(stimfunc_A[0:100, 0], 'r', label='Event_A')\n", + "Event_B = plt.plot(stimfunc_B[0:100, 0], 'g', label='Event_B')\n", + "Response = plt.plot(response, 'b', label='Response')\n", + "plt.legend(loc=1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*2.7 Specify which voxels in the brain contain signal*\n", + "\n", + "fmrisim provides tools to specify certain voxels in the brain that contain signal. The generate_signal function can produce regions of activity in a brain of different shapes, such as cubes, loops and spheres. Alternatively a volume could be loaded in that specifies the signal voxels (e.g. for ROI analyses). The value of each voxel can be specified here, or set to be a random value." + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "coordinates = np.array([[24, 24, 24]])\n", + "signal_volume = fmrisim.generate_signal(dimensions=dim[0:3],\n", + " feature_type=['cube'],\n", + " feature_coordinates=coordinates,\n", + " feature_size=[feature_size],\n", + " signal_magnitude=[1],\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "window.mpl = {};\n", + "\n", + "\n", + "mpl.get_websocket_type = function() {\n", + " if (typeof(WebSocket) !== 'undefined') {\n", + " return WebSocket;\n", + " } else if (typeof(MozWebSocket) !== 'undefined') {\n", + " return MozWebSocket;\n", + " } else {\n", + " alert('Your browser does not have WebSocket support.' +\n", + " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", + " 'Firefox 4 and 5 are also supported but you ' +\n", + " 'have to enable WebSockets in about:config.');\n", + " };\n", + "}\n", + "\n", + "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", + " this.id = figure_id;\n", + "\n", + " this.ws = websocket;\n", + "\n", + " this.supports_binary = (this.ws.binaryType != undefined);\n", + "\n", + " if (!this.supports_binary) {\n", + " var warnings = document.getElementById(\"mpl-warnings\");\n", + " if (warnings) {\n", + " warnings.style.display = 'block';\n", + " warnings.textContent = (\n", + " \"This browser does not support binary websocket messages. \" +\n", + " \"Performance may be slow.\");\n", + " }\n", + " }\n", + "\n", + " this.imageObj = new Image();\n", + "\n", + " this.context = undefined;\n", + " this.message = undefined;\n", + " this.canvas = undefined;\n", + " this.rubberband_canvas = undefined;\n", + " this.rubberband_context = undefined;\n", + " this.format_dropdown = undefined;\n", + "\n", + " this.image_mode = 'full';\n", + "\n", + " this.root = $('
');\n", + " this._root_extra_style(this.root)\n", + " this.root.attr('style', 'display: inline-block');\n", + "\n", + " $(parent_element).append(this.root);\n", + "\n", + " this._init_header(this);\n", + " this._init_canvas(this);\n", + " this._init_toolbar(this);\n", + "\n", + " var fig = this;\n", + "\n", + " this.waiting = false;\n", + "\n", + " this.ws.onopen = function () {\n", + " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", + " fig.send_message(\"send_image_mode\", {});\n", + " if (mpl.ratio != 1) {\n", + " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", + " }\n", + " fig.send_message(\"refresh\", {});\n", + " }\n", + "\n", + " this.imageObj.onload = function() {\n", + " if (fig.image_mode == 'full') {\n", + " // Full images could contain transparency (where diff images\n", + " // almost always do), so we need to clear the canvas so that\n", + " // there is no ghosting.\n", + " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", + " }\n", + " fig.context.drawImage(fig.imageObj, 0, 0);\n", + " };\n", + "\n", + " this.imageObj.onunload = function() {\n", + " this.ws.close();\n", + " }\n", + "\n", + " this.ws.onmessage = this._make_on_message_function(this);\n", + "\n", + " this.ondownload = ondownload;\n", + "}\n", + "\n", + "mpl.figure.prototype._init_header = function() {\n", + " var titlebar = $(\n", + " '
');\n", + " var titletext = $(\n", + " '
');\n", + " titlebar.append(titletext)\n", + " this.root.append(titlebar);\n", + " this.header = titletext[0];\n", + "}\n", + "\n", + "\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._init_canvas = function() {\n", + " var fig = this;\n", + "\n", + " var canvas_div = $('
');\n", + "\n", + " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", + "\n", + " function canvas_keyboard_event(event) {\n", + " return fig.key_event(event, event['data']);\n", + " }\n", + "\n", + " canvas_div.keydown('key_press', canvas_keyboard_event);\n", + " canvas_div.keyup('key_release', canvas_keyboard_event);\n", + " this.canvas_div = canvas_div\n", + " this._canvas_extra_style(canvas_div)\n", + " this.root.append(canvas_div);\n", + "\n", + " var canvas = $('');\n", + " canvas.addClass('mpl-canvas');\n", + " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", + "\n", + " this.canvas = canvas[0];\n", + " this.context = canvas[0].getContext(\"2d\");\n", + "\n", + " var backingStore = this.context.backingStorePixelRatio ||\n", + "\tthis.context.webkitBackingStorePixelRatio ||\n", + "\tthis.context.mozBackingStorePixelRatio ||\n", + "\tthis.context.msBackingStorePixelRatio ||\n", + "\tthis.context.oBackingStorePixelRatio ||\n", + "\tthis.context.backingStorePixelRatio || 1;\n", + "\n", + " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", + "\n", + " var rubberband = $('');\n", + " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", + "\n", + " var pass_mouse_events = true;\n", + "\n", + " canvas_div.resizable({\n", + " start: function(event, ui) {\n", + " pass_mouse_events = false;\n", + " },\n", + " resize: function(event, ui) {\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " stop: function(event, ui) {\n", + " pass_mouse_events = true;\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " });\n", + "\n", + " function mouse_event_fn(event) {\n", + " if (pass_mouse_events)\n", + " return fig.mouse_event(event, event['data']);\n", + " }\n", + "\n", + " rubberband.mousedown('button_press', mouse_event_fn);\n", + " rubberband.mouseup('button_release', mouse_event_fn);\n", + " // Throttle sequential mouse events to 1 every 20ms.\n", + " rubberband.mousemove('motion_notify', mouse_event_fn);\n", + "\n", + " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", + " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", + "\n", + " canvas_div.on(\"wheel\", function (event) {\n", + " event = event.originalEvent;\n", + " event['data'] = 'scroll'\n", + " if (event.deltaY < 0) {\n", + " event.step = 1;\n", + " } else {\n", + " event.step = -1;\n", + " }\n", + " mouse_event_fn(event);\n", + " });\n", + "\n", + " canvas_div.append(canvas);\n", + " canvas_div.append(rubberband);\n", + "\n", + " this.rubberband = rubberband;\n", + " this.rubberband_canvas = rubberband[0];\n", + " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", + " this.rubberband_context.strokeStyle = \"#000000\";\n", + "\n", + " this._resize_canvas = function(width, height) {\n", + " // Keep the size of the canvas, canvas container, and rubber band\n", + " // canvas in synch.\n", + " canvas_div.css('width', width)\n", + " canvas_div.css('height', height)\n", + "\n", + " canvas.attr('width', width * mpl.ratio);\n", + " canvas.attr('height', height * mpl.ratio);\n", + " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", + "\n", + " rubberband.attr('width', width);\n", + " rubberband.attr('height', height);\n", + " }\n", + "\n", + " // Set the figure to an initial 600x600px, this will subsequently be updated\n", + " // upon first draw.\n", + " this._resize_canvas(600, 600);\n", + "\n", + " // Disable right mouse context menu.\n", + " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", + " return false;\n", + " });\n", + "\n", + " function set_focus () {\n", + " canvas.focus();\n", + " canvas_div.focus();\n", + " }\n", + "\n", + " window.setTimeout(set_focus, 100);\n", + "}\n", + "\n", + "mpl.figure.prototype._init_toolbar = function() {\n", + " var fig = this;\n", + "\n", + " var nav_element = $('
')\n", + " nav_element.attr('style', 'width: 100%');\n", + " this.root.append(nav_element);\n", + "\n", + " // Define a callback function for later on.\n", + " function toolbar_event(event) {\n", + " return fig.toolbar_button_onclick(event['data']);\n", + " }\n", + " function toolbar_mouse_event(event) {\n", + " return fig.toolbar_button_onmouseover(event['data']);\n", + " }\n", + "\n", + " for(var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " // put a spacer in here.\n", + " continue;\n", + " }\n", + " var button = $('');\n", + " button.click(method_name, toolbar_event);\n", + " button.mouseover(tooltip, toolbar_mouse_event);\n", + " nav_element.append(button);\n", + " }\n", + "\n", + " // Add the status bar.\n", + " var status_bar = $('');\n", + " nav_element.append(status_bar);\n", + " this.message = status_bar[0];\n", + "\n", + " // Add the close button to the window.\n", + " var buttongrp = $('
');\n", + " var button = $('');\n", + " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", + " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", + " buttongrp.append(button);\n", + " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", + " titlebar.prepend(buttongrp);\n", + "}\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(el){\n", + " var fig = this\n", + " el.on(\"remove\", function(){\n", + "\tfig.close_ws(fig, {});\n", + " });\n", + "}\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(el){\n", + " // this is important to make the div 'focusable\n", + " el.attr('tabindex', 0)\n", + " // reach out to IPython and tell the keyboard manager to turn it's self\n", + " // off when our div gets focus\n", + "\n", + " // location in version 3\n", + " if (IPython.notebook.keyboard_manager) {\n", + " IPython.notebook.keyboard_manager.register_events(el);\n", + " }\n", + " else {\n", + " // location in version 2\n", + " IPython.keyboard_manager.register_events(el);\n", + " }\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._key_event_extra = function(event, name) {\n", + " var manager = IPython.notebook.keyboard_manager;\n", + " if (!manager)\n", + " manager = IPython.keyboard_manager;\n", + "\n", + " // Check for shift+enter\n", + " if (event.shiftKey && event.which == 13) {\n", + " this.canvas_div.blur();\n", + " // select the cell after this one\n", + " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", + " IPython.notebook.select(index + 1);\n", + " }\n", + "}\n", + "\n", + "mpl.figure.prototype.handle_save = function(fig, msg) {\n", + " fig.ondownload(fig, null);\n", + "}\n", + "\n", + "\n", + "mpl.find_output_cell = function(html_output) {\n", + " // Return the cell and output element which can be found *uniquely* in the notebook.\n", + " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", + " // IPython event is triggered only after the cells have been serialised, which for\n", + " // our purposes (turning an active figure into a static one), is too late.\n", + " var cells = IPython.notebook.get_cells();\n", + " var ncells = cells.length;\n", + " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", + " data = data.data;\n", + " }\n", + " if (data['text/html'] == html_output) {\n", + " return [cell, data, j];\n", + " }\n", + " }\n", + " }\n", + " }\n", + "}\n", + "\n", + "// Register the function which deals with the matplotlib target/channel.\n", + "// The kernel may be null if the page has been refreshed.\n", + "if (IPython.notebook.kernel != null) {\n", + " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", + "}\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "(-0.5, 63.5, 63.5, -0.5)" + ] + }, + "execution_count": 20, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "plt.figure()\n", + "plt.imshow(noise[:, :, 24, 0], cmap=plt.cm.gray)\n", + "plt.axis('off')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*3.1 Create temporal noise*\n", + "\n", + "The temporal noise of fMRI data is comprised of multiple components: drift, autoregression, task related motion and physiological noise. To estimate drift, a sine wave with a default period of 300s, is used. Although other simulators use a combination of discrete cosines for estimating drift (Welvaert, et al., 2011), our testing suggests that this provides poor estimates in long scans (>200s). This drift is then multiplied by a three-dimensional volume of Gaussian random fields of a specific FWHM. Autoregression noise is estimated by creating a time course of Gaussian noise values that are weighted by previous values of the time course. This autoregressive time course is multiplied by a brain shaped volume of Gaussian random fields. Physiological noise is modeled by sine waves comprised of heart rate (1.17Hz) and respiration rate (0.2Hz) (Biswal, et al., 1996) with random phase. This time course is also multiplied by brain shaped spatial noise. Finally, task related noise is simulated by adding Gaussian or Rician noise to time points where there are events (according to the event time course) and in turn this is multiplied by a brain shaped spatial noise volume. These four noise components are then mixed together in proportion to the size of their corresponding noise values. This aggregated volume is then Z scored and the SFNR is used to estimate the appropriate standard deviation of these values across time. \n", + "\t\n", + "*3.2 Create system noise*\n", + " \n", + "In addition to temporal noise from fluctuations in the scanner there is also machine noise that causes fluctuations in all voxels. When SNR is low, Rician noise is a good estimate of background noise data (Gudbjartsson, & Patz, 1995). From our testing, when SNR is higher than 30 then noise with an exponential distribution better describes the data. The SNR value that is supplied determines the standard deviation of this machine noise.\t\n", + "\n", + "*3.3 Combine noise and template*\n", + " \n", + "The template volume is used to estimate the appropriate baseline distribution of MR values. This estimate is then combined with the temporal noise and the system noise to make an estimate of the noise. \n", + "\n", + "*3.4 Combine signal and noise*\n", + "\n", + "Since the brain signal is expected to be small and sparse relative to the noise, it is assumed sufficient to simply add the volume containing signal with the volume modeling noise to make the simulated brain. \n" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "brain = signal + noise" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### **4. Analyse data**\n", + "\n", + "Several tools are available for multivariate analysis in BrainIAK. These greatly speed up computation and are critical in some cases, such as a whole brain searchlight. However, for this example data we will only look at data in the ROI that we know contains signal." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*4.1 Pull out data for each trial*\n", + "\n", + "Identify which voxels are in the signal ROI by using the coordinates provided earlier. To identify the relevant timepoints, assume that the peak of the neural response occurs 4 - 6s after each event onset. Take the TR corresponding to this peak response as the TR for that trial. In longer event/block designs you might instead average over each event." + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "hrf_lag = 4 # Assumed time from stimulus onset to HRF peak\n", + "\n", + "# Get the lower and upper bounds of the ROI\n", + "lb = (coordinates - ((feature_size - 1) / 2)).astype('int')[0]\n", + "ub = (coordinates + ((feature_size - 1) / 2) + 1).astype('int')[0]\n", + "\n", + "trials_A = brain[lb[0]:ub[0], lb[1]:ub[1], lb[2]:ub[2], (onsets_A + hrf_lag / tr).astype('int')]\n", + "trials_B = brain[lb[0]:ub[0], lb[1]:ub[1], lb[2]:ub[2], (onsets_B + hrf_lag / tr).astype('int')]\n", + "\n", + "trials_A = trials_A.reshape((voxels, trials_A.shape[3]))\n", + "trials_B = trials_B.reshape((voxels, trials_B.shape[3]))" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [ + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "window.mpl = {};\n", + "\n", + "\n", + "mpl.get_websocket_type = function() {\n", + " if (typeof(WebSocket) !== 'undefined') {\n", + " return WebSocket;\n", + " } else if (typeof(MozWebSocket) !== 'undefined') {\n", + " return MozWebSocket;\n", + " } else {\n", + " alert('Your browser does not have WebSocket support.' +\n", + " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", + " 'Firefox 4 and 5 are also supported but you ' +\n", + " 'have to enable WebSockets in about:config.');\n", + " };\n", + "}\n", + "\n", + "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", + " this.id = figure_id;\n", + "\n", + " this.ws = websocket;\n", + "\n", + " this.supports_binary = (this.ws.binaryType != undefined);\n", + "\n", + " if (!this.supports_binary) {\n", + " var warnings = document.getElementById(\"mpl-warnings\");\n", + " if (warnings) {\n", + " warnings.style.display = 'block';\n", + " warnings.textContent = (\n", + " \"This browser does not support binary websocket messages. \" +\n", + " \"Performance may be slow.\");\n", + " }\n", + " }\n", + "\n", + " this.imageObj = new Image();\n", + "\n", + " this.context = undefined;\n", + " this.message = undefined;\n", + " this.canvas = undefined;\n", + " this.rubberband_canvas = undefined;\n", + " this.rubberband_context = undefined;\n", + " this.format_dropdown = undefined;\n", + "\n", + " this.image_mode = 'full';\n", + "\n", + " this.root = $('
');\n", + " this._root_extra_style(this.root)\n", + " this.root.attr('style', 'display: inline-block');\n", + "\n", + " $(parent_element).append(this.root);\n", + "\n", + " this._init_header(this);\n", + " this._init_canvas(this);\n", + " this._init_toolbar(this);\n", + "\n", + " var fig = this;\n", + "\n", + " this.waiting = false;\n", + "\n", + " this.ws.onopen = function () {\n", + " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", + " fig.send_message(\"send_image_mode\", {});\n", + " if (mpl.ratio != 1) {\n", + " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", + " }\n", + " fig.send_message(\"refresh\", {});\n", + " }\n", + "\n", + " this.imageObj.onload = function() {\n", + " if (fig.image_mode == 'full') {\n", + " // Full images could contain transparency (where diff images\n", + " // almost always do), so we need to clear the canvas so that\n", + " // there is no ghosting.\n", + " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", + " }\n", + " fig.context.drawImage(fig.imageObj, 0, 0);\n", + " };\n", + "\n", + " this.imageObj.onunload = function() {\n", + " this.ws.close();\n", + " }\n", + "\n", + " this.ws.onmessage = this._make_on_message_function(this);\n", + "\n", + " this.ondownload = ondownload;\n", + "}\n", + "\n", + "mpl.figure.prototype._init_header = function() {\n", + " var titlebar = $(\n", + " '
');\n", + " var titletext = $(\n", + " '
');\n", + " titlebar.append(titletext)\n", + " this.root.append(titlebar);\n", + " this.header = titletext[0];\n", + "}\n", + "\n", + "\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._init_canvas = function() {\n", + " var fig = this;\n", + "\n", + " var canvas_div = $('
');\n", + "\n", + " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", + "\n", + " function canvas_keyboard_event(event) {\n", + " return fig.key_event(event, event['data']);\n", + " }\n", + "\n", + " canvas_div.keydown('key_press', canvas_keyboard_event);\n", + " canvas_div.keyup('key_release', canvas_keyboard_event);\n", + " this.canvas_div = canvas_div\n", + " this._canvas_extra_style(canvas_div)\n", + " this.root.append(canvas_div);\n", + "\n", + " var canvas = $('');\n", + " canvas.addClass('mpl-canvas');\n", + " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", + "\n", + " this.canvas = canvas[0];\n", + " this.context = canvas[0].getContext(\"2d\");\n", + "\n", + " var backingStore = this.context.backingStorePixelRatio ||\n", + "\tthis.context.webkitBackingStorePixelRatio ||\n", + "\tthis.context.mozBackingStorePixelRatio ||\n", + "\tthis.context.msBackingStorePixelRatio ||\n", + "\tthis.context.oBackingStorePixelRatio ||\n", + "\tthis.context.backingStorePixelRatio || 1;\n", + "\n", + " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", + "\n", + " var rubberband = $('');\n", + " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", + "\n", + " var pass_mouse_events = true;\n", + "\n", + " canvas_div.resizable({\n", + " start: function(event, ui) {\n", + " pass_mouse_events = false;\n", + " },\n", + " resize: function(event, ui) {\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " stop: function(event, ui) {\n", + " pass_mouse_events = true;\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " });\n", + "\n", + " function mouse_event_fn(event) {\n", + " if (pass_mouse_events)\n", + " return fig.mouse_event(event, event['data']);\n", + " }\n", + "\n", + " rubberband.mousedown('button_press', mouse_event_fn);\n", + " rubberband.mouseup('button_release', mouse_event_fn);\n", + " // Throttle sequential mouse events to 1 every 20ms.\n", + " rubberband.mousemove('motion_notify', mouse_event_fn);\n", + "\n", + " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", + " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", + "\n", + " canvas_div.on(\"wheel\", function (event) {\n", + " event = event.originalEvent;\n", + " event['data'] = 'scroll'\n", + " if (event.deltaY < 0) {\n", + " event.step = 1;\n", + " } else {\n", + " event.step = -1;\n", + " }\n", + " mouse_event_fn(event);\n", + " });\n", + "\n", + " canvas_div.append(canvas);\n", + " canvas_div.append(rubberband);\n", + "\n", + " this.rubberband = rubberband;\n", + " this.rubberband_canvas = rubberband[0];\n", + " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", + " this.rubberband_context.strokeStyle = \"#000000\";\n", + "\n", + " this._resize_canvas = function(width, height) {\n", + " // Keep the size of the canvas, canvas container, and rubber band\n", + " // canvas in synch.\n", + " canvas_div.css('width', width)\n", + " canvas_div.css('height', height)\n", + "\n", + " canvas.attr('width', width * mpl.ratio);\n", + " canvas.attr('height', height * mpl.ratio);\n", + " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", + "\n", + " rubberband.attr('width', width);\n", + " rubberband.attr('height', height);\n", + " }\n", + "\n", + " // Set the figure to an initial 600x600px, this will subsequently be updated\n", + " // upon first draw.\n", + " this._resize_canvas(600, 600);\n", + "\n", + " // Disable right mouse context menu.\n", + " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", + " return false;\n", + " });\n", + "\n", + " function set_focus () {\n", + " canvas.focus();\n", + " canvas_div.focus();\n", + " }\n", + "\n", + " window.setTimeout(set_focus, 100);\n", + "}\n", + "\n", + "mpl.figure.prototype._init_toolbar = function() {\n", + " var fig = this;\n", + "\n", + " var nav_element = $('
')\n", + " nav_element.attr('style', 'width: 100%');\n", + " this.root.append(nav_element);\n", + "\n", + " // Define a callback function for later on.\n", + " function toolbar_event(event) {\n", + " return fig.toolbar_button_onclick(event['data']);\n", + " }\n", + " function toolbar_mouse_event(event) {\n", + " return fig.toolbar_button_onmouseover(event['data']);\n", + " }\n", + "\n", + " for(var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " // put a spacer in here.\n", + " continue;\n", + " }\n", + " var button = $('');\n", + " button.click(method_name, toolbar_event);\n", + " button.mouseover(tooltip, toolbar_mouse_event);\n", + " nav_element.append(button);\n", + " }\n", + "\n", + " // Add the status bar.\n", + " var status_bar = $('');\n", + " nav_element.append(status_bar);\n", + " this.message = status_bar[0];\n", + "\n", + " // Add the close button to the window.\n", + " var buttongrp = $('
');\n", + " var button = $('');\n", + " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", + " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", + " buttongrp.append(button);\n", + " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", + " titlebar.prepend(buttongrp);\n", + "}\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(el){\n", + " var fig = this\n", + " el.on(\"remove\", function(){\n", + "\tfig.close_ws(fig, {});\n", + " });\n", + "}\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(el){\n", + " // this is important to make the div 'focusable\n", + " el.attr('tabindex', 0)\n", + " // reach out to IPython and tell the keyboard manager to turn it's self\n", + " // off when our div gets focus\n", + "\n", + " // location in version 3\n", + " if (IPython.notebook.keyboard_manager) {\n", + " IPython.notebook.keyboard_manager.register_events(el);\n", + " }\n", + " else {\n", + " // location in version 2\n", + " IPython.keyboard_manager.register_events(el);\n", + " }\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._key_event_extra = function(event, name) {\n", + " var manager = IPython.notebook.keyboard_manager;\n", + " if (!manager)\n", + " manager = IPython.keyboard_manager;\n", + "\n", + " // Check for shift+enter\n", + " if (event.shiftKey && event.which == 13) {\n", + " this.canvas_div.blur();\n", + " // select the cell after this one\n", + " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", + " IPython.notebook.select(index + 1);\n", + " }\n", + "}\n", + "\n", + "mpl.figure.prototype.handle_save = function(fig, msg) {\n", + " fig.ondownload(fig, null);\n", + "}\n", + "\n", + "\n", + "mpl.find_output_cell = function(html_output) {\n", + " // Return the cell and output element which can be found *uniquely* in the notebook.\n", + " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", + " // IPython event is triggered only after the cells have been serialised, which for\n", + " // our purposes (turning an active figure into a static one), is too late.\n", + " var cells = IPython.notebook.get_cells();\n", + " var ncells = cells.length;\n", + " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", + " data = data.data;\n", + " }\n", + " if (data['text/html'] == html_output) {\n", + " return [cell, data, j];\n", + " }\n", + " }\n", + " }\n", + " }\n", + "}\n", + "\n", + "// Register the function which deals with the matplotlib target/channel.\n", + "// The kernel may be null if the page has been refreshed.\n", + "if (IPython.notebook.kernel != null) {\n", + " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", + "}\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 24, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "distance_matrix = sp_distance.squareform(sp_distance.pdist(np.vstack([trials_A.transpose(), trials_B.transpose()])))\n", + "\n", + "mds = manifold.MDS(n_components=2, dissimilarity='precomputed') # Fit the mds object\n", + "coords = mds.fit(distance_matrix).embedding_ # Find the mds coordinates\n", + "plt.figure()\n", + "plt.scatter(coords[:, 0], coords[:, 1], c=['red'] * trials_A.shape[1] + ['green'] * trials_B.shape[1])\n", + "plt.axis('off')\n", + "plt.title('Low Dimensional Representation of conditions A and B')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*4.2 Test for univariate effect*\n", + "\n", + "Do a t test to compare the means of the voxels between these two conditions to determine if there is a difference" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Mean difference between condition A and B: -0.52\n", + "pvalue: 0.677\n" + ] + } + ], + "source": [ + "mean_difference = (np.mean(trials_A,0) - np.mean(trials_B,0))\n", + "ttest = stats.ttest_1samp(mean_difference, 0)\n", + "\n", + "print('Mean difference between condition A and B: ' + str(mean_difference.mean())[0:5])\n", + "print('pvalue: '+ str(ttest.pvalue)[0:5])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*4.3 Test for a multivariate effect*\n", + "\n", + "Use SVM from scikit-learn to estimate the classification accuracy between the conditions" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Classification accuracy between condition A and B: 0.833\n" + ] + } + ], + "source": [ + "input_mat = np.vstack([trials_A.transpose(), trials_B.transpose()])\n", + "input_labels = trials_A.shape[1] * [1] + trials_B.shape[1] * [0]\n", + "\n", + "\n", + "# Set up the classifier\n", + "X_train, X_test, y_train, y_test = sklearn.model_selection.train_test_split(\n", + " input_mat, input_labels, test_size=0.2, random_state=0)\n", + "\n", + "clf = sklearn.svm.SVC(kernel='linear', C=1).fit(X_train, y_train)\n", + "\n", + "score = clf.score(X_test, y_test)\n", + "print('Classification accuracy between condition A and B: ' + str(score)[0:5])" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "### **References**\n", + "Biswal, B., et al. (1996) Reduction of physiological fluctuations in fMRI using digital filters. Magnetic Resonance in Medicine 35, 107-113\n", + "\n", + "Friedman, L. and Glover, G.H. (2006) Report on a multicenter fMRI quality assurance protocol. Journal of Magnetic Resonance Imaging 23, 827-839\n", + "\n", + "Friston, K.J., et al. (1998) Event-related fMRI: characterizing differential responses. Neuroimage 7, 30-40\n", + "\n", + "Gudbjartsson, H. and Patz, S. (1995) The Rician distribution of noisy MRI data. Magnetic resonance in medicine 34, 910-914\n", + "\n", + "Welvaert, M., et al. (2011) neuRosim: An R package for generating fMRI data. Journal of Statistical Software 44, 1-18\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.1" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/tests/utils/test_fmrisim.py b/tests/utils/test_fmrisim.py index 2de4f8df5..9aa3834f8 100644 --- a/tests/utils/test_fmrisim.py +++ b/tests/utils/test_fmrisim.py @@ -75,17 +75,18 @@ def test_generate_stimfunction(): total_time=duration, ) - assert len(stimfunction) == duration * 1000, "stimfunction incorrect " \ - "length" + assert stimfunction.shape[0] == duration * 1000, "stimfunc incorrect " \ + "length" eventNumber = np.sum(event_durations * len(onsets)) * 1000 assert np.sum(stimfunction) == eventNumber, "Event number" # Create the signal function - signal_function = sim.double_gamma_hrf(stimfunction=stimfunction, - tr_duration=tr_duration, - ) - assert len(signal_function) == len(stimfunction) / (tr_duration * 1000), \ - "The length did not change" + signal_function = sim.convolve_hrf(stimfunction=stimfunction, + tr_duration=tr_duration, + ) + + stim_dur = stimfunction.shape[0] / (tr_duration * 1000) + assert signal_function.shape[0] == stim_dur, "The length did not change" onsets = [10] stimfunction = sim.generate_stimfunction(onsets=onsets, @@ -93,9 +94,9 @@ def test_generate_stimfunction(): total_time=duration, ) - signal_function = sim.double_gamma_hrf(stimfunction=stimfunction, - tr_duration=tr_duration, - ) + signal_function = sim.convolve_hrf(stimfunction=stimfunction, + tr_duration=tr_duration, + ) assert np.sum(signal_function < 0) > 0, "No values below zero" @@ -128,13 +129,13 @@ def test_apply_signal(): total_time=duration, ) - signal_function = sim.double_gamma_hrf(stimfunction=stimfunction, - tr_duration=tr_duration, - ) + signal_function = sim.convolve_hrf(stimfunction=stimfunction, + tr_duration=tr_duration, + ) # Convolve the HRF with the stimulus sequence signal = sim.apply_signal(signal_function=signal_function, - volume_static=volume, + volume_signal=volume, ) assert signal.shape == (dimensions[0], dimensions[1], dimensions[2], @@ -142,7 +143,7 @@ def test_apply_signal(): "wrong size" signal = sim.apply_signal(signal_function=stimfunction, - volume_static=volume, + volume_signal=volume, ) assert np.any(signal == signal_magnitude), "The stimfunction is not binary" @@ -177,25 +178,27 @@ def test_generate_noise(): total_time=duration, ) - signal_function = sim.double_gamma_hrf(stimfunction=stimfunction, - tr_duration=tr_duration, - ) + signal_function = sim.convolve_hrf(stimfunction=stimfunction, + tr_duration=tr_duration, + ) # Convolve the HRF with the stimulus sequence signal = sim.apply_signal(signal_function=signal_function, - volume_static=volume, + volume_signal=volume, ) # Generate the mask of the signal - mask = sim.mask_brain(signal, mask_threshold=0.1) + mask, template = sim.mask_brain(signal, mask_threshold=0.1) assert min(mask[mask > 0]) > 0.1, "Mask thresholding did not work" + assert len(np.unique(template) > 2), "Template creation did not work" stimfunction_tr = stimfunction[::int(tr_duration * 1000)] # Create the noise volumes (using the default parameters) noise = sim.generate_noise(dimensions=dimensions, stimfunction_tr=stimfunction_tr, tr_duration=tr_duration, + template=template, mask=mask, ) @@ -207,13 +210,14 @@ def test_generate_noise(): noise = sim.generate_noise(dimensions=dimensions, stimfunction_tr=stimfunction_tr, tr_duration=tr_duration, + template=template, mask=mask, - noise_dict={'temporal_noise': 0, 'sfnr': 10000}, + noise_dict={'sfnr': 10000, 'snr': 10000}, ) - temporal_noise = np.std(noise[mask[:, :, :, 0] > 0], 1).mean() + system_noise = np.std(noise[mask > 0], 1).mean() - assert temporal_noise <= 0.1, "Noise strength could not be manipulated" + assert system_noise <= 0.1, "Noise strength could not be manipulated" def test_mask_brain(): @@ -235,8 +239,8 @@ def test_mask_brain(): ) # Mask the volume to be the same shape as a brain - mask = sim.mask_brain(volume)[:, :, :, 0] - brain = volume * (mask > 0) + mask, _ = sim.mask_brain(volume) + brain = volume * mask assert np.sum(brain != 0) == np.sum(volume != 0), "Masking did not work" assert brain[0, 0, 0] == 0, "Masking did not work" @@ -253,8 +257,8 @@ def test_mask_brain(): ) # Mask the volume to be the same shape as a brain - mask = sim.mask_brain(volume)[:, :, :, 0] - brain = volume * (mask > 0) + mask, _ = sim.mask_brain(volume) + brain = volume * mask assert np.sum(brain != 0) < np.sum(volume != 0), "Masking did not work" @@ -272,7 +276,7 @@ def test_calc_noise(): # Preset the noise dict nd_orig = {'auto_reg_sigma': 0.6, 'drift_sigma': 0.4, - 'temporal_noise': 5, + 'snr': 30, 'sfnr': 30, 'max_activity': 1000, 'fwhm': 4, @@ -285,21 +289,36 @@ def test_calc_noise(): ) # Mask the volume to be the same shape as a brain - mask = sim.mask_brain(dimensions_tr) + mask, template = sim.mask_brain(dimensions_tr, mask_threshold=0.2) stimfunction_tr = stimfunction[::int(tr_duration * 1000)] noise = sim.generate_noise(dimensions=dimensions_tr[0:3], stimfunction_tr=stimfunction_tr, tr_duration=tr_duration, + template=template, mask=mask, noise_dict=nd_orig, ) + # Check that noise_system is being calculated correctly + spatial_sd = 5 + temporal_sd = 5 + noise_system = sim._generate_noise_system(dimensions_tr, + spatial_sd, + temporal_sd) + + precision = abs(noise_system[0, 0, 0, :].std() - spatial_sd) + assert precision < spatial_sd, 'noise_system calculated incorrectly' + + precision = abs(noise_system[:, :, :, 0].std() - temporal_sd) + assert precision < spatial_sd, 'noise_system calculated incorrectly' + # Calculate the noise - nd_calc = sim.calc_noise(noise, mask) + nd_calc = sim.calc_noise(volume=noise, + mask=mask) # How precise are these estimates - precision = abs(nd_calc['temporal_noise'] - nd_orig['temporal_noise']) - assert precision < 1, 'temporal_noise calculated incorrectly' + precision = abs(nd_calc['snr'] - nd_orig['snr']) + assert precision < nd_orig['snr'], 'snr calculated incorrectly' precision = abs(nd_calc['sfnr'] - nd_orig['sfnr']) - assert precision < 5, 'sfnr calculated incorrectly' + assert precision < nd_orig['sfnr'], 'sfnr calculated incorrectly'