diff --git a/brainiak/utils/fmrisim.py b/brainiak/utils/fmrisim.py index af3857feb..a079abd6b 100644 --- a/brainiak/utils/fmrisim.py +++ b/brainiak/utils/fmrisim.py @@ -2217,7 +2217,8 @@ def mask_brain(volume, if mask_self is True: mask_raw = volume elif template_name is None: - mask_raw = np.load(resource_stream(__name__, "grey_matter_mask.npy")) + mfn = resource_stream(__name__, "sim_parameters/grey_matter_mask.npy") + mask_raw = np.load(mfn) else: mask_raw = np.load(template_name) diff --git a/brainiak/utils/fmrisim_real_time_generator.py b/brainiak/utils/fmrisim_real_time_generator.py new file mode 100644 index 000000000..c3b20912a --- /dev/null +++ b/brainiak/utils/fmrisim_real_time_generator.py @@ -0,0 +1,605 @@ +# Generate simulated fMRI data with a few parameters that might be relevant +# for real time analysis +""" +This code can be run as a function in python or from the command line: +python fmrisim_real-time_generator --inputDir fmrisim_files/ --outputDir data + +The input arguments are: +Required: +inputDir - Specify input data dir where the parameters for fmrisim are +outputDir - Specify output data dir where the data should be saved + +Optional (can be modified by flags from the command line): +data_dict contains: +numTRs - Specify the number of time points +multivariate_patterns - Is the difference between conditions univariate (0) +or multivariate (1) +different_ROIs - Are there different ROIs for each condition (1) or is +it in the same ROI (0). If it is the same ROI and you are using univariate +differences, the second condition will have a smaller evoked response than +the other. +event_duration - How long, in seconds, is each event +scale_percentage - What is the percent signal change +trDuration - How many seconds per volume +save_dicom - Do you want to save data as a dicom (1) or numpy (0) +save_realtime - Do you want to save the data in real time (1) or as fast as +possible (0)? +isi - What is the time between each event (in seconds) +burn_in - How long before the first event (in seconds) +""" +import os +import time +import argparse +import datetime +import nibabel # type: ignore +import numpy as np # type: ignore +import pydicom as dicom +from brainiak.utils import fmrisim as sim # type: ignore +import logging +from pkg_resources import resource_stream +from nibabel.nifti1 import Nifti1Image +import gzip + +__all__ = ["generate_data"] + +logger = logging.getLogger(__name__) + +script_datetime = datetime.datetime.now() + + +def _generate_ROIs(ROI_file, + stimfunc, + noise, + scale_percentage, + data_dict): + """Make signal activity for an ROI of data + Creates the specified evoked response time course, calibrated to the + expected signal change, for a given ROI + + Parameters + ---------- + + ROI_file : str + Path to the file of the ROI being loaded in + + stimfunc : 1 dimensional array + Time course of evoked response. Output from + fmrisim.generate_stimfunction + + noise : 4 dimensional array + Volume of noise generated from fmrisim.generate_noise. Although this + is needed as an input, this is only so that the percent signal change + can be calibrated. This is not combined with the signal generated. + + scale_percentage : float + What is the percent signal change for the evoked response + + data_dict : dict + A dictionary to specify the parameters used for making data, + specifying the following keys + numTRs - int - Specify the number of time points + multivariate_patterns - bool - Is the difference between conditions + univariate (0) or multivariate (1) + different_ROIs - bool - Are there different ROIs for each condition ( + 1) or is it in the same ROI (0). If it is the same ROI and you are + using univariate differences, the second condition will have a + smaller evoked response than the other. + event_duration - int - How long, in seconds, is each event + scale_percentage - float - What is the percent signal change + trDuration - float - How many seconds per volume + save_dicom - bool - Save to data as a dicom (1) or numpy (0) + save_realtime - bool - Do you want to save the data in real time (1) + or as fast as possible (0)? + isi - float - What is the time between each event (in seconds) + burn_in - int - How long before the first event (in seconds) + + Returns + ---------- + + signal : 4 dimensional array + Volume of signal in the specified ROI (noise has not yet been added) + + """ + + # Create the signal in the ROI as specified. + + logger.info('Loading', ROI_file) + + # Load in the template data (it may already be loaded if doing a test) + if isinstance(ROI_file, str): + nii = nibabel.load(ROI_file) + ROI = nii.get_data() + else: + ROI = ROI_file + + # Find all the indices that contain signal + idx_list = np.where(ROI == 1) + + idxs = np.zeros([len(idx_list[0]), 3]) + for idx_counter in list(range(0, len(idx_list[0]))): + idxs[idx_counter, 0] = int(idx_list[0][idx_counter]) + idxs[idx_counter, 1] = int(idx_list[1][idx_counter]) + idxs[idx_counter, 2] = int(idx_list[2][idx_counter]) + + idxs = idxs.astype('int8') + + # How many voxels per ROI + voxels = int(ROI.sum()) + + # Create a pattern of activity across the two voxels + if data_dict['multivariate_pattern'] is True: + pattern = np.random.rand(voxels).reshape((voxels, 1)) + else: # Just make a univariate increase + pattern = np.tile(1, voxels).reshape((voxels, 1)) + + # Multiply each pattern by each voxel time course + weights = np.tile(stimfunc, voxels) * pattern.T + + # Convolve the onsets with the HRF + temporal_res = 1 / data_dict['trDuration'] + signal_func = sim.convolve_hrf(stimfunction=weights, + tr_duration=data_dict['trDuration'], + temporal_resolution=temporal_res, + scale_function=1, + ) + + # Change the type of noise + noise = noise.astype('double') + + # Create a noise function (same voxels for signal function as for noise) + noise_function = noise[idxs[:, 0], idxs[:, 1], idxs[:, 2], :].T + + # Compute the signal magnitude for the data + sf_scaled = sim.compute_signal_change(signal_function=signal_func, + noise_function=noise_function, + noise_dict=data_dict['noise_dict'], + magnitude=[scale_percentage], + method='PSC', + ) + + # Combine the signal time course with the signal volume + signal = sim.apply_signal(sf_scaled, + ROI, + ) + + # Return signal needed + return signal + + +def _write_dicom(output_name, + data, + image_number=0): + """Write the data to a dicom file + Saves the data for one TR to a dicom. + + Dicom files are difficult to set up correctly, this file will likely + crash when trying to open it using dcm2nii. However, if it is loaded in + python (e.g., dicom.dcmread) then pixel_array contains the relevant + voxel data + + Parameters + ---------- + + output_name : str + Output name for volume being created + + data : 3 dimensional array + Volume of data to be saved + + image_number : int + Number dicom to be saved. This is critical for setting up dicom file + header information. + + """ + + # Convert data from float to in + dataInts = data.astype(np.int16) + + # Populate required values for file meta information + file_meta = dicom.Dataset() + file_meta.MediaStorageSOPClassUID = '1.2' # '1.2.840.10008.5.1.4.1.1.2' + file_meta.MediaStorageSOPInstanceUID = "1.2.3" + file_meta.ImplementationClassUID = "1.2.3.4" + file_meta.TransferSyntaxUID = '1.2.840.10008.1.2' + + # Create the FileDataset + ds = dicom.FileDataset(output_name, + {}, + file_meta=file_meta, + preamble=b"\0" * 128) + + # Set image dimensions + frames, rows, cols = dataInts.shape + ds.Rows = rows + ds.Columns = cols + ds.NumberOfFrames = frames + ds.SamplesPerPixel = 1 + ds.BitsAllocated = 16 + ds.BitsStored = 16 + ds.PixelRepresentation = 0 + ds.InstanceNumber = image_number + ds.ImagePositionPatient = [0, 0, 0] + ds.ImageOrientationPatient = [.01, 0, 0, 0, 0, 0] + + # Add the data elements -- not trying to set all required here. Check DICOM + # standard + ds.PatientName = "sim" + ds.PatientID = "sim" + + # Set the transfer syntax + ds.is_little_endian = True + ds.is_implicit_VR = True + + # Set creation date/time + image_datetime = script_datetime + datetime.timedelta(seconds=image_number) + timeStr = image_datetime.strftime('%H%M%S') + ds.ContentDate = image_datetime.strftime('%Y%m%d') + ds.ContentTime = timeStr + + # Add the data + ds.PixelData = dataInts.tobytes() + + ds.save_as(output_name) + + +def _get_input_names(data_dict): + """Get names from dict + Read in the data_dict to return the relevant file names + + Parameters + ---------- + + data_dict : dict + A dictionary to specify the parameters used for making data, + specifying the following keys + numTRs - int - Specify the number of time points + multivariate_patterns - bool - Is the difference between conditions + univariate (0) or multivariate (1) + different_ROIs - bool - Are there different ROIs for each condition ( + 1) or is it in the same ROI (0). If it is the same ROI and you are + using univariate differences, the second condition will have a + smaller evoked response than the other. + event_duration - int - How long, in seconds, is each event + scale_percentage - float - What is the percent signal change + trDuration - float - How many seconds per volume + save_dicom - bool - Save to data as a dicom (1) or numpy (0) + save_realtime - bool - Do you want to save the data in real time (1) + or as fast as possible (0)? + isi - float - What is the time between each event (in seconds) + burn_in - int - How long before the first event (in seconds) + + Returns + ---------- + + ROI_A_file : str + Path to ROI for condition A + + ROI_B_file : str + Path to ROI for condition B + + template_path : str + Path to template file for data + + noise_dict_file : str + Path to file containing parameters for noise simulation + + """ + + # Load in the ROIs + if data_dict['ROI_A_file'] is None: + vol = resource_stream(__name__, "sim_parameters/ROI_A.nii.gz").read() + ROI_A_file = Nifti1Image.from_bytes(gzip.decompress(vol)).get_data() + else: + ROI_A_file = data_dict['ROI_A_file'] + + if data_dict['ROI_B_file'] is None: + vol = resource_stream(__name__, "sim_parameters/ROI_B.nii.gz").read() + ROI_B_file = Nifti1Image.from_bytes(gzip.decompress(vol)).get_data() + else: + ROI_B_file = data_dict['ROI_B_file'] + + # Get the path to the template + if data_dict['template_path'] is None: + vol = resource_stream(__name__, + "sim_parameters/sub_template.nii.gz").read() + template_path = Nifti1Image.from_bytes(gzip.decompress(vol)).get_data() + else: + template_path = data_dict['template_path'] + + # Load in the noise dict if supplied + if data_dict['noise_dict_file'] is None: + file = resource_stream(__name__, + 'sim_parameters/sub_noise_dict.txt').read() + noise_dict_file = file + else: + noise_dict_file = data_dict['noise_dict_file'] + + # Return the paths + return ROI_A_file, ROI_B_file, template_path, noise_dict_file + + +def generate_data(outputDir, + data_dict): + """Generate simulated fMRI data + Use a few parameters that might be relevant for real time analysis + + Parameters + ---------- + + inputDir : str + Specify input data dir where the parameters for fmrisim are + + outputDir : str + Specify output data dir where the data should be saved + + data_dict : dict + A dictionary to specify the parameters used for making data, + specifying the following keys + numTRs - int - Specify the number of time points + multivariate_patterns - bool - Is the difference between conditions + univariate (0) or multivariate (1) + different_ROIs - bool - Are there different ROIs for each condition ( + 1) or is it in the same ROI (0). If it is the same ROI and you are + using univariate differences, the second condition will have a + smaller evoked response than the other. + event_duration - int - How long, in seconds, is each event + scale_percentage - float - What is the percent signal change + trDuration - float - How many seconds per volume + save_dicom - bool - Save to data as a dicom (1) or numpy (0) + save_realtime - bool - Do you want to save the data in real time (1) + or as fast as possible (0)? + isi - float - What is the time between each event (in seconds) + burn_in - int - How long before the first event (in seconds) + + """ + + # If the folder doesn't exist then make it + os.system('mkdir -p %s' % outputDir) + + logger.info('Load template of average voxel value') + + # Get the file names needed for loading in the data + ROI_A_file, ROI_B_file, template_path, noise_dict_file = _get_input_names( + data_dict) + + # Load in the template data (it may already be loaded if doing a test) + if isinstance(template_path, str): + template_nii = nibabel.load(template_path) + template = template_nii.get_data() + else: + template = template_path + + dimensions = np.array(template.shape[0:3]) + + logger.info('Create binary mask and normalize the template range') + mask, template = sim.mask_brain(volume=template, + mask_self=True, + ) + + # Write out the mask as a numpy file + outFile = os.path.join(outputDir, 'mask.npy') + np.save(outFile, mask.astype(np.uint8)) + + # Load the noise dictionary + logger.info('Loading noise parameters') + + # If this isn't a string, assume it is a resource stream file + if type(noise_dict_file) is str: + with open(noise_dict_file, 'r') as f: + noise_dict = f.read() + else: + # Read the resource stream object + noise_dict = noise_dict_file.decode() + + noise_dict = eval(noise_dict) + noise_dict['matched'] = 0 # Increases processing time + + # Add it here for easy access + data_dict['noise_dict'] = data_dict + + logger.info('Generating noise') + temp_stimfunction = np.zeros((data_dict['numTRs'], 1)) + noise = sim.generate_noise(dimensions=dimensions, + stimfunction_tr=temp_stimfunction, + tr_duration=int(data_dict['trDuration']), + template=template, + mask=mask, + noise_dict=noise_dict, + ) + + # Create the stimulus time course of the conditions + total_time = int(data_dict['numTRs'] * data_dict['trDuration']) + onsets_A = [] + onsets_B = [] + curr_time = data_dict['burn_in'] + while curr_time < (total_time - data_dict['event_duration']): + + # Flip a coin for each epoch to determine whether it is A or B + if np.random.randint(0, 2) == 1: + onsets_A.append(curr_time) + else: + onsets_B.append(curr_time) + + # Increment the current time + curr_time += data_dict['event_duration'] + data_dict['isi'] + + # How many timepoints per second of the stim function are to be generated? + temporal_res = 1 / data_dict['trDuration'] + + # Create a time course of events + event_durations = [data_dict['event_duration']] + stimfunc_A = sim.generate_stimfunction(onsets=onsets_A, + event_durations=event_durations, + total_time=total_time, + temporal_resolution=temporal_res, + ) + + stimfunc_B = sim.generate_stimfunction(onsets=onsets_B, + event_durations=event_durations, + total_time=total_time, + temporal_resolution=temporal_res, + ) + + # Create a labels timecourse + outFile = os.path.join(outputDir, 'labels.npy') + np.save(outFile, (stimfunc_A + (stimfunc_B * 2))) + + # How is the signal implemented in the different ROIs + signal_A = _generate_ROIs(ROI_A_file, + stimfunc_A, + noise, + data_dict['scale_percentage'], + data_dict) + if data_dict['different_ROIs'] is True: + + signal_B = _generate_ROIs(ROI_B_file, + stimfunc_B, + noise, + data_dict['scale_percentage'], + data_dict) + + else: + + # Halve the evoked response if these effects are both expected in the + # same ROI + if data_dict['multivariate_pattern'] is False: + signal_B = _generate_ROIs(ROI_A_file, + stimfunc_B, + noise, + data_dict['scale_percentage'] * 0.5, + data_dict) + else: + signal_B = _generate_ROIs(ROI_A_file, + stimfunc_B, + noise, + data_dict['scale_percentage'], + data_dict) + + # Combine the two signal timecourses + signal = signal_A + signal_B + + logger.info('Generating TRs in real time') + for idx in range(data_dict['numTRs']): + + # Create the brain volume on this TR + brain = noise[:, :, :, idx] + signal[:, :, :, idx] + + # Convert file to integers to mimic what you get from MR + brain_int32 = brain.astype(np.int32) + + # Store as dicom or nifti? + if data_dict['save_dicom'] is True: + # Save the volume as a DICOM file, with each TR as its own file + output_file = os.path.join(outputDir, 'rt_' + format(idx, '03d') + + '.dcm') + _write_dicom(output_file, brain_int32, idx+1) + else: + # Save the volume as a numpy file, with each TR as its own file + output_file = os.path.join(outputDir, 'rt_' + format(idx, '03d') + + '.npy') + np.save(output_file, brain_int32) + + logger.info("Generate {}".format(output_file)) + + # Sleep until next TR + if data_dict['save_realtime'] == 1: + time.sleep(data_dict['trDuration']) + + +if __name__ == '__main__': + # Receive the inputs + argParser = argparse.ArgumentParser( + 'Specify input arguments. Some arguments are parameters that require ' + 'an input is provided (noted by "Param"), others are flags that when ' + 'provided will change according to the flag (noted by "Flag")') + argParser.add_argument('--output-dir', '-o', default=None, type=str, + help='Param. Output directory for simulated data') + argParser.add_argument('--ROI-A-file', default=None, type=str, + help='Param. Full path to file for cond. A ROI') + argParser.add_argument('--ROI-B-file', default=None, type=str, + help='Param. Full path to file for cond. B ROI') + argParser.add_argument('--template-path', default=None, type=str, + help='Param. Full path to file for brain template') + argParser.add_argument('--noise-dict-file', default=None, type=str, + help='Param. Full path to file setting noise ' + 'params') + argParser.add_argument('--numTRs', '-n', default=200, type=int, + help='Param. Number of time points') + argParser.add_argument('--event-duration', '-d', default=10, type=int, + help='Param. Number of seconds per event') + argParser.add_argument('--scale-percentage', '-s', default=0.5, type=float, + help='Param. Percent signal change') + argParser.add_argument('--multivariate-pattern', '-m', default=False, + action='store_true', + help='Flag. Signal is different between conditions ' + 'in a multivariate, versus univariate, way') + argParser.add_argument('--different-ROIs', '-r', default=False, + action='store_true', help='Flag. Use different ' + 'ROIs for each condition') + argParser.add_argument('--save-dicom', default=False, + action='store_true', help='Flag. Output files in ' + 'DICOM format rather ' + 'than numpy') + argParser.add_argument('--save-realtime', default=False, + action='store_true', help='Flag. Save data as if ' + 'it was coming in at ' + 'the acquisition rate') + args = argParser.parse_args() + + # Essential arguments + outputDir = args.output_dir + + if outputDir is None: + logger.info("Must specify an output directory using -o") + exit(-1) + + data_dict = {} + + # User controlled settings + + # Specify the path to the files used for defining ROIs. + data_dict['ROI_A_file'] = args.ROI_A_file + data_dict['ROI_B_file'] = args.ROI_B_file + + # Specify where the template + data_dict['template_path'] = args.template_path + data_dict['noise_dict_file'] = args.noise_dict_file + + # Specify the number of time points + data_dict['numTRs'] = args.numTRs + + # How long is each event/block you are modelling (assumes 6s rest between) + data_dict['event_duration'] = float(args.event_duration) + + # What is the percent signal change being simulated + data_dict['scale_percentage'] = args.scale_percentage + + # Are there different ROIs for each condition (True) or is it in the same + # ROI (False). + # If it is the same ROI and you are using univariate differences, + # the second condition will have a smaller evoked response than the other. + data_dict['different_ROIs'] = args.different_ROIs + + # Is this a multivariate pattern (1) or a univariate pattern + data_dict['multivariate_pattern'] = args.multivariate_pattern + + # Do you want to save data as a dicom (True) or numpy (False) + data_dict['save_dicom'] = args.save_dicom + + # Do you want to save the data in real time (1) or as fast as possible (0)? + data_dict['save_realtime'] = args.save_realtime + + # Default settings + + # How long does each acquisition take + data_dict['trDuration'] = 2 + + # What is the time between each event (in seconds) + data_dict['isi'] = 6 + + # How long before the first event (in seconds) + data_dict['burn_in'] = 6 + + # Run the function if running from command line + generate_data(outputDir, + data_dict) diff --git a/examples/utils/sim_parameters/ROI_A.nii.gz b/brainiak/utils/sim_parameters/ROI_A.nii.gz similarity index 100% rename from examples/utils/sim_parameters/ROI_A.nii.gz rename to brainiak/utils/sim_parameters/ROI_A.nii.gz diff --git a/examples/utils/sim_parameters/ROI_B.nii.gz b/brainiak/utils/sim_parameters/ROI_B.nii.gz similarity index 100% rename from examples/utils/sim_parameters/ROI_B.nii.gz rename to brainiak/utils/sim_parameters/ROI_B.nii.gz diff --git a/brainiak/utils/grey_matter_mask.npy b/brainiak/utils/sim_parameters/grey_matter_mask.npy similarity index 100% rename from brainiak/utils/grey_matter_mask.npy rename to brainiak/utils/sim_parameters/grey_matter_mask.npy diff --git a/examples/utils/sim_parameters/sub_noise_dict.txt b/brainiak/utils/sim_parameters/sub_noise_dict.txt similarity index 100% rename from examples/utils/sim_parameters/sub_noise_dict.txt rename to brainiak/utils/sim_parameters/sub_noise_dict.txt diff --git a/examples/utils/sim_parameters/sub_template.nii.gz b/brainiak/utils/sim_parameters/sub_template.nii.gz similarity index 100% rename from examples/utils/sim_parameters/sub_template.nii.gz rename to brainiak/utils/sim_parameters/sub_template.nii.gz diff --git a/examples/utils/fmrisim_real-time_generator.py b/examples/utils/fmrisim_real-time_generator.py deleted file mode 100644 index bf7dd782d..000000000 --- a/examples/utils/fmrisim_real-time_generator.py +++ /dev/null @@ -1,412 +0,0 @@ -# Generate simulated fMRI data with a few parameters that might be relevant -# for real time analysis -# This code can be run as a function in python or from the command line: -# python fmrisim_real-time_generator --inputDir fmrisim_files/ --outputDir data/ -# -# The input arguments are: -# Required: -# inputDir - Specify input data dir where the parameters for fmrisim are -# outputDir - Specify output data dir where the data should be saved -# -# Optional (can be modified by flags from the command line): -# data_dict contains: -# numTRs - Specify the number of time points -# multivariate_patterns - Is the difference between conditions univariate -# (0) or multivariate (1) -# different_ROIs - Are there different ROIs for each condition (1) or is -# it in the same ROI (0). If it is the same ROI and you are using univariate -# differences, the second condition will have a smaller evoked response than -# the other. -# event_duration - How long, in seconds, is each event -# scale_percentage - What is the percent signal change -# trDuration - How many seconds per volume -# save_dicom - Do you want to save data as a dicom (1) or numpy (0) -# save_realtime - Do you want to save the data in real time (1) or as -# fast as possible (0)? -# isi - What is the time between each event (in seconds) -# burn_in - How long before the first event (in seconds) - -import os -import glob -import time -import random -import argparse -import datetime -import nibabel # type: ignore -import numpy as np # type: ignore -import pydicom as dicom -from brainiak.utils import fmrisim as sim # type: ignore -import sys - -script_datetime = datetime.datetime.now() - -def generate_ROIs(ROI_file, - stimfunc, - noise, - scale_percentage, - data_dict): - # Create the signal in the ROI as specified. - - print('Loading', ROI_file) - - nii = nibabel.load(ROI_file) - ROI = nii.get_data() - - # Find all the indices that contain signal - idx_list = np.where(ROI == 1) - - idxs = np.zeros([len(idx_list[0]), 3]) - for idx_counter in list(range(0, len(idx_list[0]))): - idxs[idx_counter, 0] = int(idx_list[0][idx_counter]) - idxs[idx_counter, 1] = int(idx_list[1][idx_counter]) - idxs[idx_counter, 2] = int(idx_list[2][idx_counter]) - - idxs = idxs.astype('int8') - - # How many voxels per ROI - voxels = int(ROI.sum()) - - # Create a pattern of activity across the two voxels - if data_dict['multivariate_pattern'] is True: - pattern = np.random.rand(voxels).reshape((voxels, 1)) - else: # Just make a univariate increase - pattern = np.tile(1, voxels).reshape((voxels, 1)) - - # Multiply each pattern by each voxel time course - weights = np.tile(stimfunc, voxels) * pattern.T - - # Convolve the onsets with the HRF - temporal_res = 1 / data_dict['trDuration'] - signal_func = sim.convolve_hrf(stimfunction=weights, - tr_duration=data_dict['trDuration'], - temporal_resolution=temporal_res, - scale_function=1, - ) - - # Change the type of noise - noise = noise.astype('double') - - # Create a noise function (same voxels for signal function as for noise) - noise_function = noise[idxs[:, 0], idxs[:, 1], idxs[:, 2], :].T - - # Compute the signal magnitude for the data - sf_scaled = sim.compute_signal_change(signal_function=signal_func, - noise_function=noise_function, - noise_dict=data_dict['noise_dict'], - magnitude=[scale_percentage], - method='PSC', - ) - - # Combine the signal time course with the signal volume - signal = sim.apply_signal(sf_scaled, - ROI, - ) - - # Return signal needed - return signal - - -def write_dicom(output_name, - data, - image_number=0): - # Write the data to a dicom file. - # Dicom files are difficult to set up correctly, this file will likely - # crash when trying to open it using dcm2nii. However, if it is loaded in - # python (e.g., dicom.dcmread) then pixel_array contains the relevant - # voxel data - - # Convert data from float to in - dataInts = data.astype(np.int16) - - # Populate required values for file meta information - file_meta = dicom.Dataset() - file_meta.MediaStorageSOPClassUID = '1.2' # '1.2.840.10008.5.1.4.1.1.2' - file_meta.MediaStorageSOPInstanceUID = "1.2.3" - file_meta.ImplementationClassUID = "1.2.3.4" - file_meta.TransferSyntaxUID = '1.2.840.10008.1.2' - - # Create the FileDataset - ds = dicom.FileDataset(output_name, - {}, - file_meta=file_meta, - preamble=b"\0" * 128) - - # Set image dimensions - frames, rows, cols = dataInts.shape - ds.Rows = rows - ds.Columns = cols - ds.NumberOfFrames = frames - ds.SamplesPerPixel = 1 - ds.BitsAllocated = 16 - ds.BitsStored = 16 - ds.PixelRepresentation = 0 - ds.InstanceNumber = image_number - ds.ImagePositionPatient = [0, 0, 0] - ds.ImageOrientationPatient = [.01, 0, 0, 0, 0, 0] - - # Add the data elements -- not trying to set all required here. Check DICOM - # standard - ds.PatientName = "sim" - ds.PatientID = "sim" - - # Set the transfer syntax - ds.is_little_endian = True - ds.is_implicit_VR = True - - # Set creation date/time - image_datetime = script_datetime + datetime.timedelta(seconds=image_number) - timeStr = image_datetime.strftime('%H%M%S') - ds.ContentDate = image_datetime.strftime('%Y%m%d') - ds.ContentTime = timeStr - - # Add the data - ds.PixelData = dataInts.tobytes() - - ds.save_as(output_name) - - -def generate_data(inputDir, - outputDir, - data_dict): - # Generate simulated fMRI data with a few parameters that might be - # relevant for real time analysis - # inputDir - Specify input data dir where the parameters for fmrisim are - # outputDir - Specify output data dir where the data should be saved - # data_dict contains: - # numTRs - Specify the number of time points - # multivariate_patterns - Is the difference between conditions - # univariate (0) or multivariate (1) - # different_ROIs - Are there different ROIs for each condition (1) or - # is it in the same ROI (0). If it is the same ROI and you are using - # univariate differences, the second condition will have a smaller evoked - # response than the other. - # event_duration - How long, in seconds, is each event - # scale_percentage - What is the percent signal change - # trDuration - How many seconds per volume - # save_dicom - Do you want to save data as a dicom (1) or numpy (0) - # save_realtime - Do you want to save the data in real time (1) or as - # fast as possible (0)? - # isi - What is the time between each event (in seconds) - # burn_in - How long before the first event (in seconds) - - # If the folder doesn't exist then make it - if os.path.isdir(outputDir) is False: - os.makedirs(outputDir, exist_ok=True) - - print('Load template of average voxel value') - templateFile = os.path.join(inputDir, 'sub_template.nii.gz') - template_nii = nibabel.load(templateFile) - template = template_nii.get_data() - - dimensions = np.array(template.shape[0:3]) - - print('Create binary mask and normalize the template range') - mask, template = sim.mask_brain(volume=template, - mask_self=True, - ) - - # Write out the mask as a numpy file - outFile = os.path.join(outputDir, 'mask.npy') - np.save(outFile, mask.astype(np.uint8)) - - # Load the noise dictionary - print('Loading noise parameters') - noiseFile = os.path.join(inputDir, 'sub_noise_dict.txt') - with open(noiseFile, 'r') as f: - noise_dict = f.read() - noise_dict = eval(noise_dict) - noise_dict['matched'] = 0 # Increases processing time - - # Add it here for easy access - data_dict['noise_dict'] = data_dict - - print('Generating noise') - temp_stimfunction = np.zeros((data_dict['numTRs'], 1)) - noise = sim.generate_noise(dimensions=dimensions, - stimfunction_tr=temp_stimfunction, - tr_duration=int(data_dict['trDuration']), - template=template, - mask=mask, - noise_dict=noise_dict, - ) - - # Create the stimulus time course of the conditions - total_time = int(data_dict['numTRs'] * data_dict['trDuration']) - onsets_A = [] - onsets_B = [] - curr_time = data_dict['burn_in'] - while curr_time < (total_time - data_dict['event_duration']): - - # Flip a coin for each epoch to determine whether it is A or B - if np.random.randint(0, 2) == 1: - onsets_A.append(curr_time) - else: - onsets_B.append(curr_time) - - # Increment the current time - curr_time += data_dict['event_duration'] + data_dict['isi'] - - # How many timepoints per second of the stim function are to be generated? - temporal_res = 1 / data_dict['trDuration'] - - # Create a time course of events - event_durations = [data_dict['event_duration']] - stimfunc_A = sim.generate_stimfunction(onsets=onsets_A, - event_durations=event_durations, - total_time=total_time, - temporal_resolution=temporal_res, - ) - - stimfunc_B = sim.generate_stimfunction(onsets=onsets_B, - event_durations=event_durations, - total_time=total_time, - temporal_resolution=temporal_res, - ) - - # Create a labels timecourse - outFile = os.path.join(outputDir, 'labels.npy') - np.save(outFile, (stimfunc_A + (stimfunc_B * 2))) - - roiA_file = os.path.join(inputDir, 'ROI_A.nii.gz') - roiB_file = os.path.join(inputDir, 'ROI_B.nii.gz') - - # How is the signal implemented in the different ROIs - signal_A = generate_ROIs(roiA_file, - stimfunc_A, - noise, - data_dict['scale_percentage'], - data_dict) - if data_dict['different_ROIs'] is True: - - signal_B = generate_ROIs(roiB_file, - stimfunc_B, - noise, - data_dict['scale_percentage'], - data_dict) - - else: - - # Halve the evoked response if these effects are both expected in the same ROI - if data_dict['multivariate_pattern'] is False: - signal_B = generate_ROIs(roiA_file, - stimfunc_B, - noise, - data_dict['scale_percentage'] * 0.5, - data_dict) - else: - signal_B = generate_ROIs(roiA_file, - stimfunc_B, - noise, - data_dict['scale_percentage'], - data_dict) - - # Combine the two signal timecourses - signal = signal_A + signal_B - - print('Generating TRs in real time') - for idx in range(data_dict['numTRs']): - - # Create the brain volume on this TR - brain = noise[:, :, :, idx] + signal[:, :, :, idx] - - # Convert file to integers to mimic what you get from MR - brain_int32 = brain.astype(np.int32) - - # Store as dicom or nifti? - if data_dict['save_dicom'] is True: - # Save the volume as a DICOM file, with each TR as its own file - output_file = os.path.join(outputDir, 'rt_' + format(idx, '03d') + '.dcm') - write_dicom(output_file, brain_int32, idx+1) - else: - # Save the volume as a numpy file, with each TR as its own file - output_file = os.path.join(outputDir, 'rt_' + format(idx, '03d') + '.npy') - np.save(output_file, brain_int32) - - print("Generate {}".format(output_file)) - - # Sleep until next TR - if data_dict['save_realtime'] == 1: - time.sleep(data_dict['trDuration']) - - -if __name__ == '__main__': - # Receive the inputs - argParser = argparse.ArgumentParser( - 'Specify input arguments. Some arguments are parameters that require ' - 'an input is provided (noted by "Param"), others are flags that when ' - 'provided will change according to the flag (noted by "Flag")') - argParser.add_argument('--inputDir', '-i', default=None, type=str, - help='Param. Input directory for fmrisim parameters') - argParser.add_argument('--outputDir', '-o', default=None, type=str, - help='Param. Output directory for simulated data') - argParser.add_argument('--numTRs', '-n', default=200, type=int, - help='Param. Number of time points') - argParser.add_argument('--eventDuration', '-d', default=10, type=int, - help='Param. Number of seconds per event') - argParser.add_argument('--signalScale', '-s', default=0.5, type=float, - help='Param. Percent signal change') - argParser.add_argument('--useMultivariate', '-m', default=False, - action='store_true', - help='Flag. Signal is different between conditions ' - 'in a multivariate, versus univariate, way') - argParser.add_argument('--useDifferentROIs', '-r', default=False, - action='store_true', help='Flag. Use different ' - 'ROIs for each condition') - argParser.add_argument('--saveAsDicom', default=False, action='store_true', - help='Flag. Output files in DICOM format rather ' - 'than numpy') - argParser.add_argument('--saveRealtime', default=False, action='store_true', - help='Flag. Save data as if it was coming in at ' - 'the acquisition rate') - args = argParser.parse_args() - - inputDir = args.inputDir - outputDir = args.outputDir - - if inputDir is None or outputDir is None: - print("Must specify an input and output directory using -i and -o") - exit(-1) - - data_dict = {} - - ## User controlled settings - - # Specify the number of time points - data_dict['numTRs'] = args.numTRs - - # How long is each event/block you are modelling (assumes 6s rest between) - data_dict['event_duration'] = float(args.eventDuration) - - # What is the percent signal change being simulated - data_dict['scale_percentage'] = args.signalScale - - # Are there different ROIs for each condition (True) or is it in the same - # ROI (False). - # If it is the same ROI and you are using univariate differences, - # the second condition will have a smaller evoked response than the other. - data_dict['different_ROIs'] = args.useDifferentROIs - - # Is this a multivariate pattern (1) or a univariate pattern - data_dict['multivariate_pattern'] = args.useMultivariate - - # Do you want to save data as a dicom (True) or numpy (False) - data_dict['save_dicom'] = args.saveAsDicom - - # Do you want to save the data in real time (1) or as fast as possible (0)? - data_dict['save_realtime'] = args.saveRealtime - - ## Default settings - - # How long does each acquisition take - data_dict['trDuration'] = 2 - - # What is the time between each event (in seconds) - data_dict['isi'] = 6 - - # How long before the first event (in seconds) - data_dict['burn_in'] = 6 - - # Run the function if running from command line - generate_data(inputDir, - outputDir, - data_dict) diff --git a/examples/utils/sim_parameters/mask.npy b/examples/utils/sim_parameters/mask.npy deleted file mode 100644 index d3f1757fe..000000000 Binary files a/examples/utils/sim_parameters/mask.npy and /dev/null differ diff --git a/setup.py b/setup.py index 1833945e6..237e5f1ed 100644 --- a/setup.py +++ b/setup.py @@ -139,6 +139,7 @@ def finalize_options(self): 'nibabel', 'joblib', 'wheel', # See https://github.com/astropy/astropy-helpers/issues/501 + 'pydicom', ], author='Princeton Neuroscience Institute and Intel Corporation', author_email='mihai.capota@intel.com', @@ -150,7 +151,7 @@ def finalize_options(self): ext_modules=ext_modules, cmdclass={'build_ext': BuildExt}, packages=find_packages(), - package_data={'brainiak.utils': ['grey_matter_mask.npy']}, + include_package_data=True, python_requires='>=3.5', zip_safe=False, ) diff --git a/tests/utils/test_fmrisim_real_time.py b/tests/utils/test_fmrisim_real_time.py new file mode 100644 index 000000000..f348a2931 --- /dev/null +++ b/tests/utils/test_fmrisim_real_time.py @@ -0,0 +1,159 @@ +# Copyright 2016 Intel Corporation +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""fmrisim real-time generator test script + + Authors: Cameron Ellis (Princeton) 2020 +""" +import numpy as np +from brainiak.utils import fmrisim_real_time_generator as gen +import pytest +import os +import time +import glob +from pkg_resources import resource_stream +from typing import Dict +from nibabel.nifti1 import Nifti1Image +import gzip + +# Test that it crashes without inputs +with pytest.raises(TypeError): + gen.generate_data() # type: ignore + +data_dict = {} # type: Dict +vol = resource_stream(gen.__name__, "sim_parameters/ROI_A.nii.gz").read() +data_dict['ROI_A_file'] = Nifti1Image.from_bytes(gzip.decompress( + vol)).get_data() +vol = resource_stream(gen.__name__, "sim_parameters/ROI_B.nii.gz").read() +data_dict['ROI_B_file'] = Nifti1Image.from_bytes(gzip.decompress( + vol)).get_data() +vol = resource_stream(gen.__name__, + "sim_parameters/sub_template.nii.gz").read() +data_dict['template_path'] = Nifti1Image.from_bytes(gzip.decompress( + vol)).get_data() +noise_dict_file = resource_stream(gen.__name__, + "sim_parameters/sub_noise_dict.txt").read() +data_dict['noise_dict_file'] = noise_dict_file +data_dict['numTRs'] = 30 +data_dict['event_duration'] = 2 +data_dict['scale_percentage'] = 1 +data_dict['different_ROIs'] = True +data_dict['multivariate_pattern'] = False +data_dict['save_dicom'] = False +data_dict['save_realtime'] = False +data_dict['trDuration'] = 2 +data_dict['isi'] = 4 +data_dict['burn_in'] = 6 + + +# Run default test +def test_default(tmp_path, dd=data_dict): + + # Run the simulation + gen.generate_data(str(tmp_path), + dd) + + # Check that there are 32 files where there should be (30 plus label and + # mask) + assert len(os.listdir(str(tmp_path))) == 32, "Incorrect file number" + + # Check that the data is the right shape + input_template = dd['template_path'] + input_shape = input_template.shape + output_vol = np.load(tmp_path / 'rt_000.npy') + output_shape = output_vol.shape + assert input_shape == output_shape, 'Output shape is incorrect' + + # Check the labels have the correct count + labels = np.load(tmp_path / 'labels.npy') + + assert np.sum(labels > 0) == 9, 'Incorrect number of events' + + +def test_signal_size(tmp_path, dd=data_dict): + + # Change it to only use ROI A + dd['different_ROIs'] = False + + # Make the signal large + dd['scale_percentage'] = 100 + + # Run the simulation + gen.generate_data(str(tmp_path), + dd) + + # Load in the ROI masks + ROI_A = dd['ROI_A_file'] + ROI_B = dd['ROI_B_file'] + + # Load in the data just simulated + ROI_A_mean = [] + ROI_B_mean = [] + for TR_counter in range(dd['numTRs']): + + # Load the data + vol_name = 'rt_%03d.npy' % TR_counter + vol = np.load(tmp_path / vol_name) + + # Mask the data + ROI_A_mean += [np.mean(vol[ROI_A == 1])] + ROI_B_mean += [np.mean(vol[ROI_B == 1])] + + assert np.std(ROI_A_mean) > np.std(ROI_B_mean), 'Signal not scaling' + + +def test_multivariate(tmp_path, dd=data_dict): + + dd['multivariate_pattern'] = True + dd['different_ROIs'] = False + + # Make the signal large + dd['scale_percentage'] = 100 + + # Run the simulation + gen.generate_data(str(tmp_path), + dd) + + # Load in the ROI masks + ROI_A = dd['ROI_A_file'] + ROI_B = dd['ROI_B_file'] + + # Test this volume + vol = np.load(str(tmp_path / 'rt_007.npy')) + + ROI_A_std = np.std(vol[ROI_A == 1]) + ROI_B_std = np.std(vol[ROI_B == 1]) + + assert ROI_A_std > ROI_B_std, 'Multivariate not making variable signal' + + +def test_save_dicoms_realtime(tmp_path, dd=data_dict): + + dd['save_dicom'] = True + dd['save_realtime'] = True + + start_time = time.time() + + # Run the simulation + gen.generate_data(str(tmp_path), + dd) + + end_time = time.time() + + # Check it took 2s per TR + assert (end_time - start_time) > 60, 'Realtime ran fast' + + # Check correct file number + file_path = str(tmp_path / '*.dcm') + assert len(glob.glob(file_path)) == 30, "Wrong dicom file num"