diff --git a/.gitignore b/.gitignore index d34f46b..e5f0993 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,5 @@ **/*.pyc **/*.pyo *.egg-info/* +build/ + diff --git a/build/lib/rsHRF/CLI.py b/build/lib/rsHRF/CLI.py new file mode 100644 index 0000000..792af34 --- /dev/null +++ b/build/lib/rsHRF/CLI.py @@ -0,0 +1,447 @@ +import sys +import numpy as np +import os.path as op +import json +from argparse import ArgumentParser +from bids.layout import BIDSLayout +from pathlib import Path +from rsHRF import spm_dep, fourD_rsHRF, utils + +import warnings +warnings.filterwarnings("ignore") + +with open(op.join(op.dirname(op.realpath(__file__)), "VERSION"), "r") as fh: + __version__ = fh.read().strip('\n') + +def get_parser(): + parser = ArgumentParser(description='retrieves the onsets of pseudo-events triggering a ' + 'haemodynamic response from resting state fMRI BOLD ' + 'voxel-wise signal') + + group_input = parser.add_mutually_exclusive_group(required=True) + + group_input.add_argument('--ts', action='store', type=op.abspath, + help='the absolute path to a single data file') + + group_input.add_argument('--input_file', action='store', type=op.abspath, + help='the absolute path to a single data file') + + group_input.add_argument('--bids_dir', nargs='?', action='store', type=op.abspath, + help='the root folder of a BIDS valid dataset ' + '(sub-XXXXX folders should be found at the ' + 'top level in this folder).') + + group_input.add_argument('--GUI', action='store_true', + help='to execute the toolbox in GUI mode') + + parser.add_argument('--output_dir', action='store', type=op.abspath, + help='the output path for the outcomes of processing') + + parser.add_argument('--n_jobs', action='store', type=int, default=-1, + help='the number of parallel processing elements') + + parser.add_argument('-V', '--version', action='version', version='rsHRF version {}'.format(__version__)) + + parser.add_argument('--analysis_level', help='Level of the analysis that will be performed. ' + 'Multiple participant level analyses can be run independently ' + '(in parallel) using the same output_dir.', choices=['participant'], nargs='?') + + parser.add_argument('--participant_label', + help='The label(s) of the participant(s) that should be analyzed. The label ' + 'corresponds to sub- from the BIDS spec ' + '(so it does not include "sub-"). If this parameter is not ' + 'provided all subjects should be analyzed. Multiple ' + 'participants can be specified with a space separated list.', + nargs="+") + + parser.add_argument('--bids_filter_file', action='store', type=op.abspath, + help='a JSON file describing custom BIDS input filters using PyBIDS. ' + 'For further details, please check out http://bids-apps.neuroimaging.io/rsHRF/') + + group_mask = parser.add_mutually_exclusive_group(required=False) + + group_mask.add_argument('--atlas', action='store', type=op.abspath, + help='the absolute path to a single atlas file') + + group_mask.add_argument('--brainmask', action='store_true', + help='to enable the use of mask files present in the BIDS ' + 'directory itself') + + group_para = parser.add_argument_group('Parameters') + + group_para.add_argument('--estimation', action='store', + choices=['canon2dd', 'sFIR', 'FIR', 'fourier', 'hanning', 'gamma'], + help='Choose the estimation procedure from ' + 'canon2dd (canonical shape with 2 derivatives), ' + 'sFIR (smoothed Finite Impulse Response), ' + 'FIR (Finite Impulse Response), ' + 'fourier (Fourier Basis Set), ' + 'hanning (Fourier Basis w Hanning), ' + 'gamma (Gamma Basis Set)') + + group_para.add_argument('--passband', action='store', type=float, nargs=2, metavar=('LOW_FREQ','HIGH_FREQ'), + default=[0.01, 0.08], + help='set intervals for bandpass filter, default is 0.01 - 0.08') + + group_para.add_argument('--passband_deconvolve', action='store', type=float, nargs=2, metavar=('LOW_FREQ', 'HIGH_FREQ'), + default=[0.0, sys.float_info.max], + help='set intervals for bandpass filter (used while deconvolving BOLD), default is no-filtering') + + group_para.add_argument('-TR', action='store', type=float, default=-1, + help='set TR parameter') + + group_para.add_argument('-T', action='store', type=int, default=3, + help='set T parameter') + + group_para.add_argument('-T0', action='store', type=int, default=1, + help='set T0 parameter') + + group_para.add_argument('-TD_DD', action='store', type=int, default=2, + help='set TD_DD parameter') + + group_para.add_argument('-AR_lag', action='store', type=int, default=1, + help='set AR_lag parameter') + + group_para.add_argument('--thr', action='store', type=float, default=1, + help='set thr parameter') + + group_para.add_argument('--temporal_mask', action='store', type=op.abspath, + help='the path for the (temporal) mask file.\n The mask file should be a ".dat" file, consisting of a binary string of the same length as the signal') + + group_para.add_argument('--order', action='store', type=int, default=3, + help='set the number of basis vectors') + + group_para.add_argument('--len', action='store', type=int, default=24, + help='set len parameter') + + group_para.add_argument('--min_onset_search', action='store', type=int, default=4, + help='set min_onset_search parameter') + + group_para.add_argument('--max_onset_search', action='store', type=int, default=8, + help='set max_onset_search parameter') + + group_para.add_argument('--localK', action='store', type=int, + help='set localK') + + group_para.add_argument('--wiener', action='store_true', + help='to perform iterative wiener deconvolution') + + return parser + + +def run_rsHRF(): + parser = get_parser() + args = parser.parse_args() + arg_groups = {} + for group in parser._action_groups: + group_dict = {a.dest: getattr(args, a.dest, None) for a in group._group_actions } + arg_groups[group.title] = group_dict + para = arg_groups['Parameters'] + nargs = len(sys.argv) + temporal_mask = [] + + if (not args.GUI) and (args.output_dir is None): + parser.error('--output_dir is required when executing in command-line interface') + + if (not args.GUI) and (args.estimation is None): + parser.error('--estimation rule is required when executing in command-line interface') + + if (args.GUI): + if (nargs == 2): + try: + from .rsHRF_GUI import run + run.run() + except ModuleNotFoundError: + parser.error('--GUI should not be used inside a Docker container') + else: + parser.error('--no other arguments should be supplied with --GUI') + + if (args.input_file is not None or args.ts is not None) and args.analysis_level: + parser.error('analysis_level cannot be used with --input_file or --ts, do not supply it') + + if (args.input_file is not None or args.ts is not None) and args.participant_label: + parser.error('participant_labels are not to be used with --input_file or --ts, do not supply it') + + if args.input_file is not None and args.brainmask: + parser.error('--brainmask cannot be used with --input_file, use --atlas instead') + + if args.ts is not None and (args.brainmask or args.atlas): + parser.error('--atlas or --brainmask cannot be used with --ts, do not supply it') + + if args.bids_dir is not None and not (args.brainmask or args.atlas): + parser.error('--atlas or --brainmask needs to be supplied with --bids_dir') + + if args.bids_dir is not None and not args.analysis_level: + parser.error('analysis_level needs to be supplied with bids_dir, choices=[participant]') + + if args.input_file is not None and (not args.input_file.endswith(('.nii', '.nii.gz', '.gii', '.gii.gz'))): + parser.error('--input_file should end with .gii, .gii.gz, .nii or .nii.gz') + + if args.atlas is not None and (not args.atlas.endswith(('.nii', '.nii.gz','.gii', '.gii.gz'))): + parser.error('--atlas should end with .gii, .gii.gz, .nii or .nii.gz') + + if args.ts is not None and (not args.ts.endswith(('.txt'))): + parser.error('--ts file should end with .txt') + + if args.temporal_mask is not None and (not args.temporal_mask.endswith(('.dat'))): + parser.error('--temporal_mask ile should end with ".dat"') + + if args.temporal_mask is not None: + f = open(args.temporal_mask,'r') + for line in f: + for each in line: + if each in ['0','1']: + temporal_mask.append(int(each)) + + if args.estimation == 'sFIR' or args.estimation == 'FIR': + para['T'] = 1 + + if args.ts is not None: + file_type = op.splitext(args.ts) + if para['TR'] <= 0: + parser.error('Please supply a valid TR using -TR argument') + else: + TR = para['TR'] + para['dt'] = para['TR'] / para['T'] + para['lag'] = np.arange(np.fix(para['min_onset_search'] / para['dt']), + np.fix(para['max_onset_search'] / para['dt']) + 1, + dtype='int') + fourD_rsHRF.demo_rsHRF(args.ts, None, args.output_dir, para, args.n_jobs, file_type, mode='time-series', temporal_mask=temporal_mask, wiener=args.wiener) + + if args.input_file is not None: + if args.atlas is not None: + if (args.input_file.endswith(('.nii', '.nii.gz')) and args.atlas.endswith(('.gii', '.gii.gz'))) or (args.input_file.endswith(('.gii', '.gii.gz')) and args.atlas.endswith(('.nii', '.nii.gz'))): + parser.error('--atlas and input_file should be of the same type [NIfTI or GIfTI]') + + # carry analysis with input_file and atlas + file_type = op.splitext(args.input_file) + if file_type[-1] == ".gz": + file_type = op.splitext(file_type[-2])[-1] + file_type[-1] + else: + file_type = file_type[-1] + if ".nii" in file_type: + TR = (spm_dep.spm.spm_vol(args.input_file).header.get_zooms())[-1] + else: + if para['TR'] == -1: + parser.error('Please supply a valid TR using -TR argument') + else: + TR = para['TR'] + if TR <= 0: + if para['TR'] <= 0: + parser.error('Please supply a valid TR using -TR argument') + else: + if para['TR'] == -1: + para['TR'] = TR + elif para['TR'] <= 0: + print('Invalid TR supplied, using implicit TR: {0}'.format(TR)) + para['TR'] = TR + para['dt'] = para['TR'] / para['T'] + para['lag'] = np.arange(np.fix(para['min_onset_search'] / para['dt']), + np.fix(para['max_onset_search'] / para['dt']) + 1, + dtype='int') + fourD_rsHRF.demo_rsHRF(args.input_file, args.atlas, args.output_dir, para, args.n_jobs, file_type, mode='input', temporal_mask=temporal_mask, wiener=args.wiener) + + if args.bids_dir is not None: + utils.bids.write_derivative_description(args.bids_dir, args.output_dir) + bids_dir = Path(args.bids_dir) + fname = bids_dir / 'dataset_description.json' + + if fname.exists(): + desc = json.loads(Path(fname).read_text()) + if 'DataType' in desc : + if desc['DataType'] != 'derivative': + parser.error('Input data is not a derivative dataset' + ' (DataType in dataset_description.json is not equal to "derivative")') + + else : + parser.error('DataType is not defined in the dataset_description.json file. Please make sure DataType is defined. ' + 'Information on the dataset_description.json file can be found online ' + '(https://bids-specification.readthedocs.io/en/stable/03-modality-agnostic-files.html' + '#derived-dataset-and-pipeline-description)') + else : + parser.error('Could not find dataset_description.json file. Please make sure the BIDS data ' + 'structure is present and correct. Datasets can be validated online ' + 'using the BIDS Validator (http://incf.github.io/bids-validator/).') + + + if args.bids_dir is not None and args.atlas is not None: + # carry analysis with bids_dir and 1 atlas + layout = BIDSLayout(args.bids_dir, validate=False, config =['bids', 'derivatives']) + + if args.participant_label: + input_subjects = args.participant_label + subjects_to_analyze = layout.get_subjects(subject=input_subjects) + else: + subjects_to_analyze = layout.get_subjects() + + if not subjects_to_analyze: + parser.error('Could not find participants. Please make sure the BIDS data ' + 'structure is present and correct. Datasets can be validated online ' + 'using the BIDS Validator (http://incf.github.io/bids-validator/).') + + if not args.atlas.endswith(('.nii', '.nii.gz')): + parser.error('--atlas should end with .nii or .nii.gz') + + if args.bids_filter_file is not None: + filter_list = json.loads(Path(args.bids_filter_file).read_text()) + + default_input = {'extension': 'nii.gz', + 'datatype' : 'func', + 'desc': 'preproc', + 'task' : 'rest', + 'suffix': 'bold'} + default_input['subject']=subjects_to_analyze + default_input.update(filter_list['bold']) + + all_inputs = layout.get(return_type='filename',**default_input) + + else : + all_inputs = layout.get(return_type='filename',datatype='func', subject=subjects_to_analyze, task='rest',desc='preproc',suffix='bold', extension=['nii', 'nii.gz']) + + if not all_inputs != []: + parser.error('There are no files of type *bold.nii / *bold.nii.gz ' + 'Please make sure to have at least one file of the above type ' + 'in the BIDS specification') + else: + num_errors = 0 + for file_count in range(len(all_inputs)): + try: + TR = layout.get_metadata(all_inputs[file_count])['RepetitionTime'] + except KeyError as e: + TR = spm_dep.spm.spm_vol(all_inputs[file_count]).header.get_zooms()[-1] + para['TR'] = TR + para['dt'] = para['TR'] / para['T'] + para['lag'] = np.arange(np.fix(para['min_onset_search'] / para['dt']), + np.fix(para['max_onset_search'] / para['dt']) + 1, + dtype='int') + num_errors += 1 + try: + fourD_rsHRF.demo_rsHRF(all_inputs[file_count], args.atlas, args.output_dir, para, args.n_jobs, file_type, mode='bids w/ atlas', temporal_mask=temporal_mask, wiener=args.wiener) + num_errors -=1 + except ValueError as err: + print(err.args[0]) + except: + print("Unexpected error:", sys.exc_info()[0]) + success = len(all_inputs) - num_errors + if success == 0: + raise RuntimeError('Dimensions were inconsistent for all input-mask pairs; \n' + 'No inputs were processed!') + + if args.bids_dir is not None and args.brainmask: + # carry analysis with bids_dir and brainmask + layout = BIDSLayout(args.bids_dir, validate=False, config =['bids', 'derivatives']) + + if args.participant_label: + input_subjects = args.participant_label + subjects_to_analyze = layout.get_subjects(subject=input_subjects) + else: + subjects_to_analyze = layout.get_subjects() + + if not subjects_to_analyze: + parser.error('Could not find participants. Please make sure the BIDS data ' + 'structure is present and correct. Datasets can be validated online ' + 'using the BIDS Validator (http://incf.github.io/bids-validator/).') + + if args.bids_filter_file is not None: + filter_list = json.loads(Path(args.bids_filter_file).read_text()) + + default_input = {'extension': 'nii.gz', + 'datatype' : 'func', + 'desc': 'preproc', + 'task' : 'rest', + 'suffix': 'bold'} + default_input['subject']=subjects_to_analyze + default_input.update(filter_list['bold']) + + all_inputs = layout.get(return_type='filename',**default_input) + + default_mask={'extension': 'nii.gz', + 'datatype': 'func', + 'desc': 'brain', + 'task':'rest', + 'suffix':'mask'} + default_mask['subject']=subjects_to_analyze + default_mask.update(filter_list['mask']) + + all_masks = layout.get(return_type='filename',**default_mask) + + + else: + all_inputs = layout.get(return_type='filename',datatype='func', subject=subjects_to_analyze, task='rest',desc='preproc',suffix='bold', extension=['nii', 'nii.gz']) + all_masks = layout.get(return_type='filename', datatype='func', subject=subjects_to_analyze, task='rest',desc='brain',suffix='mask', extension=['nii', 'nii.gz']) + + if not all_inputs != []: + parser.error('There are no files of type *bold.nii / *bold.nii.gz ' + 'Please make sure to have at least one file of the above type ' + 'in the BIDS specification') + if not all_masks != []: + parser.error('There are no files of type *mask.nii / *mask.nii.gz ' + 'Please make sure to have at least one file of the above type ' + 'in the BIDS specification') + if len(all_inputs) != len(all_masks): + parser.error('The number of *bold.nii / .nii.gz and the number of ' + '*mask.nii / .nii.gz are different. Please make sure that ' + 'there is one mask for each input_file present') + + all_inputs.sort() + all_masks.sort() + + all_prefix_match = False + prefix_match_count = 0 + for i in range(len(all_inputs)): + input_prefix = all_inputs[i].split('/')[-1].split('_desc')[0] + mask_prefix = all_masks[i].split('/')[-1].split('_desc')[0] + if input_prefix == mask_prefix: + prefix_match_count += 1 + else: + all_prefix_match = False + break + if prefix_match_count == len(all_inputs): + all_prefix_match = True + + if not all_prefix_match: + parser.error('The mask and input files should have the same prefix for correspondence. ' + 'Please consider renaming your files') + else: + num_errors = 0 + for file_count in range(len(all_inputs)): + file_type = all_inputs[file_count].split('bold')[1] + if file_type == ".nii" or file_type == ".nii.gz": + try: + TR = layout.get_metadata(all_inputs[file_count])['RepetitionTime'] + except KeyError as e: + TR = spm_dep.spm.spm_vol(all_inputs[file_count]).header.get_zooms()[-1] + para['TR'] = TR + else: + spm_dep.spm.spm_vol(all_inputs[file_count]) + TR = spm_dep.spm.spm_vol(all_inputs[file_count]).get_arrays_from_intent("NIFTI_INTENT_TIME_SERIES")[0].meta.get_metadata()["TimeStep"] + para['TR'] = float(TR) * 0.001 + + + para['dt'] = para['TR'] / para['T'] + para['lag'] = np.arange(np.fix(para['min_onset_search'] / para['dt']), + np.fix(para['max_onset_search'] / para['dt']) + 1, + dtype='int') + num_errors += 1 + try: + fourD_rsHRF.demo_rsHRF(all_inputs[file_count], all_masks[file_count], args.output_dir, para, args.n_jobs, mode='bids', temporal_mask=temporal_mask, wiener=args.wiener) + num_errors -=1 + except ValueError as err: + print(err.args[0]) + except: + print("Unexpected error:", sys.exc_info()[0]) + success = len(all_inputs) - num_errors + if success == 0: + raise RuntimeError('Dimensions were inconsistent for all input-mask pairs; \n' + 'No inputs were processed!') + + + +def main(): + warnings.filterwarnings("ignore") + run_rsHRF() + + +if __name__ == '__main__': + raise RuntimeError("CLI.py should not be run directly;\n" + "Please `pip install` rsHRF and use the `rsHRF` command") \ No newline at end of file diff --git a/build/lib/rsHRF/VERSION b/build/lib/rsHRF/VERSION new file mode 100755 index 0000000..fa5512a --- /dev/null +++ b/build/lib/rsHRF/VERSION @@ -0,0 +1 @@ +1.5.8 \ No newline at end of file diff --git a/build/lib/rsHRF/__about__.py b/build/lib/rsHRF/__about__.py new file mode 100644 index 0000000..33bc2ee --- /dev/null +++ b/build/lib/rsHRF/__about__.py @@ -0,0 +1,12 @@ +"""Base module variables.""" +from os import path as op + +with open(op.join(op.dirname(op.realpath(__file__)), "VERSION"), "r") as fh: + __version__ = fh.read().strip('\n') + +__packagename__ = 'rsHRF' +__url__ = 'https://github.com/BIDSapps/rsHRF' + +DOWNLOAD_URL = ( + 'https://github.com/BIDS-Apps/{name}/'.format( + name=__packagename__)) \ No newline at end of file diff --git a/build/lib/rsHRF/__init__.py b/build/lib/rsHRF/__init__.py new file mode 100644 index 0000000..1b18cf2 --- /dev/null +++ b/build/lib/rsHRF/__init__.py @@ -0,0 +1,10 @@ +import rsHRF.spm_dep +import rsHRF.processing +import rsHRF.canon +import rsHRF.sFIR +import rsHRF.parameters +import rsHRF.fourD_rsHRF +import rsHRF.CLI +__all__ = ["spm_dep", "processing", "canon", "utils", + "sFIR", "parameters", "basis_functions", + "fourD_rsHRF.py", "CLI.py", "iterative_wiener_deconv"] diff --git a/build/lib/rsHRF/__main__.py b/build/lib/rsHRF/__main__.py new file mode 100644 index 0000000..a1d7bba --- /dev/null +++ b/build/lib/rsHRF/__main__.py @@ -0,0 +1,3 @@ +from rsHRF import CLI + +CLI.main() diff --git a/build/lib/rsHRF/basis_functions/__init__.py b/build/lib/rsHRF/basis_functions/__init__.py new file mode 100644 index 0000000..b3b472e --- /dev/null +++ b/build/lib/rsHRF/basis_functions/__init__.py @@ -0,0 +1 @@ +from . import basis_functions diff --git a/build/lib/rsHRF/basis_functions/basis_functions.py b/build/lib/rsHRF/basis_functions/basis_functions.py new file mode 100644 index 0000000..d0d50b2 --- /dev/null +++ b/build/lib/rsHRF/basis_functions/basis_functions.py @@ -0,0 +1,69 @@ +import math +import numpy as np +from scipy.stats import gamma +from rsHRF import canon, spm_dep, sFIR + +import warnings +warnings.filterwarnings("ignore") + +""" +BASIS FUNCTION COMPUTATION +""" + +def get_basis_function(bold_sig_shape, para): + N, nvar = bold_sig_shape + dt = para['dt'] # 'time bin for basis functions {secs}'; + l = para['len'] + if "gamma" in para['estimation']: + pst = np.arange(0,l+0.01,dt) #the 0.01 is because the rounding is different between python and matlab + bf = gamma_bf(pst,para['order']) + elif 'fourier' in para['estimation'] or 'hanning' in para['estimation']: + pst = np.arange(0,l+0.01,dt) #the 0.01 is because the rounding is different between python and matlab + pst = pst/max(pst) + bf = fourier_bf(pst,para) + elif 'canon' in para['estimation']: + bf = canon2dd_bf(bold_sig_shape, para) + bf = spm_dep.spm.spm_orth(np.asarray(bf)) + return bf + +def canon2dd_bf(data_shape, para): + """ + Returns canon basis functions + """ + N, nvar = data_shape + bf = canon.canon_hrf2dd.wgr_spm_get_canonhrf(para) + if 'Volterra' in para: + if para['Volterra'] == 2: + bf2 = np.einsum('i...,j...->ij...',bf.T,bf.T).reshape(-1, bf.shape[0]).T + bf = np.column_stack((bf, bf2)) + return bf + +def fourier_bf(pst, para): + """ + Returns Fourier (Hanning) basis functions + """ + if "hanning" in para['estimation']: + g = (1 - np.cos(2*math.pi*pst)) / 2 + else: + g = np.ones(len(pst)) + sin_ = lambda x : np.sin(x*2*math.pi) + cos_ = lambda x : np.cos(x*2*math.pi) + sin_ = np.vectorize(sin_) + cos_ = np.vectorize(cos_) + arr = np.arange(1, para['order']+1) + s = sin_(np.einsum('i,j->ij',arr,pst)) + c = cos_(np.einsum('i,j->ij',arr,pst)) + s = np.multiply(g,s) + c = np.multiply(g,c) + g = np.expand_dims(g, axis=1).T + return np.concatenate((g, s, c), axis=0).T + +def gamma_bf(u,h): + """ + Returns Gamma basis functions + """ + arr = np.arange(2, h+2) + f = np.vectorize(gamma.pdf, signature='(n),()->(n)') + m = np.power(2, arr) + return f(u,m).T + diff --git a/build/lib/rsHRF/canon/__init__.py b/build/lib/rsHRF/canon/__init__.py new file mode 100644 index 0000000..1674993 --- /dev/null +++ b/build/lib/rsHRF/canon/__init__.py @@ -0,0 +1 @@ +from . import canon_hrf2dd diff --git a/build/lib/rsHRF/canon/canon_hrf2dd.py b/build/lib/rsHRF/canon/canon_hrf2dd.py new file mode 100644 index 0000000..78826cf --- /dev/null +++ b/build/lib/rsHRF/canon/canon_hrf2dd.py @@ -0,0 +1,29 @@ +import numpy as np +from ..spm_dep import spm + +import warnings +warnings.filterwarnings("ignore") + +def wgr_spm_get_canonhrf(xBF): + dt = xBF['dt'] + fMRI_T = xBF['T'] + p = np.array([6, 16, 1, 1, 6, 0, 32], dtype=float) + p[len(p) - 1] = xBF['len'] + bf = spm.spm_hrf(dt, p, fMRI_T) + bf = bf[:, np.newaxis] + # time-derivative + if xBF['TD_DD']: + dp = 1 + p[5] = p[5] + dp + D = (bf[:, 0] - spm.spm_hrf(dt, p, fMRI_T)) / dp + D = D[:, np.newaxis] + bf = np.append(bf, D, axis=1) + p[5] = p[5] - dp + # dispersion-derivative + if xBF['TD_DD'] == 2: + dp = 0.01 + p[2] = p[2] + dp + D = (bf[:, 0] - spm.spm_hrf(dt, p, fMRI_T)) / dp + D = D[:, np.newaxis] + bf = np.append(bf, D, axis=1) + return bf diff --git a/build/lib/rsHRF/fourD_rsHRF.py b/build/lib/rsHRF/fourD_rsHRF.py new file mode 100644 index 0000000..d178f04 --- /dev/null +++ b/build/lib/rsHRF/fourD_rsHRF.py @@ -0,0 +1,201 @@ +import os +import matplotlib +matplotlib.use('agg') +import numpy as np +import nibabel as nib +import scipy.io as sio +import matplotlib.pyplot as plt +from bids.layout import BIDSLayout, parse_file_entities +from scipy import stats, signal +from scipy.sparse import lil_matrix +from rsHRF import spm_dep, processing, parameters, basis_functions, utils, iterative_wiener_deconv + +import warnings +warnings.filterwarnings("ignore") + +def demo_rsHRF(input_file, mask_file, output_dir, para, p_jobs, file_type=".nii", mode="bids", wiener=False, temporal_mask=[]): + # book-keeping w.r.t parameter values + if 'localK' not in para or para['localK'] == None: + if para['TR']<=2: + para['localK'] = 1 + else: + para['localK'] = 2 + # creating the output-directory if not already present + if not os.path.isdir(output_dir): + os.mkdir(output_dir) + # for four-dimensional input + if mode != 'time-series': + if mode == 'bids' or mode == 'bids w/ atlas': + name = input_file.split('/')[-1].split('.')[0] + v1 = spm_dep.spm.spm_vol(input_file) + else: + name = input_file.split('/')[-1].split('.')[0] + v1 = spm_dep.spm.spm_vol(input_file) + if mask_file != None: + if mode == 'bids': + mask_name = mask_file.split('/')[-1].split('.')[0] + v = spm_dep.spm.spm_vol(mask_file) + else: + mask_name = mask_file.split('/')[-1].split('.')[0] + v = spm_dep.spm.spm_vol(mask_file) + if file_type == ".nii" or file_type == ".nii.gz": + brain = spm_dep.spm.spm_read_vols(v) + else: + brain = v.agg_data().flatten(order='F') + if ((file_type == ".nii" or file_type == ".nii.gz") and \ + v1.header.get_data_shape()[:-1] != v.header.get_data_shape()) or \ + ((file_type == ".gii" or file_type == ".gii.gz") and \ + v1.agg_data().shape[0]!= v.agg_data().shape[0]): + raise ValueError ('Inconsistency in input-mask dimensions' + '\n\tinput_file == ' + name + file_type + '\n\tmask_file == ' + mask_name + file_type) + else: + if file_type == ".nii" or file_type == ".nii.gz" : + data = v1.get_fdata() + else: + data = v1.agg_data() + else: + print('No atlas provided! Generating mask file...') + if file_type == ".nii" or file_type == ".nii.gz" : + data = v1.get_fdata() + brain = np.nanvar(data.reshape(-1, data.shape[3]), -1, ddof=0) + else: + data = v1.agg_data() + brain = np.nanvar(data, -1, ddof=0) + print('Done') + voxel_ind = np.where(brain > 0)[0] + mask_shape = data.shape[:-1] + nobs = data.shape[-1] + data1 = np.reshape(data, (-1, nobs), order='F').T + bold_sig = stats.zscore(data1[:, voxel_ind], ddof=1) + # for time-series input + else: + name = input_file.split('/')[-1].split('.')[0] + data1 = (np.loadtxt(input_file, delimiter=",")) + if data1.ndim == 1: + data1 = np.expand_dims(data1, axis=1) + nobs = data1.shape[0] + bold_sig = stats.zscore(data1, ddof=1) + if len(temporal_mask) > 0 and len(temporal_mask) != nobs: + raise ValueError ('Inconsistency in temporal_mask dimensions.\n' + 'Size of mask: ' + str(len(temporal_mask)) + '\n' + 'Size of time-series: ' + str(nobs)) + bold_sig = np.nan_to_num(bold_sig) + bold_sig_deconv = processing. \ + rest_filter. \ + rest_IdealFilter(bold_sig, para['TR'], para['passband_deconvolve']) + bold_sig = processing. \ + rest_filter. \ + rest_IdealFilter(bold_sig, para['TR'], para['passband']) + data_deconv = np.zeros(bold_sig.shape) + event_number = np.zeros((1, bold_sig.shape[1])) + print('Retrieving HRF ...') + #Estimate HRF for the fourier / hanning / gamma / cannon basis functions + if not (para['estimation'] == 'sFIR' or para['estimation'] == 'FIR'): + bf = basis_functions.basis_functions.get_basis_function(bold_sig.shape, para) + beta_hrf, event_bold = utils.hrf_estimation.compute_hrf(bold_sig, para, temporal_mask, p_jobs, bf=bf) + hrfa = np.dot(bf, beta_hrf[np.arange(0, bf.shape[1]), :]) + #Estimate HRF for FIR and sFIR + else: + para['T'] = 1 + beta_hrf, event_bold = utils.hrf_estimation.compute_hrf(bold_sig, para, temporal_mask, p_jobs) + hrfa = beta_hrf[:-1,:] + nvar = hrfa.shape[1] + PARA = np.zeros((3, nvar)) + for voxel_id in range(nvar): + hrf1 = hrfa[:, voxel_id] + PARA[:, voxel_id] = \ + parameters.wgr_get_parameters(hrf1, para['TR'] / para['T']) + print('Done') + print('Deconvolving HRF ...') + if para['T'] > 1: + hrfa_TR = signal.resample_poly(hrfa, 1, para['T']) + else: + hrfa_TR = hrfa + for voxel_id in range(nvar): + hrf = hrfa_TR[:, voxel_id] + if not wiener: + H = np.fft.fft( + np.append(hrf, + np.zeros((nobs - max(hrf.shape), 1))), axis=0) + M = np.fft.fft(bold_sig_deconv[:, voxel_id]) + data_deconv[:, voxel_id] = \ + np.fft.ifft(H.conj() * M / (H * H.conj() + .1*np.mean((H * H.conj())))) + else: + data_deconv[:, voxel_id] = iterative_wiener_deconv.rsHRF_iterative_wiener_deconv(bold_sig_deconv[:, voxel_id], hrf) + event_number[:, voxel_id] = np.amax(event_bold[voxel_id].shape) + print('Done') + print('Saving Output ...') + # setting the output-path + if mode == 'bids' or mode == 'bids w/ atlas': + layout_output = BIDSLayout(output_dir) + entities = parse_file_entities(input_file) + sub_save_dir = layout_output.build_path(entities).rsplit('/',1)[0] + else: + sub_save_dir = output_dir + if not os.path.isdir(sub_save_dir): + os.makedirs(sub_save_dir, exist_ok=True) + dic = {'para': para, 'hrfa': hrfa, 'event_bold': event_bold, 'PARA': PARA} + ext = '_hrf.mat' + if mode == "time-series": + dic["event_number"] = event_number + dic["data_deconv"] = data_deconv + ext = '_hrf_deconv.mat' + name = name.rsplit('_bold', 1)[0] + sio.savemat(os.path.join(sub_save_dir, name + ext), dic) + HRF_para_str = ['height', 'T2P', 'FWHM'] + if mode != "time-series": + mask_data = np.zeros(mask_shape).flatten(order='F') + for i in range(3): + fname = os.path.join(sub_save_dir, + name + '_' + HRF_para_str[i]) + mask_data[voxel_ind] = PARA[i, :] + mask_data = mask_data.reshape(mask_shape, order='F') + spm_dep.spm.spm_write_vol(v1, mask_data, fname, file_type) + mask_data = mask_data.flatten(order='F') + fname = os.path.join(sub_save_dir, name + '_eventnumber') + mask_data[voxel_ind] = event_number + mask_data = mask_data.reshape(mask_shape, order='F') + spm_dep.spm.spm_write_vol(v1, mask_data, fname, file_type) + mask_data = np.zeros(data.shape) + dat3 = np.zeros(data.shape[:-1]).flatten(order='F') + for i in range(nobs): + fname = os.path.join(sub_save_dir, name + '_deconv') + dat3[voxel_ind] = data_deconv[i, :] + dat3 = dat3.reshape(data.shape[:-1], order='F') + if file_type == ".nii" or file_type == ".nii.gz" : + mask_data[:, :, :, i] = dat3 + else: + mask_data[:, i] = dat3 + dat3 = dat3.flatten(order='F') + spm_dep.spm.spm_write_vol(v1, mask_data, fname, file_type) + pos = 0 + while pos < hrfa_TR.shape[1]: + if np.any(hrfa_TR[:,pos]): + break + pos += 1 + event_plot = lil_matrix((1, nobs)) + if event_bold.size: + event_plot[:, event_bold[pos]] = 1 + else: + print("No Events Detected!") + return 0 + event_plot = np.ravel(event_plot.toarray()) + plt.figure() + plt.plot(para['TR'] * np.arange(1, np.amax(hrfa_TR[:, pos].shape) + 1), + hrfa_TR[:, pos], linewidth=1) + plt.xlabel('time (s)') + plt.savefig(os.path.join(sub_save_dir, name + '_hrf_plot.png')) + plt.figure() + plt.plot(para['TR'] * np.arange(1, nobs + 1), + np.nan_to_num(stats.zscore(bold_sig[:, pos], ddof=1)), + linewidth=1) + plt.plot(para['TR'] * np.arange(1, nobs + 1), + np.nan_to_num(stats.zscore(data_deconv[:, pos], ddof=1)), + color='r', linewidth=1) + markerline, stemlines, baseline = \ + plt.stem(para['TR'] * np.arange(1, nobs + 1), event_plot) + plt.setp(baseline, 'color', 'k', 'markersize', 1) + plt.setp(stemlines, 'color', 'k') + plt.setp(markerline, 'color', 'k', 'markersize', 3, 'marker', 'd') + plt.legend(['BOLD', 'Deconvolved BOLD', 'Events'], loc='best') + plt.xlabel('time (s)') + plt.savefig(os.path.join(sub_save_dir, name + '_deconvolution_plot.png')) + print('Done') + return 0 diff --git a/build/lib/rsHRF/iterative_wiener_deconv.py b/build/lib/rsHRF/iterative_wiener_deconv.py new file mode 100644 index 0000000..774ed7a --- /dev/null +++ b/build/lib/rsHRF/iterative_wiener_deconv.py @@ -0,0 +1,39 @@ +# import pyyawt +import pywt +import numpy as np +from rsHRF.processing import knee + +def rsHRF_iterative_wiener_deconv(y, h, Iterations=1000): + N = y.shape[0] + nh = max(h.shape) + h = np.append(h, np.zeros((N - nh, 1))) + H = np.fft.fft(h, axis=0) + Y = np.fft.fft(y, axis=0) + coeffs = pywt.wavedec(abs(y), 'db2', level=1) + sigma = np.median(np.abs(coeffs[-1])) / 0.6745 + Phh = np.square(abs(H)) + sqrdtempnorm = ((((np.linalg.norm(y-np.mean(y), 2)**2) - (N-1)*(sigma**2))) / (np.linalg.norm(h,1))**2) + Nf = (sigma**2)*N + tempreg = Nf/sqrdtempnorm + Pxx0 = np.square(abs(np.multiply(Y, (np.divide(np.conj(H), (Phh + N*tempreg)))))) + Pxx = Pxx0 + Sf = Pxx.reshape(-1, 1) + for i in range(0, Iterations): + M = np.divide(np.multiply(np.multiply(np.conjugate(H), Pxx), Y), np.add(np.multiply(np.square(abs(H)), Pxx), Nf)) + PxxY = np.divide(np.multiply(Pxx, Nf), np.add(np.multiply(np.square(abs(H)), Pxx), Nf)) + Pxx = np.add(PxxY, np.square(abs(M))) + Sf = np.concatenate((Sf, Pxx.reshape(-1, 1)), axis=1) + dSf = np.diff(Sf, 1, 1) + dSfmse = np.mean(np.square(dSf), axis=1) + _, idx = knee.knee_pt(dSfmse) + idm = np.argmin(dSfmse) + ratio = np.abs(dSfmse[idx] - dSfmse[idm])/(np.abs(np.max(dSfmse) - np.min(dSfmse))) + if ratio > 0.5: + id0 = idm + else: + id0 = idx + # Safe indexing to avoid IndexError + safe_index = min(id0 + 1, Sf.shape[1] - 1) + Pxx = Sf[:, safe_index] + WienerFilterEst = np.divide(np.multiply(np.conj(H), Pxx), np.add(np.multiply(np.square(abs(H)), Pxx), Nf)) + return np.real(np.fft.ifft(np.multiply(WienerFilterEst, Y))) \ No newline at end of file diff --git a/build/lib/rsHRF/parameters.py b/build/lib/rsHRF/parameters.py new file mode 100644 index 0000000..df81114 --- /dev/null +++ b/build/lib/rsHRF/parameters.py @@ -0,0 +1,44 @@ +import numpy as np +import warnings + +warnings.filterwarnings("ignore") + +def wgr_get_parameters(hdrf, dt): + """ + Find Model Parameters + h - Height + p - Time to peak (in units of dt where dt = TR/para.T) + w - Width at half-peak + """ + param = np.zeros((3, 1)) + + if(np.any(hdrf)): + n = np.fix(np.amax(hdrf.shape) * 0.8) + + p = np.argmax(np.absolute(hdrf[np.arange(0, n, dtype='int')])) + h = hdrf[p] + + if h > 0: + v = (hdrf >= (h / 2)) + else: + v = (hdrf <= (h / 2)) + v = v.astype(int) + b = np.argmin(np.diff(v)) + v[b + 1:] = 0 + w = np.sum(v) + + cnt = p - 1 + g = hdrf[1:] - hdrf[0:-1] + + while cnt > 0 and np.abs(g[cnt]) < 0.001: + h = hdrf[cnt - 1] + p = cnt + cnt = cnt - 1 + + param[0] = h + param[1] = (p + 1) * dt + param[2] = w * dt + + else: + print('.') + return param.ravel() diff --git a/build/lib/rsHRF/processing/__init__.py b/build/lib/rsHRF/processing/__init__.py new file mode 100644 index 0000000..9751036 --- /dev/null +++ b/build/lib/rsHRF/processing/__init__.py @@ -0,0 +1,2 @@ +from . import knee +from . import rest_filter diff --git a/build/lib/rsHRF/processing/knee.py b/build/lib/rsHRF/processing/knee.py new file mode 100644 index 0000000..d1cf95d --- /dev/null +++ b/build/lib/rsHRF/processing/knee.py @@ -0,0 +1,90 @@ +import numpy as np +import warnings + +warnings.filterwarnings("ignore") + +def knee_pt(y): + res_x, _id = knee_pt_helper(y) + _idm = np.argmin(y) + ratio = np.abs(y[_id]-y[_idm])/np.abs(np.max(y) - np.min(y)) + if ratio > 0.5: + idx_of_result = _idm + else: + idx_of_result = _id + return res_x, idx_of_result + +def knee_pt_helper(y, x=None): + x_was_none = False + use_absolute_dev_p = True + res_x = np.nan + idx_of_result = np.nan + + if not isinstance(y, np.ndarray): + print('knee_pt: y must be a numpy 1D vector') + return res_x, idx_of_result + else: + if y.ndim >= 2: + print('knee_pt: y must be 1 dimensional') + return res_x, idx_of_result + if np.size(y) == 0: + print('knee_pt: y can not be an empty vector') + return res_x, idx_of_result + else: + if x is None: + x_was_none = True + x = np.arange(1, np.amax(y.shape) + 1, dtype=int) + if x.shape != y.shape: + print('knee_pt: y and x must have the same dimensions') + return res_x, idx_of_result + if y.size < 3: + res_x, idx_of_result = np.min(y), np.argmin(y) + return res_x, idx_of_result + if np.all(np.diff(x) >= 0) and (not x_was_none): + idx = np.argsort(x) + y = np.sort(y) + x = np.sort(x) + else: + idx = np.arange(0, np.amax(x.shape)) + sigma_xy = np.cumsum(np.multiply(x, y), axis=0) + sigma_x = np.cumsum(x, axis=0) + sigma_y = np.cumsum(y, axis=0) + sigma_xx = np.cumsum(np.multiply(x, x), axis=0) + n = np.arange(1, np.amax(y.shape) + 1).conj().T + det = np.multiply(n, sigma_xx) - np.multiply(sigma_x, sigma_x) + mfwd = (np.multiply(n, sigma_xy) - + np.multiply(sigma_x, sigma_y)) / det + bfwd = -1 * ((np.multiply(sigma_x, sigma_xy) - + np.multiply(sigma_xx, sigma_y)) / det) + + sigma_xy = np.cumsum(np.multiply(x[::-1], y[::-1]), axis=0) + sigma_x = np.cumsum(x[::-1], axis=0) + sigma_y = np.cumsum(y[::-1], axis=0) + sigma_xx = np.cumsum(np.multiply(x[::-1], x[::-1]), axis=0) + n = np.arange(1, np.amax(y.shape) + 1).conj().T + det = np.multiply(n, sigma_xx) - np.multiply(sigma_x, sigma_x) + mbck = ((np.multiply(n, sigma_xy) - + np.multiply(sigma_x, sigma_y)) / det)[::-1] + bbck = (-1 * + ((np.multiply(sigma_x, sigma_xy) - + np.multiply(sigma_xx, sigma_y)) / det))[::-1] + + error_curve = np.full(y.shape, np.nan) + for breakpt in range(1, np.amax((y - 1).shape)): + delsfwd = (np.multiply(mfwd[breakpt], x[0:breakpt + 1]) + + bfwd[breakpt]) - y[0:breakpt + 1] + delsbck = (np.multiply(mbck[breakpt], x[breakpt:]) + + bbck[breakpt]) - y[breakpt:] + if use_absolute_dev_p: + error_curve[breakpt] = \ + np.sum(np.abs(delsfwd)) + np.sum(np.abs(delsbck)) + else: + error_curve[breakpt] = \ + np.sqrt(np.sum(np.multiply(delsfwd, delsfwd))) + \ + np.sqrt(np.sum(np.multiply(delsbck, delsbck))) + try: + loc = np.nanargmin(error_curve) + except ValueError as e: + loc = 0 + res_x = x[loc] + idx_of_result = idx[loc] + return res_x, idx_of_result diff --git a/build/lib/rsHRF/processing/rest_filter.py b/build/lib/rsHRF/processing/rest_filter.py new file mode 100644 index 0000000..9cb1e1b --- /dev/null +++ b/build/lib/rsHRF/processing/rest_filter.py @@ -0,0 +1,34 @@ +import numpy as np +import numpy.matlib +import warnings + +warnings.filterwarnings("ignore") + +def rest_IdealFilter(x, TR, Bands, m=5000): + nvar = x.shape[1] + nbin = int(np.ceil(nvar/m)) + for i in range(1, nbin + 1): + if i != nbin: + ind_X = [j for j in range((i-1)*m, i*m)] + else: + ind_X = [j for j in range((i-1)*m, nvar)] + x1 = x[:, ind_X] + x1 = conn_filter(TR,Bands,x1) + np.mean(x1, keepdims=True).repeat(x1.shape[0], axis=0) + x[:,ind_X] = x1 + return x + +def conn_filter(rt, filter, x): + Nx = x.shape[0] + fy = np.fft.fft(np.concatenate((x, np.flipud(x)), axis=0), axis=0) + f = np.arange(fy.shape[0]) + f = f.reshape(1, -1) + f = np.min((f, fy.shape[0]-f), axis=0) + low = filter[0]*(rt*fy.shape[0]) + high = filter[1]*(rt*fy.shape[0]) + idx_low = np.argwhere(np.any(f < low, axis=0)) + idx_high = np.argwhere(np.any(f >= high, axis=0)) + idx = np.concatenate((idx_low, idx_high)).reshape(-1) + fy[idx,:] = 0. + y = np.real(np.fft.ifft(fy, axis=0)) + y = y[0:Nx,:] + return y \ No newline at end of file diff --git a/build/lib/rsHRF/rsHRF_GUI/__init__.py b/build/lib/rsHRF/rsHRF_GUI/__init__.py new file mode 100644 index 0000000..d8668c3 --- /dev/null +++ b/build/lib/rsHRF/rsHRF_GUI/__init__.py @@ -0,0 +1,6 @@ +from . import core +from . import misc +from . import datatypes +from . import gui_windows + +__all__ = ["core", "misc", "datatypes", "gui_windows"] diff --git a/build/lib/rsHRF/rsHRF_GUI/core/__init__.py b/build/lib/rsHRF/rsHRF_GUI/core/__init__.py new file mode 100644 index 0000000..a8d584a --- /dev/null +++ b/build/lib/rsHRF/rsHRF_GUI/core/__init__.py @@ -0,0 +1 @@ +from .core import Core \ No newline at end of file diff --git a/build/lib/rsHRF/rsHRF_GUI/core/core.py b/build/lib/rsHRF/rsHRF_GUI/core/core.py new file mode 100644 index 0000000..ad48722 --- /dev/null +++ b/build/lib/rsHRF/rsHRF_GUI/core/core.py @@ -0,0 +1,356 @@ +import os +import numpy as np +import nibabel as nib +import scipy.io as sio +from bids.layout import BIDSLayout +from scipy import stats, signal +from scipy.sparse import lil_matrix +from ... import spm_dep, processing, parameters, basis_functions, utils + +from ..datatypes.timeseries.hrf import HRF +from ..datatypes.timeseries.bold_raw import BOLD_Raw +from ..datatypes.timeseries.bold_preprocessed import BOLD_Preprocessed +from ..datatypes.timeseries.bold_deconv import Bold_Deconv +from ..datatypes.misc.parameters import Parameters +from ..datatypes.misc.subject import Subject +from ..datatypes.misc.store import Store +from ..misc.status import Status + +import warnings +warnings.filterwarnings("ignore") + + +class Core(): + """ + All of the processing occurs here + """ + + def __init__(self): + # constructor + self.parameters = Parameters() # rsHRF parameters + self.dataStore = Store() # stores the information for the subjects + + # getters for data-store + def updateParameters(self, dic={}): + # updates the rsHRF-parameters + return self.parameters.set_parameters(dic) + + def get_parameters(self): + # returns the rsHRF-parameters + return self.parameters.get_parameters() + + def get_time_series(self, curr): + # gets the time-series + return self.dataStore.get_time_series(curr) + + def get_plotables(self, curr): + # gets all the time-series objects that can be plotted + return self.dataStore.get_plotables(curr) + + def get_store_info(self, curr): + # gets the information regarding the time-series objects + return self.dataStore.get_info(curr) + + def get_subjects(self): + # gets all the subjects present in the data-store + return self.dataStore.get_subjects() + + def get_data_labels(self, curr): + # gets all the labels present in the data-store + return self.dataStore.get_data_labels(curr) + + # saving data + def save_info(self, curr, out): + return self.dataStore.save_info(curr, out) + + def makeInput(self, inp): + """ + Obtains the input (file-paths and input-mode) from the main window + and updates the subjects/time-series objects appropriately + """ + # obtains the Raw BOLD data and preprocessed data, and updates the data-store + # each input comprises of an input_file, mask_file and a file_type, and the mode + input_file, mask_file, file_type, mode = inp + # if the mode is file (that is input is a stand-alone NIfTI or GIfTI file) + if mode == "file": + return self.makeInputFile(input_file, mask_file, file_type, 'file') + # if the input is present in a directory organized in BIDS format + elif "bids" in mode: + # obtaining all the input files + layout = BIDSLayout(input_file, validate=False, + config=['bids', 'derivatives']) + subjects_to_analyze = layout.get_subjects() + # if no subjects were found + if not subjects_to_analyze: + return Status(False, error='Could not find participants. Please make sure the BIDS data ' + 'structure is present and correct. Datasets can be validated online ' + 'using the BIDS Validator (http://incf.github.io/bids-validator/).') + all_inputs = layout.get(return_type='filename', datatype='func', subject=subjects_to_analyze, + task='rest', desc='preproc', suffix='bold', extension=['nii', 'nii.gz']) + # if no *bold.nii or *bold.nii.gz files were found + if not all_inputs != []: + return Status(False, error='There are no files of type *bold.nii / *bold.nii.gz ' + 'Please make sure to have at least one file of the above type ' + 'in the BIDS specification') + # if the atlas file is provided + if mode == "bids w/ atlas": + for file_count in range(len(all_inputs)): + # making input for every combination of input-file and atlas + temp = (self.makeInputFile( + all_inputs[file_count], mask_file, ".nii", 'bids w/ atlas')) + if not temp.get_state(): + return temp + # returns the status + return Status(True, info="Preprocessed all input", dic={"Number of Input Subjects: ": len(all_inputs)}) + # if the atlas file is not provided + elif mode == "bids": + # obtaining the masks + all_masks = layout.get(return_type='filename', datatype='func', subject=subjects_to_analyze, + task='rest', desc='brain', suffix='mask', extension=['nii', 'nii.gz']) + # if no mask files were found + if not all_masks != []: + return Status(False, error='There are no files of type *mask.nii / *mask.nii.gz ' + 'Please make sure to have at least one file of the above type ' + 'in the BIDS specification') + # if the number of mask files != number of input files + if len(all_inputs) != len(all_masks): + return Status(False, error='The number of *bold.nii / .nii.gz and the number of ' + '*mask.nii / .nii.gz are different. Please make sure that ' + 'there is one mask for each input_file present') + all_inputs.sort() + all_masks.sort() + # matching the prefix to align the inputs and the masks together + all_prefix_match = False + prefix_match_count = 0 + for i in range(len(all_inputs)): + input_prefix = all_inputs[i].split('/')[-1].split('_desc')[0] + mask_prefix = all_masks[i].split('/')[-1].split('_desc')[0] + + if input_prefix == mask_prefix: + prefix_match_count += 1 + else: + all_prefix_match = False + break + if prefix_match_count == len(all_inputs): + all_prefix_match = True + # if the mask and input files have different prefix + if not all_prefix_match: + return Status(False, error='The mask and input files should have the same prefix for correspondence. ' + 'Please consider renaming your files') + else: + # making input for every combination of input-file and mask-file + for file_count in range(len(all_inputs)): + temp = (self.makeInputFile( + all_inputs[file_count], all_masks[file_count], ".nii", 'bids')) + if not temp.get_state(): + return temp + # returns the status + return Status(True, info="Preprocessed all input", dic={"Number of Input Subjects: ": len(all_inputs)}) + + def makeInputFile(self, input_file, mask_file, file_type, mode): + """ + Obtains the Raw BOLD Time-series and Preprocessed BOLD Time-series + when the input is a stand-alone NIfTI or GIfTI file + """ + # if the input is in form of a File object + if 'bids' in mode: + input_file = input_file + # if the mask is in form of a File object + if mode == 'bids': + mask_file = mask_file + # getting the subject index + try: + subject_index = input_file.split("/")[-1].split("_")[0][4:] + except: + return Status(False, error="Input file should begin with 'sub-'") + # obtaining the header for the mask + try: + v = spm_dep.spm.spm_vol(mask_file) + except: + return Status(False, error="Invalid Mask File!") + # obtaining the header for the input + try: + v1 = spm_dep.spm.spm_vol(input_file) + except: + return Status(False, error="Invalid Input File!") + # obtaining the mask data + if file_type == ".nii" or file_type == ".nii.gz": + brain = spm_dep.spm.spm_read_vols(v) + elif file_type == ".gii" or file_type == ".gii.gz": + brain = v.agg_data().flatten(order='F') + else: + return Status(False, error="Invalid Input File Type!") + # brain voxel indices with as obtained by mask-data + voxel_ind = np.where(brain > 0)[0] + # checking the dimensions of the mask-file and input-file + if ((file_type == ".nii" or file_type == ".nii.gz") and v1.header.get_data_shape()[:-1] != v.header.get_data_shape()) or ((file_type == ".gii" or file_type == ".gii.gz") and v.agg_data().shape[0] != v.agg_data().shape[0]): + return Status(False, error='The dimension of your mask is different than the one of your fMRI data!') + # checking whether the subject is already present in the datastore + if not self.dataStore.is_present(subject_index): + # instantiating a new object if it is not present + subject = Subject(index=subject_index) + self.dataStore.add_subject(subject) + else: + # obtaining the object if it is present + subject = self.dataStore.get_subject_by_index(subject_index) + # checks if the BOLD time-series is already present + # if the time-series is not present + if not subject.is_present("BOLD", input_file): + if file_type == ".nii" or file_type == ".nii.gz": + data = v1.get_fdata() + nobs = data.shape[3] + data1 = np.reshape(data, (-1, nobs), order='F').T + elif file_type == ".gii" or file_type == ".gii.gz": + data = v1.agg_data() + _, nobs = data.shape + data1 = np.reshape(data, (-1, nobs), order='F').T + bold_ts = BOLD_Raw(label="BOLD", ts=np.array( + data1), para=self.parameters, subject_index=subject.get_subject_index()) + # the Raw BOLD Time-series is associated to a particular input-file + bold_ts.set_inputfile(input_file) + subject.add_BOLD_raw(bold_ts) # adding the time-series object + else: + bold_ts = subject.is_present("BOLD", input_file, getts=True) + data1 = bold_ts.get_ts() + # checks if the Preprocessed-BOLD time-series is already present + if not subject.is_present("Preprocessed-BOLD", (self.parameters, mask_file, bold_ts)): + # pre-processing the time-series + bold_sig = stats.zscore(data1[:, voxel_ind], ddof=1) + bold_sig = np.nan_to_num(bold_sig) + # instantiating the time-series objects + bold_pre_ts = BOLD_Preprocessed(label="Preprocessed-BOLD", ts=np.array( + bold_sig), para=self.parameters, subject_index=subject.get_subject_index()) + # the Preprocessed BOLD Time-series is associated to a particular Raw BOLD Time-Series + bold_pre_ts.set_BOLD_Raw(bold_ts) + bold_pre_ts.set_maskfile(mask_file) # and a particular mask_file + # adding the time-series objects + subject.add_BOLD_pre(bold_pre_ts) + # returning the status + return Status(True, info="Preprocessed Input!", dic={"Number of Time-Slices: ": bold_ts.get_shape()[0], "Number of Brain-Voxels: ": bold_ts.get_shape()[1]}) + return Status(False, error="Time series already exists. No new additions were made.") + + def retrieveHRF(self, bold_pre_ts, get_pos=False): + """ + Retrieves the resting-state Hemodynamic Response Function + with Preprocessed BOLD Time-series as input, and self.parameters + as the parameters + """ + subject_index = bold_pre_ts.get_subject_index( + ) # gets the subject-index associated to the preprocessed BOLD Time-series + subject = self.dataStore.get_subject_by_index( + subject_index) # gets the subject from the index + # if the HRF has already been retrieved for this particular set of inputs + if not subject.is_present("HRF", (self.parameters, bold_pre_ts)): + # inputs for retrieving the HRF + bold_sig = bold_pre_ts.get_ts() + bold_sig = processing. \ + rest_filter. \ + rest_IdealFilter(bold_sig, self.parameters.get_TR( + ), self.parameters.get_passband_deconvolve()) + para = self.parameters.get_parameters() + if not (para['estimation'] == 'sFIR' or para['estimation'] == 'FIR'): + # estimate HRF for the fourier / hanning / gamma / cannon basis functions + bf = basis_functions.basis_functions.get_basis_function( + bold_sig.shape, para) # obtaining the basis set + beta_hrf, event_bold = utils.hrf_estimation.compute_hrf( + bold_sig, para, [], -1, bf) + hrfa = np.dot(bf, beta_hrf[np.arange(0, bf.shape[1]), :]) + else: + # estimate HRF for FIR and sFIR + beta_hrf, event_bold = utils.hrf_estimation.compute_hrf( + bold_sig, para, [], -1) + hrfa = beta_hrf + # instantiating the time-series objects + hrf = HRF(label="HRF", ts=hrfa, + subject_index=subject_index, para=self.parameters) + # the HRF is associated to a particular Preprocessed BOLD Time-series + hrf.set_BOLD(bold_pre_ts) + hrf.set_event_bold(event_bold) # setting the bold-events + pos = subject.add_HRF(hrf) # adding the HRF to the subject + # returning the status + if get_pos: + return Status(True, info="Retrieved HRF!", dic={"Shape of HRF: ": hrf.get_shape()}), pos + else: + return Status(True, info="Retrieved HRF!", dic={"Shape of HRF: ": hrf.get_shape()}) + if get_pos: + return Status(False, error="Time series already exists. No new additions were made."), None + else: + return Status(False, error="Time series already exists. No new additions were made.") + + def getHRFParameters(self, hrf): + """ + Retrieves the resting-state Hemodynamic Response Function Parameters + with HRF as input, and self.parameters as the parameters + """ + # checking if the parameters have already been computed for HRF + if hrf.get_HRF_para().size != 0: + return Status(False, error="Parameters have already been computed for this HRF") + # inputs for obtaining the parameters + hrfa = hrf.get_ts() + para = hrf.get_parameters().get_parameters() + nvar = hrfa.shape[1] + PARA = np.zeros((3, nvar)) + average = [0 for i in range(3)] + # obtaining the parameters for each voxel + for voxel_id in range(nvar): + hrf1 = hrfa[:, voxel_id] + PARA[:, voxel_id] = parameters.wgr_get_parameters( + hrf1, para['TR'] / para['T']) + for i in range(3): + average[i] += PARA[i, voxel_id] + for i in range(3): + average[i] /= nvar + # setting the parameters + hrf.set_para(PARA) + # returning the status + return Status(True, info="Retrieved HRF Parameters", dic={"Average Height: ": str(average[0])[:3], + "Avearage Time-To-Peak: ": str(average[1])[:3] + " seconds", "Avearge Full-Width at Half-Max: ": str(average[2])[:3]}) + + def deconvolveHRF(self, hrf): + """ + Retrieves the Deconvolved BOLD Time-series + with HRF as input, and self.parameters as the parameters + """ + subject_index = hrf.get_subject_index( + ) # gets the subject-index associated to the preprocessed BOLD Time-series + subject = self.dataStore.get_subject_by_index( + subject_index) # gets the subject from the index + para = self.parameters.get_parameters() + # if the HRF has already been retrieved for this particular set of inputs + if not subject.is_present("Deconvolved-BOLD", (self.parameters, hrf)): + # inputs for retrieving the deconvolved BOLD + hrfa = hrf.get_ts() + bold_sig = hrf.get_associated_BOLD().get_ts() + bold_sig = processing. \ + rest_filter. \ + rest_IdealFilter(bold_sig, self.parameters.get_TR( + ), self.parameters.get_passband_deconvolve()) + event_bold = hrf.get_event_bold() + nvar = hrfa.shape[1] + nobs = bold_sig.shape[0] + data_deconv = np.zeros(bold_sig.shape) + event_number = np.zeros((1, bold_sig.shape[1])) + if para['T'] > 1: + hrfa_TR = signal.resample_poly(hrfa, 1, para['T']) + else: + hrfa_TR = hrfa + # obtaining the deconvolved BOLD for each voxel + for voxel_id in range(nvar): + hrf_ = hrfa_TR[:, voxel_id] + H = np.fft.fft(np.append(hrf_, np.zeros( + (nobs - max(hrf_.shape), 1))), axis=0) + M = np.fft.fft(bold_sig[:, voxel_id]) + data_deconv[:, voxel_id] = np.fft.ifft( + H.conj() * M / (H * H.conj() + .1*np.mean((H * H.conj())))) + event_number[:, voxel_id] = np.amax(event_bold[voxel_id].shape) + # instantiating the time-series object + dd = Bold_Deconv(label="Deconvolved-BOLD", ts=data_deconv, + subject_index=subject_index, para=hrf.get_parameters()) + # the deconvoled BOLD Time-series is associated to a particular HRF + dd.set_HRF(hrf) + dd.set_event_num(event_number) # setting the BOLD events + # adding the time-series to the subject + pos = subject.add_BOLD_deconv(dd) + # returning the status + return Status(True, info="Deconvolved BOLD Signal") + return Status(False, error="Time series already exists. No new additions were made.") diff --git a/build/lib/rsHRF/rsHRF_GUI/datatypes/__init__.py b/build/lib/rsHRF/rsHRF_GUI/datatypes/__init__.py new file mode 100644 index 0000000..a811498 --- /dev/null +++ b/build/lib/rsHRF/rsHRF_GUI/datatypes/__init__.py @@ -0,0 +1,2 @@ +from . import misc +from . import timeseries \ No newline at end of file diff --git a/build/lib/rsHRF/rsHRF_GUI/datatypes/misc/__init__.py b/build/lib/rsHRF/rsHRF_GUI/datatypes/misc/__init__.py new file mode 100644 index 0000000..92df576 --- /dev/null +++ b/build/lib/rsHRF/rsHRF_GUI/datatypes/misc/__init__.py @@ -0,0 +1,3 @@ +from . import parameters +from . import subject +from . import store diff --git a/build/lib/rsHRF/rsHRF_GUI/datatypes/misc/parameters.py b/build/lib/rsHRF/rsHRF_GUI/datatypes/misc/parameters.py new file mode 100644 index 0000000..a6a5de5 --- /dev/null +++ b/build/lib/rsHRF/rsHRF_GUI/datatypes/misc/parameters.py @@ -0,0 +1,442 @@ +import sys +import numpy as np +from copy import deepcopy +from ...misc.status import Status + +class Parameters(): + """ + rsHRF-Toolbox Parameters + For more information, visit https://github.com/BIDS-Apps/rsHRF/tree/master/rsHRF + """ + def __init__(self): + # initialize default parameters + self.estimation = 'canon2dd' + self.passband = [0.01, 0.08] + self.passband_deconvolve = [0.0, sys.float_info.max] + self.TR = 2.0 + self.localK = 1 + self.T = 3 + self.T0 = 1 + self.TD_DD = 2 + self.AR_lag = 1 + self.thr = 1 + self.order = 3 + self.volterra = 0 + self.len = 24 + self.temporal_mask = [] + self.min_onset_search = 4 + self.max_onset_search = 8 + self.dt = self.TR/self.T + self.lag = np.arange(np.fix(self.min_onset_search / self.dt), + np.fix(self.max_onset_search / self.dt) + 1, + dtype='int') + + # getters + def get_estimation(self): + return self.estimation + + def get_passband(self): + return deepcopy(self.passband) + + def get_passband_deconvolve(self): + return deepcopy(self.passband_deconvolve) + + def get_TR(self): + return self.TR + + def get_localK(self): + return self.localK + + def get_T(self): + return self.T + + def get_T0(self): + return self.get_T0 + + def get_TD_DD(self): + # TD_DD is only relevant for canon2dd estimation + if self.estimation == 'canon2dd': + return self.TD_DD + else: + return None + + def get_AR_lag(self): + return self.get_AR_lag + + def get_thr(self): + return self.thr + + def get_order(self): + # order is only relevant for fourier and gamma estimation + if 'gamma' in self.get_estimation or 'fourier' in self.get_estimation: + return self.order + else: + return None + + def get_Volterra(self): + # volterra is only relevant for canon2dd estimation + if self.get_estimation == 'canon2dd': + return self.volterra + else: + return None + + def get_len(self): + return self.len + + def get_temporal_mask(self): + return tuple(self.temporal_mask) + + def get_min_onset_search(self): + return self.min_onset_search + + def get_max_onset_search(self): + return self.max_onset_search + + # setters + def set_estimation(self,s): + self.estimation = s + return Status(True) + + def set_passband(self, l): + """ + Takes a list (with two entries) as input and sets the + passband-range + Checks whether both the values are non-negative + """ + try: + l = [float(i) for i in l.split(",")] + if l[0] < 0 or l[1] < 0: + return Status(False, error="Passband frequency values cannot be negative") + except: + return Status(False, error="Bad datatype for passband range") + else: + self.passband = deepcopy(l) + return Status(True) + + def set_passband_deconvolve(self, l): + """ + Takes a list (with two entries) as input and sets the + passband-range + Checks whether both the values are non-negative + """ + try: + l = [float(i) for i in l.split(",")] + if l[0] < 0 or l[1] < 0: + return Status(False, error="Passband frequency values cannot be negative") + except: + return Status(False, error="Bad datatype for passband range") + else: + self.passband_deconvolve = deepcopy(l) + return Status(True) + + def set_TR(self, TR): + """ + Checks whether the value is postiive + and updates 'dt' as it is dependent on TR + """ + try: + TR = float(TR) + except: + return Status(False,error="Bad Datatype For TR") + else: + if TR <= 0: + return Status(False, error="BOLD Repetition Time must be greater than 0") + else: + self.TR = TR + self.update_dt() + return Status(True) + + def set_localK(self, localK): + """ + Checks whether the value of localK is positive + """ + try: + localK = int(localK) + except: + return Status(False,error="Bad Datatype For localK") + else: + if localK <= 0: + return Status(False, error="localK must be greater than 0") + else: + self.localK = localK + return Status(True) + + def set_T(self, T): + """ + Checks whether the value of T >= 1 + and updates 'dt' as it is dependent on T + """ + try: + T = int(T) + except: + return Status(False, error="Bad Datatype For T") + else: + if T < 1: + return Status(False, error="Magnification factor must not be less than 1") + else: + self.T = T + self.update_dt() + return Status(True) + + def set_T0(self, T0): + try: + T0 = int(T0) + except: + return Status(False, error="Bad Datatype For T0") + else: + self.T0 = T0 + return Status(True) + + def set_TD_DD(self, TD_DD): + """ + Sets the value of TD_DD if the estimation is canon2dd + """ + if self.estimation == 'canon': + try: + TD_DD = int(TD_DD) + except: + return Status(False, error="Magnification factor must not be less than 1") + if TD_DD not in [0,1,2]: + return Status(False,error="TD_DD can only take one of these values: 0, 1, 2") + else: + self.TD_DD = TD_DD + return Status(True) + else: + return Status(True) + + def set_AR_lag(self, AR_lag): + """ + Checks whether the AR_lag is non-negative + """ + try: + AR_lag = int(AR_lag) + except: + return Status(False, error="Bad datatype for AR_lag") + else: + if AR_lag < 0: + return Status(False, error="AR_lag must not be negative") + else: + self.AR_lag = AR_lag + return Status(True) + + def set_thr(self, thr): + """ + If estimation is FIR or sFIR, thr is a list + Otherwise, it is an int + """ + if 'FIR' not in self.estimation : + try: + thr = int(thr) + except: + return Status(False, error="Bad datatype for thr") + else: + self.thr = thr + return Status(True) + else: + try: + thr = [int(i) for i in thr.split(",")] + except: + return Status(False, error="Bad datatype for thr") + else: + self.thr = thr + return Status(True) + + def set_order(self, order): + """ + Checks whether the order lies between 1 to 60 + """ + try: + order = int(order) + except: + return Status(False, error="Bad datatype for order") + else: + if order < 0 or order > 60: + status = Status(False) + if order < 0: + status.set_error("Order must not be negative") + elif order > 60: + status.set_error("Order is too high") + return status + else: + self.order = order + return Status(True) + + def set_Volterra(self, volterra): + """ + Sets Volterra if the estimation rule is canon2dd + """ + if self.estimation == 'canon': + try: + volterra = int(volterra) + except: + return Status(False, error="Bad datatype for Volterra") + else: + self.volterra = volterra + return Status(True) + else: + return Status(True) + + def set_len(self, length): + """ + Checks whether the length of the hemodynamic response function is non-negative + """ + try: + length = int(length) + except: + return Status(False, error="Bad datatype for len") + else: + if length <= 0: + status = Status(False) + status.set_error("HRF length must not be postitive") + return status + else: + self.len = length + return Status(True) + + def set_temporal_mask(self, temporal_mask): + self.temporal_mask = deepcopy(temporal_mask) + return Status(True) + + def set_min_onset_search(self, mos): + """ + Checks whether the min-onset-search is smaller than max-onset-search + and non-negative + """ + try: + mos = int(mos) + except: + return Status(False, error="Bad datatype for min_onset_search") + else: + if mos > self.max_onset_search or mos < 0: + status = Status(False) + if mos > self.max_onset_search: + status.set_error("Min Onset Search must not be greater than Max Onset Search") + elif mos < 0: + status.set_error("Onset Search values must not be negative") + return status + else: + self.min_onset_search = mos + self.update_lag() + return Status(True) + + def set_max_onset_search(self, mos): + """ + Checks whether max-onset-search is greater than min-onset-search + and non-negative + """ + try: + mos = int(mos) + except: + return Status(False, error="Bad datatype for max_onset_search") + else: + if mos < self.min_onset_search or mos < 0: + status = Status(False) + if mos < self.min_onset_search: + status.set_error("Max Onset Search must not be lesser than Min Onset Search") + elif mos < 0: + status.set_error("Onset Search values must not be negative") + return status + else: + self.max_onset_search = mos + self.update_lag() + return Status(True) + + def update_dt(self): + """ + Re-calculating dt + """ + self.dt = self.TR / self.T + self.update_lag() + + def update_lag(self): + """ + Re-calculating lag + """ + self.lag = np.arange(np.fix(self.min_onset_search / self.dt), + np.fix(self.max_onset_search / self.dt) + 1, + dtype='int') + + def get_parameters(self): + """ + Gets all the parameters in the form of a dictionary for rsHRF computation + """ + para = {} + para['estimation'] = self.estimation + para['passband'] = deepcopy(self.passband) + para['passband_deconvolve'] = deepcopy(self.passband_deconvolve) + para['TR'] = self.TR + para['T'] = self.T + para['localK'] = self.localK + para['T0'] = self.T0 + para['AR_lag'] = self.AR_lag + para['thr'] = self.thr + para['len'] = self.len + para['temporal_mask'] = deepcopy(self.temporal_mask) + para['min_onset_search'] = self.min_onset_search + para['max_onset_search'] = self.max_onset_search + if self.estimation == 'canon2dd': + para['TD_DD'] = self.TD_DD + para['Volterra'] = self.volterra + elif 'gamma' in self.estimation or 'fourier' in self.estimation: + para['order'] = self.order + para['dt'] = self.dt + para['lag'] = self.lag + return para + + def set_parameters(self, dic): + """ + Takes a dictionary as input and sets all the rsHRF parameters accordingly + """ + for key in dic.keys(): + if key == "estimation": + out = self.set_estimation(dic[key]) + elif key == "passband": + out = self.set_passband(dic[key]) + elif key == "passband_deconvolve": + out = self.set_passband_deconvolve(dic[key]) + elif key == "T": + out = self.set_T(dic[key]) + elif key == "TR": + out = self.set_TR(dic[key]) + elif key == "localK": + out = self.set_localK(dic[key]) + elif key == "T0": + out = self.set_T0(dic[key]) + elif key == "TD_DD": + out = self.set_TD_DD(dic[key]) + elif key == "AR_lag": + out = self.set_AR_lag(dic[key]) + elif key == "thr": + out = self.set_thr(dic[key]) + elif key == "order": + out = self.set_order(dic[key]) + elif key == "Volterra": + out = self.set_thr(dic[key]) + elif key == "len": + out = self.set_len(dic[key]) + elif key == "temporal_mask": + out = self.set_temporal_mask(dic[key]) + elif key == "min_onset_search": + out = self.set_min_onset_search(dic[key]) + elif key == "max_onset_search": + out = self.set_max_onset_search(dic[key]) + if not out.get_state(): + return out + return Status(True, info="Parameters Updated Succefully") + + def compareParameters(self, p): + """ + Takes another parameter object and determines if it is equal + to the this object + """ + x = self.get_parameters() + y = p.get_parameters() + for key in x.keys(): + if key not in y.keys(): + return False + elif key == "dt" or key == "lag": + continue + else: + if x[key] != y[key]: + return False + return True diff --git a/build/lib/rsHRF/rsHRF_GUI/datatypes/misc/store.py b/build/lib/rsHRF/rsHRF_GUI/datatypes/misc/store.py new file mode 100644 index 0000000..3c8ad0c --- /dev/null +++ b/build/lib/rsHRF/rsHRF_GUI/datatypes/misc/store.py @@ -0,0 +1,80 @@ +from . import subject + +class Store(): + """ + Stores all the subjects as a dictionary, with their keys as the index + """ + def __init__(self): + self.subjects = {} + + # getters + def get_subjects(self): + """ Returns a tuple with the indices of all subjects""" + return tuple(self.subjects.keys()) + + def get_subject_by_index(self, index): + """ Takes subject-index as the input and returns the corresponding subject """ + try: + return self.subjects[index] + except: + return None + + def get_plotables(self, index): + """ Gets all the time-series that can be plotted, for a particular subject""" + return self.subjects[index].get_plotables() + + def get_data_labels(self, index): + """ Gets the labels for all the time-series data for that subject""" + return self.subjects[index].get_data_labels() + + # misc. + def get_time_series(self, s): + """ Gets the time-series object corresponding to s""" + s = s.split("_") + index = s[0] + ts_type = s[1] + ts_num = int(s[2]) + return self.subjects[index].get_time_series_by_index(ts_type, ts_num) + + def get_info(self, s): + """ Gets all the information about a time-series to be displayed to the logging-window """ + ts = self.get_time_series(s) + subject_index = ts.get_subject_index() + subject = self.subjects[subject_index] + dic = ts.get_info() + if ts.get_label() == "BOLD": + dic["Input File"] = ts.get_inputfile() + elif ts.get_label() == "Preprocessed-BOLD": + dic["Mask File"] = ts.get_maskfile() + dic["Associated Raw BOLD"] = subject_index + "_BOLD_" + subject.get_time_series_pos(ts.get_BOLD_Raw()) + elif ts.get_label() == "HRF": + dic["Associated BOLD"] = subject_index + "_Preprocessed-BOLD_" + subject.get_time_series_pos(ts.get_associated_BOLD()) + elif ts.get_label() == "Deconvolved-BOLD": + dic["Associated HRF"] = subject_index + "_HRF_" + subject.get_time_series_pos(ts.get_associated_HRF()) + return dic + + def is_present(self, subject_index): + """ Checks whether a subject-index is present is already present""" + if subject_index in self.subjects.keys(): + return True + return False + + def add_subject(self, sub): + """ Adds a new subject """ + self.subjects[sub.get_subject_index()] = sub + + def remove_subject(self, sub): + """ Removing a subject """ + try: + del self.subjects[sub.get_subject_index()] + except: + return + + def number_of_subjects(self): + """ Gets the number of subjects """ + return len(self.subjects.keys()) + + def save_info(self, s, out): + """ Saves all the time-series objects for a particular subject in the form of .mat files """ + ts = self.get_time_series(s) + return ts.save_info(out+"/sub-"+s+".mat") \ No newline at end of file diff --git a/build/lib/rsHRF/rsHRF_GUI/datatypes/misc/subject.py b/build/lib/rsHRF/rsHRF_GUI/datatypes/misc/subject.py new file mode 100644 index 0000000..956a829 --- /dev/null +++ b/build/lib/rsHRF/rsHRF_GUI/datatypes/misc/subject.py @@ -0,0 +1,180 @@ +import numpy as np + +class Subject(): + """ + Stores the information corresponding to a particular subject. + + Attrbutes: + 1. index : This is the index of the subject (as determined by BIDS convention) + 2. BOLD_raw : An array which stores the corresponding Raw BOLD time-series for the subject + 3. BOLD_pre : An array which stores the corresponding Preprocessed-BOLD time-series for the subject + 4. BOLD_deconv : An array which stores the corresponding Deconvolved-BOLD time-series for the subject + 5. HRF : An array which stores the corresponding Hemodynamic Response Function time-series for the subject + + -> All the attributes from 2-5, are arrays of TimeSeries objects + """ + def __init__(self, index): + self.index = index + self.BOLD_raw = [] + self.BOLD_pre = [] + self.BOLD_deconv = [] + self.HRF = [] + + # getters + def get_input_filename(self): + return self.input_filename + + def get_subject_index(self): + return self.index + + def get_BOLD_raw(self): + return tuple(self.BOLD_raw) + + def get_BOLD_pre(self): + return tuple(self.BOLD_pre) + + def get_BOLD_deconv(self): + return tuple(self.BOLD_deconv) + + def get_HRF(self): + return tuple(self.HRF) + + # adding to time-series objects of the existing subject + def add_BOLD_raw(self, ts): + self.BOLD_raw.append(ts) + return len(self.BOLD_raw) - 1 + + def add_BOLD_deconv(self, ts): + self.BOLD_deconv.append(ts) + return len(self.BOLD_raw) - 1 + + def add_BOLD_pre(self, ts): + self.BOLD_pre.append(ts) + return len(self.BOLD_raw) - 1 + + def add_HRF(self, ts): + self.HRF.append(ts) + return len(self.BOLD_raw) - 1 + + # misc. + def is_present(self, label, misc, getts=False): + """ Checks whether a time-series is already present + Misc takes in all the relevant information which determines the uniqueness + of a time-series + If getts = True, the function returns the time-series object if it is present """ + if label == "BOLD": + # Looking for Raw BOLD Data + for each in self.BOLD_raw: + # Determines whether the Raw BOLD data is already present + # Checks for the input-file + if misc == each.get_inputfile(): + if getts : + return each + return True + elif label == "Preprocessed-BOLD": + # Looking for Preprocessed BOLD Data + para = misc[0] + mask = misc[1] + bold = misc[2] + for each in self.BOLD_pre: + # Determines whether Preprocessed BOLD data is already present + # Checks the parameters, mask-file and RAW Bold + if para.compareParameters(each.get_parameters()) \ + and each.get_maskfile() == misc[1] \ + and bold.compareTimeSeries(each.get_BOLD_Raw()): + if getts: + return each + return True + elif label == "HRF": + # Looking for HRF Data + para = misc[0] + BOLD_pre = misc[1] + for each in self.HRF: + # Determines whether the HRF is already present + # Checks the parameters and Preprocessed BOLD + if para.compareParameters(each.get_parameters()) \ + and BOLD_pre.compareTimeSeries(each.get_associated_BOLD()): + if getts: + return each + return True + elif label == "Deconvolved-BOLD": + # Looking for Deconvolved BOLD Data + para = misc[0] + HRF = misc[1] + for each in self.BOLD_deconv: + # Determines whether the Deconvolved BOLD is already present + # Checks the associated HRF + if para.compareParameters(each.get_parameters()) \ + and HRF.compareTimeSeries(each.get_associated_HRF()): + if getts : + return each + return True + return False + + def get_time_series_pos(self, ts): + """ + Takes the time-series as input and returns its position in the array + """ + label = ts.get_label() + if label == "BOLD": + arr = self.BOLD_raw + elif label == "Preprocessed-BOLD": + arr = self.BOLD_pre + elif label == "Deconvolved-BOLD": + arr = self.BOLD_deconv + elif label == "HRF": + arr = self.HRF + else : + arr = [] + for i in range(len(arr)): + if ts.compareTimeSeries(arr[i]): + return str(i) + return None + + def get_time_series_by_index(self, ts_type, index): + """ Takes the index of a time-series and returns the time-series """ + if ts_type == "BOLD": + arr = self.BOLD_raw + elif ts_type == "Preprocessed-BOLD": + arr = self.BOLD_pre + elif ts_type == "Deconvolved-BOLD": + arr = self.BOLD_deconv + elif ts_type == "HRF": + arr = self.HRF + else: + return + return arr[index] + + def get_plotables(self): + """ + Returns an array of all the time-series objects that can be plotted for the subject + The array contains of tuples of the format : (time-series labels, time-series numpy arrays) + """ + out = [] + for i in range(len(self.BOLD_raw)): + out.append((self.index+"_BOLD_"+str(i),self.BOLD_raw[i].get_ts())) + for i in range(len(self.BOLD_pre)): + out.append((self.index+"_Preprocessed-BOLD_"+str(i),self.BOLD_pre[i].get_ts())) + for i in range(len(self.BOLD_deconv)): + out.append((self.index+"_Deconvolved-BOLD_"+str(i),self.BOLD_deconv[i].get_ts())) + for i in range(len(self.HRF)): + out.append((self.index+"_HRF_"+str(i),self.HRF[i].get_ts())) + return out + + def get_data_labels(self): + """ + Returns an array with labels for all the time-series objects for the subject + """ + out = [] + for i in range(len(self.BOLD_raw)): + out.append(self.index+"_BOLD_"+str(i)) + for i in range(len(self.BOLD_pre)): + out.append(self.index+"_Preprocessed-BOLD_"+str(i)) + for i in range(len(self.BOLD_deconv)): + out.append(self.index+"_Deconvolved-BOLD_"+str(i)) + for i in range(len(self.HRF)): + out.append(self.index+"_HRF_"+str(i)) + return out + + + diff --git a/build/lib/rsHRF/rsHRF_GUI/datatypes/timeseries/__init__.py b/build/lib/rsHRF/rsHRF_GUI/datatypes/timeseries/__init__.py new file mode 100644 index 0000000..0815fa4 --- /dev/null +++ b/build/lib/rsHRF/rsHRF_GUI/datatypes/timeseries/__init__.py @@ -0,0 +1,3 @@ +from . import bold_deconv +from . import timeseries +from . import hrf diff --git a/build/lib/rsHRF/rsHRF_GUI/datatypes/timeseries/bold_deconv.py b/build/lib/rsHRF/rsHRF_GUI/datatypes/timeseries/bold_deconv.py new file mode 100644 index 0000000..1937fc1 --- /dev/null +++ b/build/lib/rsHRF/rsHRF_GUI/datatypes/timeseries/bold_deconv.py @@ -0,0 +1,62 @@ +import numpy as np +from scipy.io import savemat +from copy import deepcopy + +from .hrf import HRF +from .timeseries import TimeSeries +from ...datatypes.misc.parameters import Parameters + +class Bold_Deconv(TimeSeries): + """ + This stores the Deconvolved BOLD Time-series + + Attributes: + 1. HRF : The HRF used to obtain the Deconvolved BOLD + 2. event_num : The event-numbers + """ + def __init__(self, label="",ts=np.array([]),subject_index="", para=Parameters()): + TimeSeries.__init__(self, label="",ts=np.array([]),subject_index="", para=Parameters()) + self.label = label + self.subject_index = subject_index + self.timeseries = ts + self.shape = ts.shape + self.parameters = deepcopy(para) + self.HRF = HRF() + self.event_num = np.array([]) + + #setters + def set_HRF(self, HRF): + self.HRF = HRF + def set_event_num(self, ev): + self.event_num = ev + + #getters + def get_event_num(self): + return self.event_num + def get_associated_HRF(self): + return self.HRF + + #misc + def compareTimeSeries(self, ts): + """ Compares another time-series with itself to determine if both are identical + Two checks are performed: + 1. Label + 2. HRF + If all the three comparisions return true, then both the HRF + time-series objects are identical + """ + if self.label == ts.get_label() \ + and self.HRF.compareTimeSeries(ts.get_associated_HRF()): + return True + else: + return False + def save_info(self, name): + """ Saves the information about the time-series in a .mat file """ + dic = {} + dic["timeseries"] = self.timeseries + dic["eventNum"] = self.event_num + para = self.parameters.get_parameters() + for each in para.keys(): + dic[each] = para[each] + savemat(name, dic) + return True \ No newline at end of file diff --git a/build/lib/rsHRF/rsHRF_GUI/datatypes/timeseries/bold_preprocessed.py b/build/lib/rsHRF/rsHRF_GUI/datatypes/timeseries/bold_preprocessed.py new file mode 100644 index 0000000..20e3bb7 --- /dev/null +++ b/build/lib/rsHRF/rsHRF_GUI/datatypes/timeseries/bold_preprocessed.py @@ -0,0 +1,69 @@ +import numpy as np +from scipy.io import savemat +from copy import deepcopy + +from ...datatypes.misc.parameters import Parameters +from .timeseries import TimeSeries +from .bold_raw import BOLD_Raw + +class BOLD_Preprocessed(TimeSeries): + """ + This stores the Preprocessed BOLD Time-series + + Attributes: + 1. BOLD_Raw : The Raw BOLD time-series object through which it was derived + 2. mask_file : The mask-file path + """ + def __init__(self, label="",ts=np.array([]),subject_index="", para=Parameters()): + TimeSeries.__init__(self, label="",ts=np.array([]),subject_index="", para=Parameters()) + self.label = label + self.subject_index = subject_index + self.timeseries = ts + self.shape = ts.shape + self.parameters = deepcopy(para) + self.BOLD_Raw = BOLD_Raw() + self.mask_file = "" + + # setters + def set_maskfile(self, mask_file): + self.mask_file = mask_file + + def set_BOLD_Raw(self, BOLD_Raw): + self.BOLD_Raw = BOLD_Raw + + # getters + def get_maskfile(self): + return self.mask_file + + def get_BOLD_Raw(self): + return self.BOLD_Raw + + # misc. + def compareTimeSeries(self, ts): + """ Compares another time-series with itself to determine if both are identical + Four checks are performed: + 1. Label + 2. Parameters + 3. Raw BOLD associated with it + 4. Mask-file + If all the three comparisions return true, then both the HRF + time-series objects are identical + """ + if self.label == ts.get_label() \ + and self.parameters == ts.get_parameters() \ + and self.BOLD_Raw.compareTimeSeries(ts.get_BOLD_Raw()) \ + and ts.get_maskfile() == self.mask_file: + return True + else: + return False + + def save_info(self, name): + """ Saves the information about the time-series in a .mat file """ + try: + dic = {} + dic["timeseries"] = self.timeseries + dic["mask_file"] = self.mask_file + savemat(name, dic) + return True + except: + return False \ No newline at end of file diff --git a/build/lib/rsHRF/rsHRF_GUI/datatypes/timeseries/bold_raw.py b/build/lib/rsHRF/rsHRF_GUI/datatypes/timeseries/bold_raw.py new file mode 100644 index 0000000..5f81a53 --- /dev/null +++ b/build/lib/rsHRF/rsHRF_GUI/datatypes/timeseries/bold_raw.py @@ -0,0 +1,55 @@ +import numpy as np +from scipy.io import savemat +from copy import deepcopy + +from ...datatypes.misc.parameters import Parameters +from .timeseries import TimeSeries + +class BOLD_Raw(TimeSeries): + """ + This stores the Raw BOLD Time-series + + Attributes: + 1. input_file : the input-file path of the Raw BOLD Time-series + """ + def __init__(self, label="",ts=np.array([]),subject_index="", para=Parameters()): + TimeSeries.__init__(self, label="",ts=np.array([]),subject_index="", para=Parameters()) + self.label = label + self.subject_index = subject_index + self.timeseries = ts + self.shape = ts.shape + self.parameters = deepcopy(para) + self.input_file = "" + + # setters + def set_inputfile(self, input_file): + self.input_file = input_file + + # getters + def get_inputfile(self): + return self.input_file + + # misc. + def compareTimeSeries(self, ts): + """ Compares another time-series with itself to determine if both are identical + Two checks are performed: + 1. Label + 2. Input-file name + If all the three comparisions return true, then both the HRF + time-series objects are identical + """ + if self.label == ts.get_label() and self.input_file == ts.get_inputfile(): + return True + else: + return False + + def save_info(self, name): + """ Saves the information about the time-series in a .mat file """ + try: + dic = {} + dic["timeseries"] = self.timeseries + dic["input_file"] = self.input_file + savemat(name, dic) + return True + except: + return False diff --git a/build/lib/rsHRF/rsHRF_GUI/datatypes/timeseries/hrf.py b/build/lib/rsHRF/rsHRF_GUI/datatypes/timeseries/hrf.py new file mode 100644 index 0000000..2d14b65 --- /dev/null +++ b/build/lib/rsHRF/rsHRF_GUI/datatypes/timeseries/hrf.py @@ -0,0 +1,76 @@ +import numpy as np +from scipy.io import savemat +from copy import deepcopy + +from ...datatypes.misc.parameters import Parameters +from .timeseries import TimeSeries + +class HRF(TimeSeries): + """ + This stores the Hemodynamic Response Function Time-series + + Attributes: + 1. BOLD = stores the associated Preprocessed-BOLD time-series object through which it was retrieved + 2. PARA = stores the HRF parameters (Full-width at half-max, Time-to-peak and height) + 3. event_bold = stores the bold-events + """ + def __init__(self, label="",ts=np.array([]),subject_index="", para=Parameters()): + TimeSeries.__init__(self, label="",ts=np.array([]),subject_index="", para=Parameters()) + self.label = label + self.subject_index = subject_index + self.timeseries = ts + self.shape = ts.shape + self.parameters = deepcopy(para) + self.BOLD = TimeSeries() + self.PARA = np.array([]) + self.event_bold = np.array([]) + + # setters + def set_para(self, PARA): + self.PARA = PARA + + def set_BOLD(self, BOLD): + self.BOLD = BOLD + + def set_event_bold(self, event_bold): + self.event_bold = event_bold + + # getters + def get_event_bold(self): + return self.event_bold + + def get_associated_BOLD(self): + return self.BOLD + + def get_HRF_para(self): + return self.PARA + + # misc. + def compareTimeSeries(self, ts): + """ Compares whether another HRF time-series is similar to it. + Three checks are performed for this: + 1. Label + 2. rsHRF Parameters + 3. Preprocessed-BOLD which was used as input + If all the three comparisions return true, then both the HRF + time-series objects are identical + """ + if self.label == ts.get_label() and self.parameters.compareParameters(ts.get_parameters()) and self.BOLD.compareTimeSeries(ts.get_associated_BOLD()): + return True + else: + return False + + def save_info(self, name): + """ Saves the information about the time-series in a .mat file """ + try: + dic = {} + dic["timeseries"] = self.timeseries + dic["PARA"] = self.PARA + dic["eventBold"] = self.event_bold + para = self.parameters.get_parameters() + for each in para.keys(): + dic[each] = para[each] + savemat(name, dic) + return True + except: + return False \ No newline at end of file diff --git a/build/lib/rsHRF/rsHRF_GUI/datatypes/timeseries/timeseries.py b/build/lib/rsHRF/rsHRF_GUI/datatypes/timeseries/timeseries.py new file mode 100644 index 0000000..5bde3f1 --- /dev/null +++ b/build/lib/rsHRF/rsHRF_GUI/datatypes/timeseries/timeseries.py @@ -0,0 +1,77 @@ +import numpy as np +from copy import deepcopy + +from ...datatypes.misc.parameters import Parameters + +class TimeSeries(): + """ + Deals with all the time-series artifacts that appear during the processing. + + These are the various time-series that are dealt with: + 1. BOLD : (Raw BOLD Data) + 2. HRF : (Hemodynamic Response Function) + 3. Preprocessed-BOLD : (Z-score and Passband filtered) + 4. Deconvolved-BOLD : (Deconvolved using the HRF) + + Attributes: + 1. label : Takes on one of the values as described above + 2. subject_index : index of the subject to which the current time-series belongs + 3. timeseries : stores the time-series as a 2-dimensional numpy array + 4. shape : shape of the time-series (voxels x time-slices) + 5. parameters : rsHRF-parameters associated in retrieving this time-series + """ + def __init__(self, label="",subject_index="",ts=np.array([]), para=Parameters()): + self.label = label + self.subject_index = subject_index + self.timeseries = deepcopy(ts) + self.shape = ts.shape + self.parameters = deepcopy(para) + + # setters + def set_ts(self, ts): + self.timeseries = deepcopy(ts) + self.shape = self.timeseries.shape + + def set_parameters(self,para): + self.parameters = deepcopy(para) + + def set_label(self,label): + self.label = label + + def set_subject_index(self,subject_index): + self.subject_index = subject_index + + # getters + def get_ts(self): + return self.timeseries + + def get_subject_index(self): + return self.subject_index + + def get_label(self): + return self.label + + def get_parameters(self): + return self.parameters + + def get_shape(self): + return self.shape + + # misc. + def get_info(self): + """ Returns the information about the time-series in the form of a dictionary """ + dic = {} + dic["Type"] = self.label + dic["Subject"] = self.subject_index + dic["Time Series Shape"] = self.shape + dic["Parameters"] = self.parameters.get_parameters() + return dic + + def compareTimeSeries(self, ts): + """ Compares another time-series with itself to determine if both are identical """ + raise NotImplementedError("This needs to be overridden in the child-classes") + + def save_info(self, name): + """ Saves the information about the time-series in a .mat file """ + raise NotImplementedError("This needs to be overridden in the child-classes") + diff --git a/build/lib/rsHRF/rsHRF_GUI/documentation/developer/class_diagram.png b/build/lib/rsHRF/rsHRF_GUI/documentation/developer/class_diagram.png new file mode 100755 index 0000000..158392c Binary files /dev/null and b/build/lib/rsHRF/rsHRF_GUI/documentation/developer/class_diagram.png differ diff --git a/build/lib/rsHRF/rsHRF_GUI/documentation/developer/class_info.md b/build/lib/rsHRF/rsHRF_GUI/documentation/developer/class_info.md new file mode 100755 index 0000000..21f29a1 --- /dev/null +++ b/build/lib/rsHRF/rsHRF_GUI/documentation/developer/class_info.md @@ -0,0 +1,350 @@ +# Class-Structure Information + +Here we have tried to define the purpose and usage of each class relavant to the rsHRF-GUI. + +--- + +## rsHRF_GUI.core.Core + +This class in concerned with all the technical data-processing w.r.t the rsHRF toolbox. + +### Class Variables +- `self.parameters` - an instance of the Parameters() object +- `self.dataStre` - an instance of the Store() object + +--- + +### Class Methods +- `updateParameters(self, dic={})` - takes a dictionary and updates the *self.parameters* variables.
+- `get_parameters(self)` - returns the rsHRF parameters as a dictionary +- `get_time_series(self, curr)` - returns the time-series object as defined by the string variable *curr*. +- `get_plotables(self, curr)` - returns the plotable artefact corresponding to the string variable *curr*. +- `get_store_info(self, curr)` - returns the information about the time-series object corresponding to the string variable *curr*. +- `get_subjects(self)` - returns a list of all the subjects that are present in the *self.dataStore*. +- `get_data_labels(self, curr)` - returns a label string corresponding to the string variable *curr*. +- `save_info(self, curr, out)` - save sthe info regarding the time-series corresponding to the variable *curr*, in the directory as specified using *out*. +- `makeInput(self, inp)` - takes different forms of input files and forms the *RAW-BOLD* and *Preprocessed-BOLD* time-series objects, and saves it with respect to their corresponding subjects in the *self.dataStore*. +- `makeInputFile(self, input_file, mask_file, file_type, mode)` - deals with one-input-instance at a time. +- `retrieveHRF(self, bold_pre_ts, get_pos=False)` - retrieves the rsHRF corresponding to a preprocessed-BOLD time-series. +- `getHRFParameters(self, hrf)` - computes the corresponding rsHRF parameters with respct to an HRF time-series. +- `deconvolveHRF(self, hrf)` - uses the retrieved rsHRF to deconvolve the preprocessed BOLD time-series. + +--- + +## rsHRF_GUI.datatypes.misc.parameters.Parameters +This class corresponds to all the rsHRF-toolbox parameters. + +### Class Variables +All the information regarding the class-variables can be found along with the main documentation of the rsHRF-toolbox. +- `self.estimation` +- `self.passband` +- `self.passband_deconvolve` +- `self.TR` +- `self.localK` +- `self.T` +- `self.T0` +- `self.TD_DD` +- `self.AR_lag` +- `self.thr` +- `self.order` +- `self.volterra` +- `self.len` +- `self.temporal_mask` +- `self.min_onset_search` +- `self.max_onset_search` +- `self.dt` +- `self.lag` + +--- + +### Class Methods +**Getters** +- `get_estimation(self)` +- `get_passband(self)` +- `get_passband_deconvolve(self)` +- `get_TR(self)` +- `get_localK(self)` +- `get_T(self)` +- `get_T0(self)` +- `get_TD_DD(self)` +- `get_AR_lag(self)` +- `get_thr(self)` +- `get_order(self)` +- `get_volterra(self)` +- `get_len(self)` +- `get_temporal_mask(self)` +- `get_min_onset_search(self)` +- `get_max_onset_search(self)` +- `get_dt(self)` +- `get_lag(self)` + +--- + +**Setters** + +All the setters take the corresponding error-handling into account. If an illegal value is provided, a suitable error gets returned.
+Also, if setting there are variables that get affected by the one being set, the change gets propogated automatically. + +- `set_estimation(self, var)` +- `set_passband(self, var)` +- `set_passband_deconvolve(self, var)` +- `set_TR(self, var)` +- `set_localK(self, var)` +- `set_T(self, var)` +- `set_T0(self, var)` +- `set_TD_DD(self, var)` +- `set_AR_lag(self, var)` +- `set_thr(self, var)` +- `set_order(self, var)` +- `set_volterra(self, var)` +- `set_len(self, var)` +- `set_temporal_mask(self, var)` +- `set_min_onset_search(self, var)` +- `set_max_onset_search(self, var)` + +--- + +These parameters are dependent parameters, and hence, these are automatically updates as and when required. + +- `update_dt(self)` +- `update_lag(self)` + +--- + +- `get_parameters(self)` - returns all the parameters as a dictionary. +- `set_parameters(self, dic)` - takes a dictionary as input and sets the parameters. +- `compareParameters(self, p)` - takes another *parameters* object and compares whether it is similar as this parameters object or not. + +--- + +## rsHRF_GUI.datatypes.misc.store.Store +This corresponds to where all the *subject* objects are dealt with. + +### Class Variables + +- `self.subjects` - a dictionary which stores the subjects against their labels (the labels are the string value as in BIDS format). + +--- + +### Class Methods + +- `get_subjects(self)` - returns the labels for all the subjects. +- `get_subject_by_index(self, index)` - takes the subject index and returns the subject object. +- `get_plotables(self, index)` - takes the subject index and returns all the plotables artefats associated to that subject. +- `get_data_labels(self, index)` - returns labels for all the relavant artefacts corresponding to the subject given by the index. +- `get_time_series(self, s)` - returns the time-series object corresponding to the string *s*. +- `get_info(self, s)` - return the information regarding the time-series object corresponding to the string *s*. +- `is_present(self, subject_index)` - checks whether the subject corresponding to an index value is already present. +- `add_subject(self, sub)` - takes a *subject* object and adds it into the store. +- `remove_subject(self, sub)` - takes a *subject* object and removes it from the store. +- `number_of_subjects(self)` - returns the number of subjects corrently in the store. +- `save_info(self, s, out)` - saves the information regarding the subject as denoted by *s*, into the directory corresponding to *out*. + +--- + +## rsHRF_GUI.datatypes.misc.subject.Subject +This corresponds to every individual subject. + +### Class Variables +- `self.index` - index of the subject (as in BIDS format) +- `self.BOLD_raw` - list of all the RAW BOLD time-series data. +- `self.BOLD_pre` - list of the preprocessed BOLD time-series data. +- `self.BOLD_deconv` - list of all the deconvolved BOLD time-series data. +- `self.HRF` - list of all the HRF time-series data. + +--- + +### Class Methods +**Getters** +- `get_input_filename(self)` +- `get_subject_index(self)` +- `get_BOLD_raw(self)` +- `get_BOLD_pre(self)` +- `get_BOLD_deconv(self)` +- `get_HRF(self)` + +--- + +**Adding Data** + +Adds time-series data to the existing collection. + +- `add_BOLD_raw(self, ts)` +- `add_BOLD_deconv(self, ts)` +- `add_BOLD_pre(self, ts)` +- `add_HRF(self, ts)` + +--- + +- `is_present(self, label, misc, getts=False)` - checks whether a time-series object is already present or not. +- `get_time_series_pas(self, ts)` - gets the position of the time-series object in the list that it is stored. +- `get_time_series_by_index(self, ts_type, index)` - gets the time-series object as specified by the *index*. +- `get_plotables(self)` - returns all the plotable time-series artefacts. +- `get_data_labels(self)` - returns the lables for all the time-series artefacts. + +--- + +## rsHRF_GUI.datatypes.timeseries.TimeSeries +This deals with the time-series objects. +### Class Variables +- `self.label` - label of the time-series +- `self.subject_index` - subject-index of the subject that the time-series is associated to. +- `self.timeseries` - the time-series stored as a numpy array. +- `self.shape` - shape of the time-series +- `self.parameters` - rsHRF parameters required while obtaining the time-series. + +--- + +### Class Methods +** Setters ** +- `set_ts(self, ts)` +- `set_parameters(self, para)` +- `set_label(self, label)` +- `set_subject_index(self, subject_index)` + +--- + +** Getters *8 +- `get_ts(self)` +- `get_subject_index(self)` +- `get_label(self)` +- `get_parameters(self)` +- `get_shape(self)` +- `get_info(self)` + +--- + +- `compareTimeSeries(self, ts)` - compares the two-timeseries to figure out whether they are identicle. +- `save_info(self, name)` - saves the information of the time-series in a .mat file. + +--- + +**Note:** The *rsHRF_GUI.datatypes.timeseries* package contains other classes: +
+1. Bold_Deconv +2. Bold_Preprocessed +3. Bold_Raw +4. HRF + +All these classes inherit from the time-series classes, and suitable class variables, getters, setters, are added, and appropriate functions are overridden. However, the basic structure for all these classes remains the same. + +--- + +## rsHRF_GUI.misc.status.Status + +This class allows for making the GUI more interactive. It provides an inherent way of communication across the different modules. + +### Class Variables +- `self.state` - a boolean value +- `self.info` - information string regarding the message +- `self.error` - error string corresponding to the message +- `self.time` - stores the time of the message +- `self.dic` - a dictionary for storing misc. information about the message + +--- + +### Class Methods +** Setters ** +- `set_state(self, s)` +- `set_info(self, s)` +- `set_error(self, e)` +- `set_time(self, t)` +- `set_dic(self, d)` + +--- + +** Getters ** +- `get_state(self)` +- `get_info(self)` +- `get_error(self)` +- `get_time(self)` +- `get_dic(self)` + +--- + +## rsHRF_GUI.misc.gui_windows.inputWindow.InputWindow + +This class is concerned with accepting the input, determining whether its in a valid format, and subsequently passing a valid input forward. + +### Class Variables +- `self.input_file` - stores the path to input file +- `self.mask_file` - stores the path to mask file +- `self.file_type` - stores the type of the file +- `self.output_dir` - stores the path to the output directory + +--- + +### Class Methods + +- `getInput(self)` - retrieves the input from the GUI and passes it to `main.py` + +--- + +## rsHRF_GUI.misc.gui_windows.loggingWindow.LoggingWindow + +### Class Variables +- `self.lineNum` - the position of the current line in the logging window +- `self.text` - ScrolledText window + +--- + +### Class Methods +- `putLog(self, plotables=None, data_info=None, status=None)` - depending upon the input, appropriately formats it and displays it on the logging-screen. + +--- + +## rsHRF_GUI.misc.gui_windows.parameterWindow.ParameterWindow + +### Class Variables +- `self.window` - Toplevel instance +- `self.parameters` - dictionary containing all the rsHRF parameters +- `self.labels` - labels corresponding to each parameter +- `self.entries` - current value of each parameter + +--- + +### Class Methods +- `updateParameters(self)` - updates the parameter dictionary +- `getParameters(self)` - gets the parameter dictionary +- `setParameters(self)` - sets the parameter dictionary +- `display(self)` - puts the elements of the parameter dictionary into *self.labels* and *self.entries* in a way that they can be displayed by the logging-window. + +--- + +## rsHRF_GUI.misc.gui_windows.plotterOptionsWindow.PlotterOptionsWindow + +### Class Variables +- `self.plotterScreen` - PlotterWindow instance +- `self.numberOfPlots` - the number of plots on the plotter-screen +- `self.plot` - string corresponding to the time-series which is to be plotted +- `self.plotVal` - voxel value of the time-series that is to be plotted +- `self.plotValStore` - decides whether the current plot is to be displayed or not +- `self.plotables` - list of strings corresponding to all the time-series that can be plotted +- `self.options` - to display the various plots that can be plotted in the drop-down menu + +--- + +### Class Methods +- `updatePlots(self, plotables)` - takes an iterable (plotables) and reconfigures the various time-series that can be plotted +- `updateWidgets(self)` - updates the widgets on the *plotter-window* +- `get_plotables(self)` - returns the plotables + +--- + +## rsHRF_GUI.misc.gui_windows.plotterWindow.PlotterWindow + +### Class Variables +- `self.numberOfPlots` - the number of plots that can be displayed at a time (right now = 3) +- `self.ts` - list of time-series that are to be plotted +- `self.plot` - displaying the plots +- `self.canvas` - FigureCanvasTkAgg instance + +--- + +### Class Methods +- `getNumberOfPlots(self)` - gives the number of plots +- `makePlots(self, ts, val, num)` - displays the appropriate plots + +--- + +## Fin. \ No newline at end of file diff --git a/build/lib/rsHRF/rsHRF_GUI/documentation/rsHRF-GUI-User-Manual.pdf b/build/lib/rsHRF/rsHRF_GUI/documentation/rsHRF-GUI-User-Manual.pdf new file mode 100755 index 0000000..a7cc1f8 Binary files /dev/null and b/build/lib/rsHRF/rsHRF_GUI/documentation/rsHRF-GUI-User-Manual.pdf differ diff --git a/build/lib/rsHRF/rsHRF_GUI/gui_windows/__init__.py b/build/lib/rsHRF/rsHRF_GUI/gui_windows/__init__.py new file mode 100644 index 0000000..c137f03 --- /dev/null +++ b/build/lib/rsHRF/rsHRF_GUI/gui_windows/__init__.py @@ -0,0 +1,6 @@ +from . import inputWindow +from . import loggingWindow +from . import parameterWindow +from . import plotterOptionsWindow +from . import plotterWindow +from . import main \ No newline at end of file diff --git a/build/lib/rsHRF/rsHRF_GUI/gui_windows/inputWindow.py b/build/lib/rsHRF/rsHRF_GUI/gui_windows/inputWindow.py new file mode 100644 index 0000000..2268c33 --- /dev/null +++ b/build/lib/rsHRF/rsHRF_GUI/gui_windows/inputWindow.py @@ -0,0 +1,105 @@ +import os +from tkinter import Toplevel, Checkbutton, IntVar, Button, filedialog, NORMAL, DISABLED, OptionMenu, StringVar, Label + + +class InputWindow(): + def __init__(self): + # input window + window = Toplevel() + window.title("Input Window") + # get screen width and height + screen_width = window.winfo_screenwidth() + screen_height = window.winfo_screenheight() + # placing the toplevel + window.geometry("350x220+%d+%d" % ((280/1900)*screen_width, ((((1040.0-220)/1000)*screen_height)-390)-280)) + # variables which shall get sent to the front end + self.input_file = () + self.mask_file = () + self.file_type = () + self.output_dir = () + # other class vairables + # 1 corresponds to BIDS input + self.inputFormatVar = IntVar() + self.inputFormatVar.set(0) + # 1 corresponds to mask file being present in the BIDS directory + self.maskFormatVar = IntVar() + self.maskFormatVar.set(0) + # selecting the estimation rule + self.estimationOption = StringVar() + self.estimationOption.set('canon2dd') + + def getInputFile(): + if self.inputFormatVar.get(): # input takes a directory + self.input_file = filedialog.askdirectory(initialdir=os.getcwd()) + maskFormat.configure(state=NORMAL) + else: + self.input_file = filedialog.askopenfilename(initialdir=os.getcwd(), title="Input File Path", filetypes=(("nifti files", "*.nii"), ("nifti files", "*.nii.gz"), ("gifti files", "*.gii"), ("gifti files", "*.gii.gz"))) + maskFormat.configure(state=DISABLED) + try: + self.file_type = self.input_file.split(os.extsep,1)[-1] + self.file_type = '.' + self.file_type + except: + self.file_type = () + try: + inputFileLabel.configure(text=self.input_file.split('/')[-1]) + except: + inputFileLabel.configure(text="") + + def maskFormatButtonState(): + if self.maskFormatVar.get(): + maskFormatButton.configure(state=DISABLED) + else: + maskFormatButton.configure(state=NORMAL) + + def inputFormatButtonState(): + if self.inputFormatVar.get(): + maskFormat.configure(state=NORMAL) + else: + maskFormat.configure(state=DISABLED) + maskFormatButtonState() + + def getMaskFile(): + self.mask_file = filedialog.askopenfilename(initialdir=os.getcwd(), title="Input File Path", filetypes=(("nifti files", "*.nii"), ("nifti files", "*.nii.gz"), ("gifti files", "*.gii"), ("gifti files", "*.gii.gz"))) + try: + maskFileLabel.configure(text=self.mask_file.split("/")[-1]) + except: + maskFileLabel.configure(text="") + + def getOutputDir(): + self.output_dir = filedialog.askdirectory(initialdir=os.getcwd()) + try: + outputPathLabel.configure(text="Output path: " + self.output_dir.split("/")[-1]) + except: + outputPathLabel.configure(text="") + + # defining widgets + inputFormat = Checkbutton(window, text="BIDS Format", variable=self.inputFormatVar, command=inputFormatButtonState) + maskFormat = Checkbutton(window, text="Mask File in BIDS", variable=self.maskFormatVar, state=DISABLED, command=maskFormatButtonState) + inputFormatButton = Button (window, text="Select Input", command=getInputFile, height=1, width=20) + maskFormatButton = Button (window, text="Select Mask File", state=NORMAL, command=getMaskFile, height=1, width=20) + outputPathButton = Button (window, text="Select Output Directory", command=getOutputDir, height=1, width=20) + estimationLabel = Label (window, text="Estimation Rule: ") + inputFileLabel = Label (window, text="") + maskFileLabel = Label (window, text="") + outputPathLabel = Label (window, text="") + estimationDropDown = OptionMenu (window, self.estimationOption, "canon2dd", "sFIR", "FIR", "gamma", "fourier", "fourier w/ hanning") + # placing widgets + inputFormat.grid (row=0,column=0,padx=(5,5),pady=(5,5)) + inputFormatButton.grid (row=0,column=1,padx=(5,5),pady=(5,5)) + inputFileLabel.grid (row=1,column=0,padx=(5,5),pady=(5,5),columnspan=2) + maskFormat.grid (row=2,column=0,padx=(5,5),pady=(5,5)) + maskFormatButton.grid (row=2,column=1,padx=(5,5),pady=(5,5)) + maskFileLabel.grid (row=3,column=0,padx=(5,5),pady=(5,5),columnspan=2) + outputPathButton.grid (row=4,column=1,padx=(5,5),pady=(5,5)) + outputPathLabel.grid (row=4,column=0,padx=(5,5),pady=(5,5)) + estimationLabel.grid (row=5,column=0,padx=(5,5),pady=(5,5)) + estimationDropDown.grid(row=5,column=1,padx=(5,5),pady=(5,5)) + + def getInput(self): + if self.inputFormatVar.get() * self.maskFormatVar.get(): + mode = "bids" + elif self.inputFormatVar.get(): + mode = "bids w/ atlas" + else: + mode = "file" + return (self.input_file, self.mask_file, self.file_type, mode, self.estimationOption.get(), self.output_dir) \ No newline at end of file diff --git a/build/lib/rsHRF/rsHRF_GUI/gui_windows/loggingWindow.py b/build/lib/rsHRF/rsHRF_GUI/gui_windows/loggingWindow.py new file mode 100644 index 0000000..7c3b2c5 --- /dev/null +++ b/build/lib/rsHRF/rsHRF_GUI/gui_windows/loggingWindow.py @@ -0,0 +1,129 @@ +from ..misc.status import Status +from datetime import datetime +from tkinter import Toplevel, BOTH, INSERT +from tkinter.scrolledtext import ScrolledText + +class LoggingWindow(): + def __init__(self): + window = Toplevel(bg='white') + window.title("Logger") + # get screen width and height + screen_width = window.winfo_screenwidth() + screen_height = window.winfo_screenheight() + window.geometry("600x400+%d+%d" % (screen_width-600, 0)) + self.lineNum = 1 + self.text = ScrolledText(window, background='black') + self.text.pack(fill=BOTH, expand = 1) + + def putLog(self, plotables=None, data_info=None, status=None): + # writing current time + current_time = datetime.now().strftime("%H:%M:%S") + + if status is not None: + if status.get_state(): + message = status.get_info() + self.text.insert(INSERT, \ + ">> " + current_time + "\t SUCCESS " + message + "\n" + ) + self.text.tag_add("success", str(self.lineNum) + ".12", str(self.lineNum) + ".20") + self.text.tag_config("success", foreground="green") + self.text.tag_add("success message", str(self.lineNum) + ".20", str(self.lineNum) + "." + str(21 + len(message))) + self.text.tag_config("success message", foreground="LightBlue1") + else: + message = status.get_error() + self.text.insert(INSERT, \ + ">> " + current_time + "\t ERROR " + message + "\n" + ) + self.text.tag_add("error", str(self.lineNum) + ".12", str(self.lineNum) + ".18") + self.text.tag_config("error", foreground="red") + self.text.tag_add("error message", str(self.lineNum) + ".18", str(self.lineNum) + "." + str(19 + len(message))) + self.text.tag_config("error message", foreground="yellow") + self.text.tag_add("time", str(self.lineNum) + ".03", str(self.lineNum) + ".11") + self.text.tag_config("time", foreground="PeachPuff2") + self.lineNum += 1 + + dic = status.get_dic() + for each in dic.keys(): + self.text.insert(INSERT, \ + ">> " + current_time + "\t" + each + str(dic[each]) + "\n" + ) + self.text.tag_add("time", str(self.lineNum) + ".03", str(self.lineNum) + ".11") + self.text.tag_config("time", foreground="PeachPuff2") + self.text.tag_add("dic_label", str(self.lineNum) + ".12", str(self.lineNum) + "." + str(12+len(each))) + self.text.tag_config("dic_label", foreground="azure") + self.lineNum += 1 + self.text.insert(INSERT, "\n") + self.lineNum += 1 + return + + elif plotables is not None: + self.text.insert(INSERT, \ + ">> " + current_time + "\t Time Series' Information \n\n" + \ + "TYPE" + "\t\t\t\t\t\t LENGTH" + "\t\t NUMBER OF VOXELS \n\n" + ) + self.text.tag_add("labels", str(self.lineNum) + ".12", str(self.lineNum) + "." + str(37)) + self.text.tag_config("labels", foreground="PaleTurquoise1") + self.text.tag_add("titles", str(self.lineNum+2) + ".00", str(self.lineNum+2) + "." + str(37)) + self.text.tag_config("titles", foreground="azure") + self.text.tag_add("time", str(self.lineNum) + ".03", str(self.lineNum) + ".11") + self.text.tag_config("time", foreground="PeachPuff2") + for each in plotables: + self.text.insert(INSERT, \ + each[0] + "\t\t\t\t\t\t " + str(each[1].shape[0]) + "\t\t " + str(each[1].shape[1]) + "\n" + ) + self.lineNum += 1 + self.text.insert(INSERT, "\n") + self.lineNum += 5 + return + + elif data_info is not None: + self.text.insert(INSERT, \ + ">> " + current_time + "\t Information \n\n") + self.text.tag_add("labels", str(self.lineNum) + ".12", str(self.lineNum) + "." + str(37)) + self.text.tag_config("labels", foreground="PaleTurquoise1") + self.text.tag_add("time", str(self.lineNum) + ".03", str(self.lineNum) + ".11") + self.text.tag_config("time", foreground="PeachPuff2") + self.lineNum += 2 + if data_info["Type"] == "BOLD": + self.text.insert(INSERT, \ + " :- " + "Type of Data: BOLD Time Series\n" + + " :- " + "Subject: " + data_info["Subject"] + "\n" + " :- " + "Input File: " + data_info["Input File"] + "\n" + ) + self.lineNum += 1 + elif data_info["Type"] == "Preprocessed-BOLD": + self.text.insert(INSERT, \ + " :- " + "Type of Data: Preprocessed-BOLD Time Series\n" + + " :- " + "Subject: " + data_info["Subject"] + "\n" + + " :- " + "Mask File: " + data_info["Mask File"] + "\n" + + " :- " + "Associated Raw BOLD: " + data_info["Associated Raw BOLD"] + "\n" + ) + self.lineNum += 2 + elif data_info["Type"] == "HRF": + self.text.insert(INSERT, \ + " :- " + "Type of Data: Hemodynamic Response Function Time Series\n" + + " :- " + "Associated BOLD: " + data_info["Associated BOLD"] + "\n" + ) + elif data_info["Type"] == "Deconvolved-BOLD": + self.text.insert(INSERT, \ + " :- " + "Type of Data: Deconvolved BOLD Time Series\n" + + " :- " + "Associated HRF: " + data_info["Associated HRF"] + "\n" + ) + self.text.tag_config("type", foreground="azure") + self.text.tag_config("line2", foreground="azure") + self.lineNum += 2 + self.text.insert(INSERT, \ + " :- " + "\t\t\tAssociated Parameters\n" + ) + self.text.tag_add("ap", str(self.lineNum) + ".03", str(self.lineNum) + ".30") + self.text.tag_config("ap", foreground="Khaki1") + self.lineNum += 1 + for each in data_info["Parameters"].keys(): + self.text.insert(INSERT, \ + "\t\t" + each + "\t\t\t=\t\t" + str(data_info["Parameters"][each]) + "\n") + self.text.tag_add("pv", str(self.lineNum) + ".00", str(self.lineNum) + "." + str(6+len(each))) + self.text.tag_config("pv", foreground="LemonChiffon2") + self.lineNum +=1 + self.text.insert(INSERT, "\n") + self.lineNum += 1 + return \ No newline at end of file diff --git a/build/lib/rsHRF/rsHRF_GUI/gui_windows/main.py b/build/lib/rsHRF/rsHRF_GUI/gui_windows/main.py new file mode 100644 index 0000000..3cbe9bc --- /dev/null +++ b/build/lib/rsHRF/rsHRF_GUI/gui_windows/main.py @@ -0,0 +1,286 @@ +import os +import time +from scipy.io import savemat +from tkinter import Tk, Button, Label, OptionMenu, StringVar + +from .inputWindow import InputWindow +from .loggingWindow import LoggingWindow +from .plotterOptionsWindow import PlotterOptionsWindow +from .parameterWindow import ParameterWindow + +from ..misc.status import Status +from ..core.core import Core + +class Main(): + def __init__(self): + # main window + root = Tk() + root.title("rsHRF Toolbox") + # get screen width and height + screen_width = root.winfo_screenwidth() + screen_height = root.winfo_screenheight() + # placing the toplevel + root.geometry("500x300+%d+%d" % (((280/1900)*screen_width + 425), (((1040.0-220)/1000)*screen_height)-670)) + # defining the toplevels + input_window = InputWindow() + parameter_window = ParameterWindow() + core = Core() + logger = LoggingWindow() + plotter = PlotterOptionsWindow() + # defining variables used by tkinter widgets + current = StringVar() + current_subject = StringVar() + output_path = os.getcwd() + self.options = ["None"] + self.subjects = ["None"] + current.set ("None") + current_subject.set("None") + # other variables + input = () # receives the input from the input window + output = {} # receives the output from the core + # initializing parameter window + parameter_window.setParameters(core.get_parameters()) + parameter_window.display() + + ''' Gets the input from the input toplevel. + The input comprises of three elements: + 1. Input file (resting-state fMRI) + 2. Brain Mask + 3. Estimation Rule + ''' + def makeInput(): + start = time.process_time() + # interacting with the input window + input = input_window.getInput() + # interacting with the back-end + parameters = parameter_window.getParameters() + parameters['estimation'] = input[-2] + core.updateParameters(parameters) # only passing the estimation-rule + output = core.makeInput(input[:-2]) # receiving input from core + output_path = input[-1] # the path to output directory + # interacting with the parameter-window + parameter_window.setParameters(core.get_parameters()) + parameter_window.display() + # updating data + updateSubjects() + # logging + time_taken = time.process_time() - start + output.set_time(str(time_taken)[:5]) + logger.putLog(status=output) + + def changeEstimationRule(): + start = time.process_time() + # interacting with the input window + input = input_window.getInput() + # interacting with the core + output = core.updateParameters({"estimation":input[-2]}) + # interacting with the parameter-window + parameter_window.setParameters(core.get_parameters()) + parameter_window.display() + # logging + time_taken = time.process_time() - start + output.set_time(str(time_taken)[:5]) + logger.putLog(status=output) + + def updateParameters(): + start = time.process_time() + # interacting with the core + output = core.updateParameters(dic=parameter_window.updateParameters()) + # interacting with the parameter-window + parameter_window.setParameters(core.get_parameters()) + parameter_window.display() + time_taken = time.process_time() - start + # logging + output.set_time(str(time_taken)[:5]) + logger.putLog(status=output) + + def retrieveHRF(get_pos=False): + if "Preprocessed-BOLD" not in current.get(): + logger.putLog(status=Status(False, error="Select a Preprocessed-BOLD Timeseries to retrieve HRF")) + else: + start = time.process_time() + core.updateParameters(dic =parameter_window.updateParameters()) + bold_ts = core.get_time_series(current.get()) + output = core.retrieveHRF(bold_ts, get_pos) + # updating data + updateValues(log=False) + # logging + time_taken = time.process_time() - start + if get_pos: + output, pos = output + else: + pos = None + output.set_time(str(time_taken)[:5]) + logger.putLog(status=output) + return pos + + def retrieveHRFParameters(pos=None): + if pos != None: + current.set(current.get().split("_")[0] + "_HRF_" + str(pos)) + if "HRF" not in current.get(): + logger.putLog(status=Status(False, error="Select an HRF Timeseries to retrieve parameters")) + else: + start = time.process_time() + # interacting with the core + output = core.getHRFParameters(core.get_time_series(current.get())) + # updating data + updateValues(log=False) + # logging + time_taken = time.process_time() - start + output.set_time(str(time_taken)[:5]) + logger.putLog(status=output) + + def deconvolveHRF(pos=None): + if pos != None: + current.set(current.get().split("_")[0] + "_HRF_" + str(pos)) + if "HRF" not in current.get(): + logger.putLog(status=Status(False, error="Select an HRF Timeseries to deconvolve BOLD")) + else: + start = time.process_time() + # interacting with the core + output = core.deconvolveHRF(core.get_time_series(current.get())) + # updating data + updateValues(log=False) + # logging + time_taken = time.process_time() - start + output.set_time(str(time_taken)[:5]) + logger.putLog(status=output) + + def runCompletePipeline(): + if "Preprocessed-BOLD" in current.get(): + pos = retrieveHRF(get_pos=True) + retrieveHRFParameters(pos) + deconvolveHRF(pos) + elif "HRF" in current.get(): + retrieveHRFParameters() + deconvolveHRF() + else: + logger.putLog(status=Status(False, error="Select an HRF or Preprocessed-BOLD Timeseries to run pipeline")) + + def updatePlots(): + # interacting with the core + try: + plotables = core.get_plotables(current_subject.get()) + except: + logger.putLog(status=Status(False, error="Please select a subject")) + else: + # logging + if plotables == None: + logger.putLog(status=Status(False, error="Please select a subject")) + else: + logger.putLog(plotables=plotables) + # interacting with the plotter + plotter.updatePlots(plotables) + + def getInfo(): + try: + if len(current.get().split("_")) == 3: + logger.putLog(data_info=core.get_store_info(current.get())) + else: + logger.putLog(status = Status(False, error="Select a valid input to get information")) + except: + logger.putLog(status = Status(False, error="Select a valid input to get information")) + + def saveValue(): + try: + if len(current.get().split("_")) == 3: + if os.path.isdir(output_path): + out = output_path + else: + out = os.cwd() + if core.save_info(current.get(), out): + logger.putLog(status = Status(True, info="File saved successfully")) + else: + logger.putLog(status = Status(False, error="Unable to save file")) + else: + logger.putLog(status = Status(False, error="Select a valid input to save")) + except: + logger.putLog(status = Status(False, error="Select a valid input to save")) + + def saveForSubject(): + try: + temp = core.get_data_labels(current_subject.get()) + except: + logger.putLog(status=Status(False, error="Please select a subject")) + else: + for each in temp: + current.set(each) + saveValue() + + def add_new(l1, l2): + if l1 == None: + return l2 + if l2 == None: + return l1 + out = [] + for each in l1: + if each not in l2: + out.append(each) + for each in l2: + out.append(each) + return out + + def updateSubjects(): + temp = list(core.get_subjects()) + temp.append ("None") + self.subjects = add_new(self.subjects, temp) + self.subjects = sorted(self.subjects) + subjectOptions = OptionMenu(root, current_subject, *self.subjects) + subjectOptions.grid (row=0,column=1,padx=(30,30),pady=(5,5)) + + def updateValues(log=True): + try: + temp = core.get_data_labels(current_subject.get()) + except: + logger.putLog(status=Status(False, error="Please select a subject")) + else: + if log: + logger.putLog(status=Status(True, info="Updated values for subject: " + current_subject.get())) + temp.append ("None") + self.options = add_new(self.options, temp) + self.options = sorted(self.options) + valueOptions = OptionMenu(root, current, *self.options) + valueOptions.grid (row=1,column=1,padx=(30,30),pady=(5,5)) + + + + # defining the widgets + getInput = Button (root, text="Get Input", command=makeInput, height = 1, width = 20) + changeParameterButton = Button (root, text="Update Parameters", command=updateParameters, height = 1, width = 20) + retrieveHRFButton = Button (root, text="Retrieve HRF", command=retrieveHRF, height = 1, width = 20) + retrieveHRFParametersButton = Button (root, text="Retrieve HRF Parameters", command=retrieveHRFParameters, height = 1, width = 20) + deconvolveHRFButton = Button (root, text="Deconvolve BOLD", command=deconvolveHRF, height = 1, width = 20) + updatePlotsButton = Button (root, text="Update Plots", command=updatePlots, height = 1, width = 20) + changeEstimationRule = Button (root, text="Set Estimation Rule", command=changeEstimationRule, height = 1, width = 20) + getValueInfo = Button (root, text="Get Info", command=getInfo, height = 1, width = 20) + storeValue = Button (root, text="Save Value", command=saveValue, height = 1, width = 20) + updateStore = Button (root, text="Update Values", command=updateValues, height = 1, width = 20) + runCompletePipeline = Button (root, text="Run Pipeline", command=runCompletePipeline, height = 1, width = 20) + saveAllValues = Button (root, text="Save All", command=saveForSubject, height = 1, width = 20) + dataLabel = Label (root, text="Stored Values: ") + subjectLabel = Label (root, text="Subjects: ") + valueOptions = OptionMenu(root, current, *self.options) + subjectOptions = OptionMenu(root, current_subject, *self.subjects) + # placing the widgets + subjectLabel.grid (row=0,column=0,padx=(30,30),pady=(5,5)) + subjectOptions.grid (row=0,column=1,padx=(30,30),pady=(5,5)) + dataLabel.grid (row=1,column=0,padx=(30,30),pady=(5,5)) + valueOptions.grid (row=1,column=1,padx=(30,30),pady=(5,5)) + getInput.grid (row=2,column=0,padx=(30,30),pady=(5,5)) + updateStore.grid (row=2,column=1,padx=(30,30),pady=(5,5)) + changeEstimationRule.grid (row=3,column=0,padx=(30,30),pady=(5,5)) + changeParameterButton.grid (row=3,column=1,padx=(30,30),pady=(5,5)) + retrieveHRFButton.grid (row=4,column=0,padx=(30,30),pady=(5,5)) + retrieveHRFParametersButton.grid(row=4,column=1,padx=(30,30),pady=(5,5)) + deconvolveHRFButton.grid (row=5,column=0,padx=(30,30),pady=(5,5)) + updatePlotsButton.grid (row=5,column=1,padx=(30,30),pady=(5,5)) + getValueInfo.grid (row=6,column=0,padx=(30,30),pady=(5,5)) + storeValue.grid (row=6,column=1,padx=(30,30),pady=(5,5)) + runCompletePipeline.grid (row=7,column=0,padx=(30,30),pady=(5,5)) + saveAllValues.grid (row=7,column=1,padx=(30,30),pady=(5,5)) + + root.mainloop() + + +if __name__ == "__main__": + Main() \ No newline at end of file diff --git a/build/lib/rsHRF/rsHRF_GUI/gui_windows/parameterWindow.py b/build/lib/rsHRF/rsHRF_GUI/gui_windows/parameterWindow.py new file mode 100644 index 0000000..f1bd56b --- /dev/null +++ b/build/lib/rsHRF/rsHRF_GUI/gui_windows/parameterWindow.py @@ -0,0 +1,54 @@ +from copy import deepcopy +from tkinter import Toplevel, Button, Entry, Label, DISABLED + +class ParameterWindow(): + def __init__(self): + self.window = Toplevel() + self.window.title("Parameters") + self.parameters = {} + self.labels = [] + self.entries = [] + # get screen width and height + screen_width = self.window.winfo_screenwidth() + screen_height = self.window.winfo_screenheight() + self.window.geometry("350x420+%d+%d" % (screen_width*(280.0/1900), (((1040.0-220)/1000)*screen_height)-390)) + + def updateParameters(self): + for i in range(len(self.labels)): + self.parameters[self.labels[i].cget('text')] = self.entries[i].get() + parameters = dict(self.parameters) + return deepcopy(parameters) + + def getParameters(self): + return deepcopy(self.parameters) + + def setParameters(self, dic): + self.parameters = deepcopy(dic) + + def display(self): + # removing existing widgets + for each in self.labels: + each.destroy() + for each in self.entries: + each.destroy() + # making new widgets + self.labels = [] + self.entries = [] + for each in self.parameters.keys(): + self.labels.append(Label(self.window, text=each)) + self.entries.append(Entry(self.window, width=15)) + if each == "passband" or each == "passband_deconvolve": + self.entries[-1].insert(0, str(self.parameters[each][0]) + "," + str(self.parameters[each][1])) + else: + self.entries[-1].insert(0, str(self.parameters[each])) + if each == "lag" or each == "dt": + self.entries[-1].configure(state=DISABLED) + row_i = 0 + while row_i < len(self.labels): + key = self.labels[row_i].cget("text") + if key == "estimation" or key == "temporal_mask": + row_i += 1 + continue + self.labels[row_i].grid(row=row_i,column=0,padx=(5,5),pady=(5,5)) + self.entries[row_i].grid(row=row_i,column=1,padx=(5,5),pady=(5,5)) + row_i += 1 diff --git a/build/lib/rsHRF/rsHRF_GUI/gui_windows/plotterOptionsWindow.py b/build/lib/rsHRF/rsHRF_GUI/gui_windows/plotterOptionsWindow.py new file mode 100644 index 0000000..31a2c22 --- /dev/null +++ b/build/lib/rsHRF/rsHRF_GUI/gui_windows/plotterOptionsWindow.py @@ -0,0 +1,76 @@ +from .plotterWindow import PlotterWindow +from tkinter import Toplevel, OptionMenu, StringVar, _setit, Entry, IntVar, Checkbutton + +class PlotterOptionsWindow(): + def __init__(self): + window = Toplevel() + window.title("Plotting Menu") + # get screen width and height + screen_width = window.winfo_screenwidth() + screen_height = window.winfo_screenheight() + window.geometry("650x200+%d+%d" % (screen_width/2.95, 100 + screen_height/2)) + self.plotterScreen = PlotterWindow() + + # can have 3 different plots at once + self.numberOfPlots = self.plotterScreen.get_numberOfPlots() + self.plot = [StringVar() for i in range(self.numberOfPlots)] + self.plotVal = [IntVar() for i in range(self.numberOfPlots)] + self.plotValStore = [0 for i in range(self.numberOfPlots)] + self.plotables = [] + self.options = ["None"] + for each in self.plot: + each.set("None") + + def plotTS(): + for i in range(0, self.numberOfPlots): + if self.plotVal[i].get() != self.plotValStore[i]: + self.plotValStore[i] = self.plotVal[i].get() + if self.plotValStore[i] == 0: + self.plotterScreen.makePlot(None, 0, i) + return + for each in self.plotables: + if each[0] == self.plot[i].get(): + self.plotterScreen.makePlot(each[1][:,int(self.selectVoxel[i].get())], 1, i) + return + self.plotterScreen.makePlot(None, 0, i) + return + + self.plotSelectDropDown = [OptionMenu(window, self.plot[i], *self.options) for i in range(self.numberOfPlots) for i in range(self.numberOfPlots)] + self.selectVoxel = [Entry(window, width=10) for i in range(self.numberOfPlots) for i in range(self.numberOfPlots)] + self.plotButton = [Checkbutton(window, text=" Plot TS " + str(i), variable=self.plotVal[i], command=plotTS) for i in range(self.numberOfPlots)] + + for each in self.selectVoxel: + each.insert(0, 0) + + for i in range(0, self.numberOfPlots): + self.plotSelectDropDown[i].grid(row=0,column=i,padx=(5,5),pady=(5,5)) + self.selectVoxel[i].grid (row=1,column=i,padx=(5,5),pady=(5,5)) + self.plotButton[i].grid (row=2,column=i,padx=(5,5),pady=(5,5)) + + def updatePlots(self, plotables): + temp_names = [] + for each in self.plotables: + temp_names.append(each[0]) + for each in plotables: + if each[0] not in temp_names: + self.plotables.append(each) + self.options = ["None"] + for each in self.plotables: + self.options.append(each[0]) + self.options[1:] = sorted(self.options[1:]) + self.updateWidgets() + + def updateWidgets(self): + + if len(self.options) == 0: + self.options.append("None") + + for each in self.plotSelectDropDown: + each['menu'].delete(0, 'end') + + for choice in self.options: + for i in range(0, self.numberOfPlots): + self.plotSelectDropDown[i]['menu'].add_command(label=choice, command=_setit(self.plot[i], choice)) + + def get_plotables(self): + return self.plotables diff --git a/build/lib/rsHRF/rsHRF_GUI/gui_windows/plotterWindow.py b/build/lib/rsHRF/rsHRF_GUI/gui_windows/plotterWindow.py new file mode 100644 index 0000000..793be61 --- /dev/null +++ b/build/lib/rsHRF/rsHRF_GUI/gui_windows/plotterWindow.py @@ -0,0 +1,37 @@ +import mpld3 +import matplotlib +matplotlib.use("TkAgg") +import numpy as np +from matplotlib.figure import Figure +from tkinter import ttk, Toplevel, Canvas, TOP, BOTH, BOTTOM +from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg, NavigationToolbar2Tk + +class PlotterWindow(): + def __init__(self): + window = Toplevel() + window.title("Screen") + # get screen width and height + screen_width = window.winfo_screenwidth() + screen_height = window.winfo_screenheight() + window.geometry("600x400+%d+%d" % (screen_width-600, screen_height/2)) + figure = Figure(figsize=(5,5), dpi=100) + self.numberOfPlots = 3 + self.ts = [[] for i in range(0, self.numberOfPlots)] + self.plot = [figure.add_subplot(111) for i in range(0, self.numberOfPlots)] + self.canvas = FigureCanvasTkAgg(figure, window) + + def get_numberOfPlots(self): + return self.numberOfPlots + + def makePlot(self, ts, val, num): + for each in self.plot: + each.clear() + if val == 0: + self.ts[num] = 0 + elif val == 1: + self.ts[num] = ts + for i in range(self.numberOfPlots): + self.plot[i].plot(self.ts[i]) + self.canvas.draw() + self.canvas.get_tk_widget().pack(side=BOTTOM, fill=BOTH, expand=True) + diff --git a/build/lib/rsHRF/rsHRF_GUI/misc/__init__.py b/build/lib/rsHRF/rsHRF_GUI/misc/__init__.py new file mode 100644 index 0000000..9ff8118 --- /dev/null +++ b/build/lib/rsHRF/rsHRF_GUI/misc/__init__.py @@ -0,0 +1 @@ +from . import status \ No newline at end of file diff --git a/build/lib/rsHRF/rsHRF_GUI/misc/status.py b/build/lib/rsHRF/rsHRF_GUI/misc/status.py new file mode 100644 index 0000000..e5b9769 --- /dev/null +++ b/build/lib/rsHRF/rsHRF_GUI/misc/status.py @@ -0,0 +1,39 @@ +class Status(): + def __init__(self, state, error="", info="", time="", dic={}): + self.state = state + self.info = info + self.error = error + self.time = time + self.dic = dic + + # setters + def set_state(self, b): + self.state = b + + def set_info(self, s): + self.info = s + + def set_error(self, e): + self.error = e + + def set_time(self, t): + self.time = t + + def set_dic(self, d): + self.dic = d + + # getters + def get_state(self): + return self.state + + def get_info(self): + return self.info + + def get_error(self): + return self.error + + def get_time(self): + return self.time + + def get_dic(self): + return self.dic \ No newline at end of file diff --git a/build/lib/rsHRF/rsHRF_GUI/run.py b/build/lib/rsHRF/rsHRF_GUI/run.py new file mode 100644 index 0000000..06d3ebe --- /dev/null +++ b/build/lib/rsHRF/rsHRF_GUI/run.py @@ -0,0 +1,4 @@ +from .gui_windows.main import Main + +def run(): + Main() \ No newline at end of file diff --git a/build/lib/rsHRF/sFIR/__init__.py b/build/lib/rsHRF/sFIR/__init__.py new file mode 100644 index 0000000..3c91074 --- /dev/null +++ b/build/lib/rsHRF/sFIR/__init__.py @@ -0,0 +1 @@ +from . import smooth_fir diff --git a/build/lib/rsHRF/sFIR/smooth_fir.py b/build/lib/rsHRF/sFIR/smooth_fir.py new file mode 100644 index 0000000..88e1ced --- /dev/null +++ b/build/lib/rsHRF/sFIR/smooth_fir.py @@ -0,0 +1,161 @@ +import numpy as np +from scipy import linalg +from ..processing import knee + +import warnings +warnings.filterwarnings("ignore") + +def wgr_regress(y, X): + n, ncolX = X.shape + Q,R,perm = linalg.qr(X, mode='economic', pivoting=True) + if R.ndim == 0: + p = 0 + elif R.ndim == 1: + p = int(abs(R[0]) > 0) + else: + if np.amin(R.shape) == 1: + p = int(abs(R[0]) > 0) + else: + p = np.sum(np.abs(np.diagonal(R)) > abs(max(n, ncolX)*np.spacing(R[0][0]))) + if p < ncolX: + R = R[0:p,0:p] + Q = Q[:,0:p] + perm = perm[0:p] + b = np.zeros((ncolX)) + if(R.shape[0] == R.shape[1]): + try: + b[perm] = linalg.solve(R,np.matmul(Q.T,y)) + except: + b[perm] = linalg.lstsq(R,np.matmul(Q.T,y)) + else: + b[perm] = linalg.lstsq(R,np.matmul(Q.T,y)) + return b + +def wgr_glsco(X, Y, sMRI = [], AR_lag=0, max_iter=20): + """ + Linear regression when disturbance terms follow AR(p) + ----------------------------------- + Model: + Yt = Xt * Beta + ut , + ut = Phi1 * u(t-1) + ... + Phip * u(t-p) + et + where et ~ N(0,s^2) + ----------------------------------- + Algorithm: + Cochrane-Orcutt iterated regression (Feasible generalized least squares) + ----------------------------------- + Usage: + Y = dependent variable (n * 1 vector) + X = regressors (n * k matrix) + AR_lag = number of lags in AR process + ----------------------------------- + Returns: + Beta = estimator corresponding to the k regressors + """ + nobs, nvar = X.shape + if sMRI == []: + Beta = wgr_regress(Y,X) + else: + sMRI = np.array(sMRI) + try: + Beta = linalg.solve(np.matmul(X.T,X)+sMRI,np.matmul(X.T,Y)) + except: + Beta = linalg.lstsq(np.matmul(X.T,X)+sMRI,np.matmul(X.T,Y)) + resid = Y - np.matmul(X,Beta) + if AR_lag == 0: + res_sum = np.cov(resid) + return res_sum, Beta + max_tol = min(1e-6, max(np.absolute(Beta)) / 1000) + for r in range(max_iter): + Beta_temp = Beta + X_AR = np.zeros((nobs - (2 * AR_lag), AR_lag)) + for m in range(AR_lag): + X_AR[:, m] = resid[AR_lag - m - 1:nobs - AR_lag - m - 1] + Y_AR = resid[AR_lag:nobs - AR_lag] + AR_para = wgr_regress(Y_AR, X_AR) + X_main = X[AR_lag:nobs, :] + Y_main = Y[AR_lag:nobs] + for m in range(AR_lag): + X_main = \ + X_main - (AR_para[m] * (X[AR_lag - m - 1:nobs - m - 1, :])) + Y_main = Y_main - (AR_para[m] * (Y[AR_lag - m - 1:nobs - m - 1])) + if sMRI == []: + Beta = wgr_regress(Y_main, X_main) + else: + try: + Beta = linalg.solve(np.matmul(X_main.T,X_main)+sMRI,np.matmul(X_main.T,Y_main)) + except: + Beta = linalg.lstsq(np.matmul(X_main.T,X_main)+sMRI,np.matmul(X_main.T,Y_main)) + resid = Y[AR_lag:nobs] - X[AR_lag:nobs, :].dot(Beta) + if(max(np.absolute(Beta - Beta_temp)) < max_tol): + break + res_sum = np.cov(resid) + return res_sum, Beta + +def Fit_sFIR2(output, length, TR, input, T, flag_sfir, AR_lag): + NN = int(np.floor(length/TR)) + _input = np.expand_dims(input[0], axis=0) + X = linalg.toeplitz(input, np.concatenate((_input, np.zeros((1, NN-1))), axis = 1)) + X = np.concatenate((X, np.ones((input.shape))), axis = 1) + if flag_sfir: + fwhm = 7 #fwhm=7 seconds smoothing - ref. Goutte + nh = NN-1 + dt = TR + _ = np.expand_dims(np.arange(1, nh+1).T, axis=1) + C = np.matmul(_,np.ones((1, nh))) + h = np.sqrt(1./(fwhm/dt)) + v = 0.1 + R = v * np.exp(-h/2 * (C-C.T)**2) + RI = linalg.inv(R) + MRI = np.zeros((nh + 1, nh + 1)) + MRI[0:nh,0:nh] = RI; + sigma = 1 + sMRI0 = sigma**2*MRI + sMRI = np.zeros((NN+1, NN+1)) + sMRI[0:NN,0:NN] = sMRI0; + if AR_lag == 0: + try: + hrf = linalg.solve((np.matmul(X.T,X)+sMRI),np.matmul(X.T,output)) + except: + hrf = linalg.lstsq((np.matmul(X.T,X)+sMRI),np.matmul(X.T,output)) + resid = output - np.matmul(X, hrf) + res_sum = np.cov(resid) + else: + res_sum, hrf = wgr_glsco(X,output,AR_lag=AR_lag,sMRI=sMRI); + else: + if AR_lag == 0: + hrf = linalg.lstsq(X,output) + hrf = hrf[0] + resid = output - np.matmul(X, hrf) + res_sum = np.cov(resid) + else: + res_sum, hrf = wgr_glsco(X,output,AR_lag=AR_lag) + return hrf, res_sum + +def wgr_FIR_estimation_HRF(u, dat, para, N): + if para['estimation'] == 'sFIR': + firmode = 1 + else: + firmode = 0 + lag = para['lag'] + nlag = np.amax(lag.shape) + len_bin = int(np.floor(para['len'] / para['TR'])) + hrf = np.zeros((len_bin+1, nlag)) + Cov_E = np.zeros((1, nlag)) + kk = 0 + for i_lag in range(1, nlag + 1): + RR = u - i_lag + RR = RR[RR >= 0] + if RR.size != 0: + design = np.zeros((N, 1)) + design[RR] = 1 + hrf_kk, e3 = Fit_sFIR2(dat, para['len'], para['TR'], design, len_bin, firmode, para['AR_lag']) + hrf[:, kk] = np.ravel(hrf_kk) + Cov_E[:, kk] = (np.ravel(e3)) + else: + Cov_E[:, kk] = np.inf + kk += 1 + placeholder, ind = knee.knee_pt(np.ravel(Cov_E)) + if ind == np.amax(Cov_E.shape) - 1: + ind = ind - 1 + rsH = hrf[:,ind+1] + return rsH, u \ No newline at end of file diff --git a/build/lib/rsHRF/spm_dep/__init__.py b/build/lib/rsHRF/spm_dep/__init__.py new file mode 100644 index 0000000..6d21f0c --- /dev/null +++ b/build/lib/rsHRF/spm_dep/__init__.py @@ -0,0 +1,2 @@ +from . import spm + diff --git a/build/lib/rsHRF/spm_dep/spm.py b/build/lib/rsHRF/spm_dep/spm.py new file mode 100644 index 0000000..c693867 --- /dev/null +++ b/build/lib/rsHRF/spm_dep/spm.py @@ -0,0 +1,153 @@ +import math +import numpy as np +import nibabel as nib +from scipy.special import gammaln + +import warnings +warnings.filterwarnings("ignore") + +def spm_vol(input_file): + """ + Get header information for images + """ + v = nib.load(input_file) + return v + + +def spm_read_vols(mapped_image_volume): + """ + Read in entire image volumes + """ + data = mapped_image_volume.get_fdata() + data = data.flatten(order='F') + return data + + +def spm_orth(X, OPT='pad'): + """ + Recursive Gram-Schmidt orthogonalisation of basis functions + @X - matrix + @OPT - 'norm' - for Euclidean normalisation + 'pad' - for zero padding of null space (default) + """ + def gs_cofficient(v1, v2): + return np.dot(v2, v1) / np.dot(v1, v1) + + def multiply(cofficient, v): + return map((lambda x: x * cofficient), v) + + def proj(v1, v2): + return multiply(gs_cofficient(v1, v2), v1) + + def gs(X, row_vecs=True, norm=True): + if not row_vecs: + X = X.T + Y = X[0:1, :].copy() + for i in range(1, X.shape[0]): + proj = np.diag((X[i, :].dot(Y.T) / + np.linalg.norm(Y, axis=1) ** 2).flat).dot(Y) + Y = np.vstack((Y, X[i, :] - proj.sum(0))) + if norm: + Y = np.diag(1 / np.linalg.norm(Y, axis=1)).dot(Y) + if row_vecs: + return Y + else: + return Y.T + + if OPT == 'norm': + return gs(X, row_vecs=False, norm=True) + elif OPT == 'pad': + return gs(X, row_vecs=False, norm=False) + else: + return X + + +def spm_hrf(RT, P=None, fMRI_T=16): + """ + @RT - scan repeat time + @P - parameters of the response function (two gamma functions) + + defaults (seconds) + % P[0] - Delay of Response (relative to onset) 6 + % P[1] - Delay of Undershoot (relative to onset) 16 + % P[2] - Dispersion of Response 1 + % P[3] - Dispersion of Undershoot 1 + % P[4] - Ratio of Response to Undershoot 6 + % P[5] - Onset (seconds) 0 + % P[6] - Length of Kernel (seconds) 32 + + hrf - hemodynamic response function + P - parameters of the response function + """ + p = np.array([6, 16, 1, 1, 6, 0, 32], dtype=float) + if P is not None: + p[0:len(P)] = P + _spm_Gpdf = lambda x, h, l: \ + np.exp(h * np.log(l) + (h - 1) * np.log(x) - (l * x) - gammaln(h)) + # modelled hemodynamic response function - {mixture of Gammas} + dt = RT / float(fMRI_T) + u = np.arange(0, int(p[6] / dt + 1)) - p[5] / dt + with np.errstate(divide='ignore'): # Known division-by-zero + hrf = _spm_Gpdf( + u, p[0] / p[2], dt / p[2] + ) - _spm_Gpdf( + u, p[1] / p[3], dt / p[3] + ) / p[4] + idx = np.arange(0, int((p[6] / RT) + 1)) * fMRI_T + hrf = hrf[idx] + hrf = np.nan_to_num(hrf) + hrf = hrf / np.sum(hrf) + return hrf + + +def spm_detrend(x, p=0): + """ + Polynomial detrending over columns + + spm_detrend removes linear and nonlinear trends + from column-wise data matrices. + + @x - data matrix + @p - order of polynomial [default : 0] + + Returns: + y - detrended data matrix + """ + m, n = x.shape + if (not m) or (not n): + y = [] + return y + + if (not p): + y = x - np.ones((m, 1), dtype='int') * x.mean(axis=0) + return y + + G = np.zeros((m, p + 1)) + for i in range(0, p + 1): + d = np.arange(1, m + 1) ** i + G[:, i] = d.flatten(1) + y = x - G.dot(np.linalg.pinv(G).dot(x)) + return y + + +def spm_write_vol(image_volume_info, image_voxels, image_name, file_type): + """ + Writes an image volume to disk + + @image_volume_info - a structure containing image volume + information (see spm_vol) + @image_voxels - a one, two or three dimensional matrix + containing the image voxels + @image_name - name of the file to save the image in + """ + if file_type == ".nii" or file_type == ".nii.gz": + data = image_voxels + affine = image_volume_info.affine + image_volume_info = nib.Nifti1Image(data, affine) + nib.save(image_volume_info, image_name + file_type) + else: + file_type = '.gii' + data = image_voxels + gi = nib.GiftiImage() + gi.add_gifti_data_array(nib.gifti.GiftiDataArray(image_voxels)) + nib.gifti.giftiio.write(gi, image_name + file_type) diff --git a/build/lib/rsHRF/unit_tests/__init__.py b/build/lib/rsHRF/unit_tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/build/lib/rsHRF/unit_tests/test_basis_functions.py b/build/lib/rsHRF/unit_tests/test_basis_functions.py new file mode 100644 index 0000000..2b57315 --- /dev/null +++ b/build/lib/rsHRF/unit_tests/test_basis_functions.py @@ -0,0 +1,30 @@ +import pytest +import numpy as np +from ..basis_functions import basis_functions + +def test_fourier_bf(): + assert type(basis_functions.fourier_bf(np.random.random(5), {'estimation': 'fourier', 'order': 3})) == type(np.random.random(1)) + assert type(basis_functions.fourier_bf(np.random.random(5), {'estimation': 'fourier-hanning', 'order': 3})) == type(np.random.random(1)) + assert basis_functions.fourier_bf(np.random.random(5), {'estimation': 'fourier', 'order': 3}).shape == (5, 2*3 + 1) + assert basis_functions.fourier_bf(np.random.random(5), {'estimation': 'fourier-hanning', 'order': 3}).shape == (5, 2*3 + 1) + assert basis_functions.fourier_bf(np.random.random(15), {'estimation': 'fourier', 'order': 3}).shape == (15, 2*3 + 1) + assert basis_functions.fourier_bf(np.random.random(15), {'estimation': 'fourier', 'order': 3}).shape == (15, 2*3 + 1) + assert basis_functions.fourier_bf(np.random.random(7), {'estimation': 'fourier', 'order': 11}).shape == (7, 2*11 + 1) + assert basis_functions.fourier_bf(np.random.random(7), {'estimation': 'fourier-hanning', 'order': 11}).shape == (7, 2*11 + 1) + assert np.allclose(basis_functions.fourier_bf(np.zeros(5), {'estimation': 'fourier', 'order': 3}), np.asarray([[1, 0, 0, 0, 1, 1, 1] for i in range(5)])) + assert np.allclose(basis_functions.fourier_bf(np.zeros(5), {'estimation': 'fourier-hanning', 'order': 3}), np.zeros((5, 2*3+1))) + assert np.allclose(basis_functions.fourier_bf(np.ones(5), {'estimation': 'fourier', 'order': 3}), np.asarray([[1, 0, 0, 0, 1, 1, 1] for i in range(5)])) + assert np.allclose(basis_functions.fourier_bf(np.ones(5), {'estimation': 'fourier-hanning', 'order': 3}), np.zeros((5, 2*3+1))) + assert np.allclose(basis_functions.fourier_bf(np.full(5, np.pi/2), {'estimation': 'fourier', 'order': 3}), np.asarray([[ 1.,-0.43030122,0.77685322,-0.97220684,-0.90268536,0.62968173,-0.23412359] for i in range(5)])) + assert np.allclose(basis_functions.fourier_bf(np.full(5, np.pi/2), {'estimation': 'fourier-hanning', 'order': 3}), np.asarray([[0.9513426809665357,-0.40936391340403033,0.7390536246669112,-0.9249018639367675,-0.8587631122906557,0.599043100699187,-0.222731764045654] for i in range(5)])) + +def test_gamma_bf(): + assert type(basis_functions.gamma_bf(np.random.random(5), 7)) == type(np.random.random(1)) + assert basis_functions.gamma_bf(np.random.random(5), 7).shape == (5, 7) + assert basis_functions.gamma_bf(np.random.random(6), 3).shape == (6, 3) + assert basis_functions.gamma_bf(np.random.random(5), 17).shape == (5, 17) + assert basis_functions.gamma_bf(np.random.random(15), 7).shape == (15, 7) + assert basis_functions.gamma_bf(np.random.random(15), 11).shape == (15, 11) + assert basis_functions.gamma_bf(np.random.random(15), 27).shape == (15, 27) + assert np.allclose(basis_functions.gamma_bf(np.zeros(5), 7), np.zeros((5, 7))) + assert np.allclose(basis_functions.gamma_bf(np.ones(5), 3), np.asarray([[6.13132402e-02, 7.29919526e-05, 2.81323432e-13] for i in range(5)])) \ No newline at end of file diff --git a/build/lib/rsHRF/unit_tests/test_canon_hrf2dd.py b/build/lib/rsHRF/unit_tests/test_canon_hrf2dd.py new file mode 100644 index 0000000..7f8a1ae --- /dev/null +++ b/build/lib/rsHRF/unit_tests/test_canon_hrf2dd.py @@ -0,0 +1,32 @@ +import pytest +import numpy as np +from ..spm_dep import spm +from ..processing import knee +from ..canon import canon_hrf2dd + +def test_wgr_spm_get_canonhrf(): + xBF = {} + xBF['dt'] = 2./3 + xBF['T'] = 3 + xBF['len'] = 24 + col1 = np.asarray([0.0, 0.00044807354048549653, 0.007361556038668112, 0.028700948964288868, 0.06209550607132158, 0.09729254798416208, 0.12429404240995112, 0.13792056661956356, 0.13802824093013008, 0.12762119074617928, 0.11077174845169241, 0.09120490808593994, 0.07165360815082951, 0.05376914935191879, 0.03832112614220298, 0.025486073112810325, 0.015110211890879629, 0.006900111913950147, 0.0005371690913944248, -0.004268194929392522, -0.007761293488920797, -0.010148007538333283, -0.01160630142209315, -0.012296736371437677, -0.012368815040058252, -0.011962915647414259, -0.011209052391768252, -0.01022411801820082, -0.009109051781949583, -0.007946887485859182, -0.006802115673078029, -0.005721375663317213, -0.004735225783873632, -0.003860617757643374, -0.0031036899586120515, -0.0024625525460218694, -0.001929827998394537]) + col2 = np.asarray([0.0, 0.00044807354048549653, 0.007342074519043911, 0.026270427564126324, 0.04604777855227284, 0.05298027426688717, 0.04436283618465861, 0.02599351455862868, 0.005549510765577248, -0.011472257205405234, -0.022704308365117443, -0.02817568454942533, -0.029165400859453156, -0.02728041542804467, -0.023923908293136474, -0.020091968757198677, -0.016373874758523064, -0.013045686593496723, -0.010184308584631133, -0.007764223130961229, -0.005723757639472168, -0.004001906734121734, -0.002552994436344191, -0.001347742979076045, -0.0003679965853120599, 0.0003999224194857678, 0.0009704879108325976, 0.0013618737912804435, 0.0015965904821198564, 0.0017004983374237365, 0.0017009756059719753, 0.001624910779748503, 0.0014969777484681833, 0.0013384308028707623, 0.0011664765672502776, 0.000994160094124114, 0.0008306404084625725]) + col3 = np.asarray([0.0, -0.0033744701883765727, -0.02853758656830626, -0.059581285703262654, -0.06276731227887467, -0.03315584298773516, 0.010290875704026259, 0.046611800411483495, 0.06497283018250533, 0.06515064496825096, 0.05287141888607183, 0.035134036642002486, 0.017473063697658053, 0.0031007110271825955, -0.006822656698791951, -0.01249576703581251, -0.014792986331505636, -0.014769241817926427, -0.01337626141896856, -0.011343705230548555, -0.009163530885765957, -0.0071259520361853165, -0.005372467154512138, -0.003946877599547011, -0.0028361208240172686, -0.001999202787460186, -0.0013856669839915764, -0.0009461276756354484, -0.000637391149634764, -0.00042422350787132274, -0.00027925231676790424, -0.0001819802281565093, -0.00011749866181071447, -7.521986433059284e-05, -4.777425670014074e-05, -3.0120010414058304e-05, -1.8859316315220087e-05]) + xBF['TD_DD'] = 0 + out = canon_hrf2dd.wgr_spm_get_canonhrf(xBF) + assert type(out) == type(np.asarray([])) + assert out.shape == (37, 1) + assert np.allclose(out[:,0], col1) + xBF['TD_DD'] = 1 + out = canon_hrf2dd.wgr_spm_get_canonhrf(xBF) + assert type(out) == type(np.asarray([])) + assert out.shape == (37, 2) + assert np.allclose(out[:,0], col1) + assert np.allclose(out[:,1], col2) + xBF['TD_DD'] = 2 + out = canon_hrf2dd.wgr_spm_get_canonhrf(xBF) + assert type(out) == type(np.asarray([])) + assert out.shape == (37, 3) + assert np.allclose(out[:,0], col1) + assert np.allclose(out[:,1], col2) + assert np.allclose(out[:,2], col3) \ No newline at end of file diff --git a/build/lib/rsHRF/unit_tests/test_iterative_wiener_deconv.py b/build/lib/rsHRF/unit_tests/test_iterative_wiener_deconv.py new file mode 100644 index 0000000..6b9bc49 --- /dev/null +++ b/build/lib/rsHRF/unit_tests/test_iterative_wiener_deconv.py @@ -0,0 +1,15 @@ +import pytest +import numpy as np +from .. import iterative_wiener_deconv + +def test_rsHRF_iterative_wiener_deconv(): + var1 = np.random.randint(100, 150) + var2 = np.random.randint(4, 20) + out = iterative_wiener_deconv.rsHRF_iterative_wiener_deconv(np.random.random(var1), np.random.random(var2)) + assert type(out) == type(np.asarray([])) + assert out.size == var1 + y = np.asarray([0.6589422147608651, 0.7239697210706654, 0.5596745029686809, 0.1518281619871923, 0.9344253850406739, 0.06696275606106217, 0.8730982497140573, 0.22794714268930805, 0.6120490212897929]) + h = np.asarray([0.8401094233417159, 0.5039895222403991, 0.008901275117447982, 0.28477784598041767]) + out_expected = np.asarray([0.02218669, 0.02221875, 0.02217417, 0.02219945, 0.02219854, 0.02217766, 0.02221443, 0.02217762, 0.02221311]) + out = iterative_wiener_deconv.rsHRF_iterative_wiener_deconv(y, h) + assert np.allclose(out, out_expected) \ No newline at end of file diff --git a/build/lib/rsHRF/unit_tests/test_knee.py b/build/lib/rsHRF/unit_tests/test_knee.py new file mode 100644 index 0000000..df3d5be --- /dev/null +++ b/build/lib/rsHRF/unit_tests/test_knee.py @@ -0,0 +1,20 @@ +import pytest +import numpy as np +from ..processing import knee + +def test_knee_pt(): + size = np.random.randint(10, 99) + y = np.random.random(size) + np.random.random(size) * 1j + out1, out2 = knee.knee_pt(y) + # assert type(out1) == np.int_ + # assert type(out2) == np.int_ + assert isinstance(out1, (int, np.integer)) + assert isinstance(out2, (int, np.integer)) + y = np.zeros(size) + np.zeros(size) * 1j + out1, out2 = knee.knee_pt(y) + assert out1 == 2 + assert out2 == 1 + y = np.asarray([(0.13638490363348788+0.4483173781686369j), (0.21565052592551615+0.5219779191597427j), (0.9979379796276104+0.07548834892584189j), (0.29470166952920507+0.36391228981190515j), (0.3464996550382653+0.0014059490581040945j), (0.4641546533051284+0.06567864084007502j), (0.46315096641592435+0.3323511950388086j), (0.17130338335642903+0.9791635674939458j), (0.9625869976661646+0.5830153856154539j), (0.8588617152267485+0.6576125931014645j)]) + out1, out2 = knee.knee_pt(y) + assert out1 == 6 + assert out2 == 0 \ No newline at end of file diff --git a/build/lib/rsHRF/unit_tests/test_parameters.py b/build/lib/rsHRF/unit_tests/test_parameters.py new file mode 100644 index 0000000..6380311 --- /dev/null +++ b/build/lib/rsHRF/unit_tests/test_parameters.py @@ -0,0 +1,11 @@ +import pytest +import numpy as np +from .. import parameters + +def test_wgr_get_parameters(): + assert type(parameters.wgr_get_parameters(np.random.random(np.random.randint(1, 100)), np.random.uniform(0, 5))) == type(np.asarray([])) + assert parameters.wgr_get_parameters(np.random.random(np.random.randint(1, 100)), np.random.uniform(0, 5)).size == 3 + assert np.allclose(parameters.wgr_get_parameters(np.zeros(np.random.randint(1, 100)), np.random.uniform(0, 5)), np.zeros(3)) + dt = np.random.uniform(0, 5) + assert np.allclose(parameters.wgr_get_parameters(np.ones(np.random.randint(100)), dt), np.asarray([1., dt, dt])) + assert np.allclose(parameters.wgr_get_parameters(np.asarray([0.54353434, 0.39763409, 0.92579636, 0.74797229, 0.89733085]), 0.1), np.asarray([0.92579636, 0.3, 0.1])) \ No newline at end of file diff --git a/build/lib/rsHRF/unit_tests/test_rest_filter.py b/build/lib/rsHRF/unit_tests/test_rest_filter.py new file mode 100644 index 0000000..fa1ddd4 --- /dev/null +++ b/build/lib/rsHRF/unit_tests/test_rest_filter.py @@ -0,0 +1,27 @@ +import pytest +import numpy as np +from ..processing import rest_filter + +def test_conn_filter(): + TR = np.random.randint(1, 10) + filter = [0.01, 0.08] + x = np.zeros((152, 1)) + y = rest_filter.conn_filter(TR, filter, x) + assert type(y) == type(np.asarray([])) + assert y.shape == (152, 1) + assert np.allclose(y, np.zeros((152, 1))) + x = np.ones((152, 1)) + y = rest_filter.conn_filter(TR, filter, x) + assert np.allclose(y, np.zeros((152, 1))) + +def test_rest_IdealFilter(): + TR = 2.0 + filter = [0.01, 0.08] + x = np.zeros((152, 1)) + y = rest_filter.rest_IdealFilter(x, TR, filter) + assert type(y) == type(np.asarray([])) + assert y.shape == (152, 1) + assert np.allclose(y, np.zeros((152, 1))) + x = np.ones((152, 1)) + y = rest_filter.rest_IdealFilter(x, TR, filter) + assert np.allclose(y, np.ones((152, 1))) \ No newline at end of file diff --git a/build/lib/rsHRF/unit_tests/test_smooth_fir.py b/build/lib/rsHRF/unit_tests/test_smooth_fir.py new file mode 100644 index 0000000..2d1ead8 --- /dev/null +++ b/build/lib/rsHRF/unit_tests/test_smooth_fir.py @@ -0,0 +1,91 @@ +import pytest +import numpy as np +from ..sFIR import smooth_fir + +def test_wgr_regress(): + val1 = np.random.randint(1, 10) + val2 = np.random.randint(11, 200) + y = np.random.random((val2)) + X = np.random.random((val2, val1)) + out = smooth_fir.wgr_regress(y, X) + assert type(out) == type(np.asarray([])) + assert out.size == val1 + y = np.zeros((val2)) + out = smooth_fir.wgr_regress(y, X) + assert np.allclose(np.zeros((val1)), out) + X = np.asarray([[0.24355643043258524, 0.7455782849363849, 0.23335172385642733, 0.020660847834777285], [0.3072110437035208, 0.18602592674503415, 0.9244063171012394, 0.4804649391334549], [0.8016594604103529, 0.9591356084764378, 0.21194228781653557, 0.9827060442384642], [0.07815744716811968, 0.3166325884280309, 0.11251053277597933, 0.6210284344051261], [0.2646152373376658, 0.5240016539430198, 0.31930431720842334, 0.0251519420757097], [0.021911974219001706, 0.01690869083226454, 0.3297262058981726, 0.6618903359199471], [0.08355709035078951, 0.07514938226594636, 0.0031208479890078022, 0.03888221191825747], [0.9970940046680711, 0.5878412240919779, 0.6558184846102507, 0.7206928654146247], [0.727649859981568, 0.4065500141583437, 0.775884796806493, 0.7015181953666869], [0.10971518494894683, 0.9959080243481756, 0.5422490052771556, 0.8266212359977652], [0.26459232860949, 0.8356742095986958, 0.1095077682766954, 0.14378156713215096], [0.259763773901589, 0.3746508316758117, 0.9001388308688932, 0.9050451624557473], [0.9856637279719391, 0.6479702871193095, 0.09198170835828423, 0.8476385625410208], [0.3815258418918812, 0.6703358792038151, 0.4383431460974705, 0.4641501816990993], [0.8473061691948293, 0.5224300782979008, 0.4121441061796459, 0.7523270716310405]]) + y = np.asarray([0.028931607771739087, 0.06696412129540685, 0.22387070667513875, 0.1646337105002752, 0.8412800842253353, 0.9425084870509369, 0.17822863038093084, 0.8479179507498892, 0.5029580038596141, 0.6444064372391561, 0.4070548187562978, 0.505922721410633, 0.15972752609496743, 0.6658366133135427, 0.8404498444201265]) + out_exp = np.asarray([0.0411757, 0.24767687, 0.42055364, 0.19608055]) + out = smooth_fir.wgr_regress(y, X) + assert(np.allclose(out, out_exp)) + +def test_wgr_glsco(): + val1 = np.random.randint(1, 10) + val2 = np.random.randint(11, 200) + Y = np.random.random((val2)) + X = np.random.random((val2, val1)) + sMRI = [] + AR_lag = 1 + max_iter = 20 + out = smooth_fir.wgr_glsco(X, Y, sMRI = sMRI, AR_lag=AR_lag, max_iter=max_iter) + assert type(out) == type(()) + out1, out2 = out + assert type(out1) == type(np.asarray([])) + assert type(out2) == type(np.asarray([])) + assert out1.size == 1 + assert out2.size == val1 + X = np.asarray([[0.43544989084462715, 0.43777641933208156, 0.643080155843435, 0.813958593825052], [0.38928213850325455, 0.7593771015592158, 0.9254380082419151, 0.9447717588006423], [0.7259162100876503, 0.9598502412184825, 0.2899975081150954, 0.061478720535798836], [0.59879813056311, 0.48257959397766126, 0.5172554453824507, 0.16492819946945403], [0.11836248792536241, 0.5160265578209018, 0.8461462946573081, 0.912137637714782], [0.8737499577883742, 0.03795309267342717, 0.606763078803996, 0.5976806041129074], [0.4718986859311063, 0.2344092091943103, 0.473912352830851, 0.23718309987498087], [0.11576771623753401, 0.9128385253518523, 0.28409971172402615, 0.11940757702240845], [0.16412366113275922, 0.38190328702582255, 0.8469984086100405, 0.6973147802496206], [0.49506282363232845, 0.5190719545916632, 0.14244310294990092, 0.6104011894845355]]) + Y = np.asarray([0.23860649900343323, 0.5973580998661552, 0.3928186546630218, 0.05778943990961405, 0.7063978648709636, 0.07370285332712434, 0.009135396063910783, 0.09269627999757779, 0.5547467325191129, 0.9656335840727327]) + out = smooth_fir.wgr_glsco(X, Y, sMRI = sMRI, AR_lag=AR_lag, max_iter=max_iter) + out1, out2 = out + out1_exp = np.asarray([0.02611677]) + out2_exp = np.asarray([ 0.0367933 , 0.51386258, -0.97340772, 1.26010343]) + assert np.allclose(np.asarray(out1), out1_exp) + assert np.allclose(np.asarray(out2), out2_exp) + +def test_Fit_sFIR2(): + val = np.random.randint(11, 200) + output = np.random.random((val)) + length = 24 + TR = 2.0 + inp = np.random.random((val, 1)) + T = 12 + flag_sfir = 1 + AR_lag = 1 + out = smooth_fir.Fit_sFIR2(output, length, TR, inp, T, flag_sfir, AR_lag) + assert type(out) == type(()) + out1, out2 = out + assert type(out1) == type(np.asarray([])) + assert out1.size == 13 + assert type(out2) == type(np.asarray([])) + assert out2.size == 1 + output = np.zeros((val)) + out = smooth_fir.Fit_sFIR2(output, length, TR, inp, T, flag_sfir, AR_lag) + out1, out2 = out + assert np.allclose(np.zeros((13)), out1) + assert np.allclose(np.zeros((1)), out2) + inp = np.asarray([[0.11208123437880424], [0.47337000532832263], [0.5496760145257539], [0.1984332652066324], [0.29117413061557396], [0.707023774370774], [0.7601515226800849], [0.7330200743158585], [0.5033685350197652], [0.6139854490565899], [0.050357274201077495], [0.15618758068425953], [0.17767321816121273], [0.8491434293211506], [0.9017948328311318], [0.3170850576159929], [0.2645298774321433], [0.3485067071439355], [0.13392947284086754], [0.5564549094193428], [0.7549219830345042], [0.41491153064393316], [0.052470190827242136], [0.7863378666190465], [0.6221018221283174], [0.5318620796466], [0.9184682835509385], [0.09511349772823907], [0.2641100933731456], [0.5577093763583515]]) + output = np.asarray([0.03613177914943355, 0.19932562473088167, 0.3728048908567443, 0.07152471033316932, 0.13231120596801615, 0.8336233605946677, 0.19626102779888965, 0.11014151567005881, 0.2165833771773754, 0.09860670113620307, 0.13462733526174642, 0.4026949077980243, 0.7309591700260996, 0.5494233508094024, 0.36049825129091106, 0.07427512868812536, 0.9117468354810448, 0.1486779741324824, 0.23574106866255196, 0.20385094450271402, 0.17652149436966758, 0.8576716535577852, 0.7281154108083755, 0.5313725666429985, 0.6759241704468735, 0.19017576791553958, 0.6072763968829056, 0.35461323171658643, 0.07176955460623347, 0.3522955937509432]) + out2_exp = 0.062454251905728904 + out1_exp = np.asarray([-0.006013666009659839, -0.021447182383702836, -0.023945735843657407, -0.018371231986189136, -0.02233574736823547, -0.020078382177641022, 0.010920928033064338, 0.051249693019481325, 0.06888676660465594, 0.06398005690442186, 0.048868321199953546, 0.07298304109849973, 0.3096652201182815]) + out = smooth_fir.Fit_sFIR2(output, length, TR, inp, T, flag_sfir, AR_lag) + out1, out2 = out + assert np.allclose(out1_exp ,out1) + assert np.allclose(out2_exp ,out2) + +def test_wgr_FIR_estimation_HRF(): + u = np.random.random((7)) + dat = np.random.random((152)) + para = {'estimation': 'sFIR', 'passband': [0.01, 0.08], 'passband_deconvolve': [0.0, 1.7976931348623157e+308], 'TR': 2.0, 'T': 1, 'T0': 3, 'TD_DD': 2, 'AR_lag': 1, 'thr': np.asarray([ 1., np.inf]), 'temporal_mask': [], 'order': 3, 'len': 24, 'min_onset_search': 4, 'max_onset_search': 8, 'localK': 1, 'wiener': False, 'dt': 0.6666666666666666, 'lag': np.asarray([ 6, 7, 8, 9, 10, 11, 12])} + N = 152 + output = smooth_fir.wgr_FIR_estimation_HRF(u, dat, para, N) + assert type(output) == type(()) + u = np.asarray([0.8265743871886244, 0.4571816356384941, 0.9167817425355507, 0.6497771799779721, 0.6873170490384577, 0.45801903592817417, 0.7458612352409498]) + dat = np.asarray([0.07965334399271984, 0.5784720482033168, 0.7241693158799498, 0.8594543092961181, 0.3333887323361263, 0.8663470970895705, 0.1692672051067885, 0.8132782708370856, 0.3545430163067911, 0.1389518058293373, 0.04008377738610536, 0.5251678917354534, 0.03331965019371441, 0.5560024201801771, 0.6492800171038505]) + output = smooth_fir.wgr_FIR_estimation_HRF(u, dat, para, N) + out1, out2 = output + assert type(out1) == type(np.asarray([])) + assert type(out2) == type(np.asarray([])) + assert np.allclose(out1, np.zeros((13))) + out2_exp = np.array([0.82657439, 0.45718164, 0.91678174, 0.64977718, 0.68731705, 0.45801904, 0.74586124]) + assert np.allclose(out2_exp ,out2) \ No newline at end of file diff --git a/build/lib/rsHRF/unit_tests/test_spm.py b/build/lib/rsHRF/unit_tests/test_spm.py new file mode 100644 index 0000000..7308fc9 --- /dev/null +++ b/build/lib/rsHRF/unit_tests/test_spm.py @@ -0,0 +1,97 @@ +import pytest +from unittest import mock +import os +import math +import numpy as np +import nibabel as nib +from scipy.special import gammaln +from ..spm_dep import spm + +SHAPE = (10, 10, 10, 10) + +def get_data(image_type): + data = np.array(np.random.random(SHAPE), dtype=np.float32) + mask = np.random.random(SHAPE[:3]) > 0.1 + if len(SHAPE) > 3: + data[mask, :] = 0 + else: + data[mask] = 0 + if image_type == 'nifti': + data = nib.Nifti1Image(data, np.eye(4)) + else: + data = nib.gifti.GiftiDataArray(data) + return data + +def test_spm_vol(): + test_file_1 = 'test.gii' + test_file_2 = 'test.gii.gz' + test_file_3 = 'test.nii' + test_file_4 = 'test.nii.gz' + test_files = [test_file_1, test_file_2, test_file_3, test_file_4] + with mock.patch('nibabel.load') as load_mock: + for test_file in test_files: + if 'nii' in test_file: + load_mock.return_value = get_data('nifti') + elif 'gii' in test_file: + load_mock.return_value = get_data('gifti') + v = spm.spm_vol(test_file) + assert ('nii' in test_file or 'gii' in test_file) + if 'nii' in test_file: + assert type(v) == type(nib.Nifti1Image(np.asarray([]), np.eye(4))) + elif 'gii' in test_file: + assert type(v) == type(nib.gifti.GiftiDataArray(np.asarray([]))) + +def test_spm_read_vols(): + nifti = get_data('nifti') + data = spm.spm_read_vols(nifti) + assert type(data) == type(np.asarray([])) + assert data.shape[0] == pow(10, 4) + +def test_spm_orth(): + tests = [(3, 4), (7, 5), (4, 12), (13, 6), (11, 11)] + for test in tests: + X = np.random.random(test) + Y = spm.spm_orth(X) + assert type(Y) == type(X) + assert Y.shape == X.shape + +def test_spm_hrf(): + tests = [.5, 1, 2, 3, 4, 1.5, 2.5, 10] + for test in tests: + hrf = spm.spm_hrf(test) + assert type(hrf) == type(np.asarray([])) + assert len(hrf.shape) == 1 + assert hrf.size in [int(33/test) - 1, int(33/test), int(33/test) + 1] + +def test_spm_detrend(): + tests = [(3, 4), (7, 5), (4, 12), (13, 6), (11, 11)] + for test in tests: + X = np.random.random(test) + Y = spm.spm_detrend(X) + assert type(Y) == type(X) + assert Y.shape == X.shape + Y = Y.T + Y_sum = np.sum(Y, axis=1) + assert np.allclose(Y_sum, np.zeros(Y_sum.shape)) + +def test_spm_write_vol(): + test_file_1 = 'test.gii' + test_file_2 = 'test.gii.gz' + test_file_3 = 'test.nii' + test_file_4 = 'test.nii.gz' + test_files = [test_file_1, test_file_2, test_file_3, test_file_4] + with mock.patch('nibabel.load') as load_mock: + for test_file in test_files: + if 'nii' in test_file: + load_mock.return_value = get_data('nifti') + elif 'gii' in test_file: + load_mock.return_value = get_data('gifti') + v1 = spm.spm_vol(test_file) + mask_data = np.zeros(SHAPE[:-1]).flatten(order='F') + fname = test_file.split('.')[0] + file_type = '.' + test_file.split('.', 1)[1] + spm.spm_write_vol(v1, mask_data, fname, file_type) + if 'gii' in file_type: + file_type = '.gii' + assert os.path.isfile(fname + file_type) + os.remove(fname + file_type) \ No newline at end of file diff --git a/build/lib/rsHRF/utils/__init__.py b/build/lib/rsHRF/utils/__init__.py new file mode 100644 index 0000000..afcda41 --- /dev/null +++ b/build/lib/rsHRF/utils/__init__.py @@ -0,0 +1,2 @@ +from . import hrf_estimation +from . import bids \ No newline at end of file diff --git a/build/lib/rsHRF/utils/bids.py b/build/lib/rsHRF/utils/bids.py new file mode 100644 index 0000000..41ca950 --- /dev/null +++ b/build/lib/rsHRF/utils/bids.py @@ -0,0 +1,53 @@ +"""Utilities to handle BIDS inputs.""" +import os +import sys +import json +from pathlib import Path + +def write_derivative_description(bids_dir, deriv_dir): + from ..__about__ import __version__, DOWNLOAD_URL + + print(__version__) + bids_dir = Path(bids_dir) + deriv_dir = Path(deriv_dir) + desc = { + 'Name': 'rsHRF - retrieve the haemodynamic response function from resting state fMRI data', + 'BIDSVersion': '1.4.0', + 'DatasetType': 'derivative', + 'GeneratedBy': [{ + 'Name': 'rsHRF', + 'Version': __version__, + 'CodeURL': DOWNLOAD_URL, + }], + 'HowToAcknowledge': + 'Please cite our paper : https://doi.org/10.1016/j.neuroimage.2021.118591 ' + , + } + + # Keys that can only be set by environment + if 'RSHRF_DOCKER_TAG' in os.environ: + desc['GeneratedBy'][0]['Container'] = { + "Type": "docker", + "Tag": f"bids/rshrf:{os.environ['FMRIPREP_DOCKER_TAG']}" + } + if 'RSHRF_SINGULARITY_URL' in os.environ: + desc['GeneratedBy'][0]['Container'] = { + "Type": "singularity", + "URI": os.getenv('FMRIPREP_SINGULARITY_URL') + } + + # Keys deriving from source dataset + orig_desc = {} + fname = bids_dir / 'dataset_description.json' + if fname.exists(): + orig_desc = json.loads(fname.read_text()) + + if 'DatasetDOI' in orig_desc: + desc['SourceDatasets'] = [{ + 'URL': f'https://doi.org/{orig_desc["DatasetDOI"]}', + 'DOI': orig_desc['DatasetDOI'] + }] + if 'License' in orig_desc: + desc['License'] = orig_desc['License'] + + Path.write_text(deriv_dir / 'dataset_description.json', json.dumps(desc, indent=4)) diff --git a/build/lib/rsHRF/utils/hrf_estimation.py b/build/lib/rsHRF/utils/hrf_estimation.py new file mode 100644 index 0000000..6b05587 --- /dev/null +++ b/build/lib/rsHRF/utils/hrf_estimation.py @@ -0,0 +1,144 @@ +import os +import shutil +import tempfile +import numpy as np +from scipy import stats +from scipy.sparse import lil_matrix +from joblib import load, dump +from joblib import Parallel, delayed +from rsHRF import processing, sFIR +from ..processing import knee + +import warnings +warnings.filterwarnings("ignore") + +""" +HRF ESTIMATION +""" + +def compute_hrf(bold_sig, para, temporal_mask, p_jobs, bf = None): + para['temporal_mask'] = temporal_mask + N, nvar = bold_sig.shape + folder = tempfile.mkdtemp() + data_folder = os.path.join(folder, 'data') + dump(bold_sig, data_folder) + data = load(data_folder, mmap_mode='r') + results = Parallel(n_jobs=p_jobs)(delayed(estimate_hrf)(data, i, para, + N, bf) for i in range(nvar)) + beta_hrf, event_bold = zip(*results) + try: + shutil.rmtree(folder) + except: + print("Failed to delete: " + folder) + # return np.array(beta_hrf).T, np.array(event_bold) + return np.array(beta_hrf).T, event_bold + # events_bold can have shapes + # (6,), (12,), (2,), (5,), (11,), (8,), (14,), (4,), (7,), (10,), (13,), (3,), (9,) + + +def estimate_hrf(bold_sig, i, para, N, bf = None): + """ + Estimate HRF + """ + dat = bold_sig[:, i] + localK = para['localK'] + if para['estimation'] == 'sFIR' or para['estimation'] == 'FIR': + #Estimate HRF for the sFIR or FIR basis functions + if np.count_nonzero(para['thr']) == 1: + para['thr'] = np.array([para['thr'], np.inf]) + thr = para['thr'] #Thr is a vector for (s)FIR + u = wgr_BOLD_event_vector(N, dat, thr, localK, para['temporal_mask']) + u = u.toarray().flatten('C').ravel().nonzero()[0] + beta_hrf, event_bold = sFIR.smooth_fir.wgr_FIR_estimation_HRF(u, dat, para, N) + else: + thr = [para['thr']] #Thr is a scalar for the basis functions + u0 = wgr_BOLD_event_vector(N, dat, thr, localK, para['temporal_mask']) + u = np.append(u0.toarray(), np.zeros((para['T'] - 1, N)), axis=0) + u = np.reshape(u, (1, - 1), order='F') + beta_hrf = wgr_hrf_fit(dat, para, u, bf) + u = u0.toarray()[0].nonzero()[0] + return beta_hrf, u + +def wgr_onset_design(u, bf, T, T0, nscans): + """ + @u - BOLD event vector (microtime). + @bf - basis set matrix + @T - microtime resolution (number of time bins per scan) + @T0 - microtime onset (reference time bin, see slice timing) + """ + ind = np.arange(0, max(u.shape)) + X = np.empty((0, len(ind))) + for p in range(bf.shape[1]): + x = np.convolve(u, bf[:, p]) + x = x[ind] + X = np.append(X, [x], axis=0) + X = X.T + """ + Resample regressors at acquisition times + """ + X = X[(np.arange(0, nscans) * T) + (T0 - 1), :] + return X + + +def wgr_glm_estimation(dat, u, bf, T, T0, AR_lag): + """ + @u - BOLD event vector (microtime). + """ + nscans = dat.shape[0] + x = wgr_onset_design(u, bf, T, T0, nscans) + X = np.append(x, np.ones((nscans, 1)), axis=1) + res_sum, Beta = sFIR.smooth_fir.wgr_glsco(X, dat, AR_lag=AR_lag) + return np.real(res_sum), Beta + +def wgr_hrf_fit(dat, xBF, u, bf): + """ + @u - BOLD event vector (microtime). + @nlag - time lag from neural event to BOLD event + """ + lag = xBF['lag'] + AR_lag = xBF['AR_lag'] + nlag = len(lag) + erm = np.zeros((1, nlag)) + beta = np.zeros((bf.shape[1] + 1, nlag)) + for i in range(nlag): + u_lag = np.append(u[0, lag[i]:], np.zeros((1, lag[i]))).T + erm[0, i], beta[:, i] = \ + wgr_glm_estimation(dat, u_lag, bf, xBF['T'], xBF['T0'], AR_lag) + x, idx = knee.knee_pt(np.ravel(erm)) + if idx == nlag-1: + idx = idx - 1 + beta_hrf = beta[:, idx+1] + beta_hrf = np.append(beta_hrf, lag[idx+1]) + return beta_hrf + +def wgr_BOLD_event_vector(N, matrix, thr, k, temporal_mask): + """ + Detect BOLD event. + event > thr & event < 3.1 + """ + data = lil_matrix((1, N)) + matrix = matrix[:, np.newaxis] + matrix = np.nan_to_num(matrix) + if 0 in np.array(temporal_mask).shape: + matrix = stats.zscore(matrix, ddof=1) + for t in range(1 + k, N - k + 1): + if matrix[t - 1, 0] > thr[0] and \ + np.all(matrix[t - k - 1:t - 1, 0] < matrix[t - 1, 0]) and \ + np.all(matrix[t - 1, 0] > matrix[t:t + k, 0]): + data[0, t - 1] = 1 + else: + tmp = temporal_mask + for i in range(len(temporal_mask)): + if tmp[i] == 1: + temporal_mask[i] = i + datm = np.mean(matrix[temporal_mask]) + datstd = np.std(matrix[temporal_mask]) + if datstd == 0: datstd = 1 + matrix = (matrix - datm)/datstd + for t in range(1 + k, N - k + 1): + if tmp[t-1]: + if matrix[t - 1, 0] > thr[0] and \ + np.all(matrix[t - k - 1:t - 1, 0] < matrix[t - 1, 0]) and \ + np.all(matrix[t - 1, 0] > matrix[t:t + k, 0]): + data[0, t - 1] = 1. + return data diff --git a/consistency_test/VoxelTest/rsHRF_python.py b/consistency_test/VoxelTest/rsHRF_python.py index 6bc968e..21458a8 100644 --- a/consistency_test/VoxelTest/rsHRF_python.py +++ b/consistency_test/VoxelTest/rsHRF_python.py @@ -66,7 +66,7 @@ v1 = spm_dep.spm.spm_vol(input_file) brain = spm_dep.spm.spm_read_vols(v) voxel_ind = np.where(brain > 0)[0] - data = v1.get_data() + data = v1.get_fdata() nobs = data.shape[3] data1 = np.reshape(data, (-1, nobs), order='F').T data1 = data1[:, voxel_ind] diff --git a/requirements.txt b/requirements.txt index 1eae084..7ba566e 100644 --- a/requirements.txt +++ b/requirements.txt @@ -8,4 +8,4 @@ pandas patsy mpld3 joblib -pyyawt \ No newline at end of file +PyWavelets \ No newline at end of file diff --git a/rsHRF/fourD_rsHRF.py b/rsHRF/fourD_rsHRF.py index 2491f08..1acc1ff 100755 --- a/rsHRF/fourD_rsHRF.py +++ b/rsHRF/fourD_rsHRF.py @@ -49,13 +49,13 @@ def demo_rsHRF(input_file, mask_file, output_dir, para, p_jobs, file_type=".nii" raise ValueError ('Inconsistency in input-mask dimensions' + '\n\tinput_file == ' + name + file_type + '\n\tmask_file == ' + mask_name + file_type) else: if file_type == ".nii" or file_type == ".nii.gz" : - data = v1.get_data() + data = v1.get_fdata() else: data = v1.agg_data() else: print('No atlas provided! Generating mask file...') if file_type == ".nii" or file_type == ".nii.gz" : - data = v1.get_data() + data = v1.get_fdata() brain = np.nanvar(data.reshape(-1, data.shape[3]), -1, ddof=0) else: data = v1.agg_data() @@ -171,7 +171,7 @@ def demo_rsHRF(input_file, mask_file, output_dir, para, p_jobs, file_type=".nii" break pos += 1 event_plot = lil_matrix((1, nobs)) - if event_bold.size: + if len(event_bold)>0:#event_bold.size: event_plot[:, event_bold[pos]] = 1 else: print("No Events Detected!") diff --git a/rsHRF/iterative_wiener_deconv.py b/rsHRF/iterative_wiener_deconv.py index 577a689..774ed7a 100755 --- a/rsHRF/iterative_wiener_deconv.py +++ b/rsHRF/iterative_wiener_deconv.py @@ -1,36 +1,39 @@ -import pyyawt +# import pyyawt +import pywt import numpy as np -from rsHRF.processing import knee +from rsHRF.processing import knee def rsHRF_iterative_wiener_deconv(y, h, Iterations=1000): - N = y.shape[0] - nh = max(h.shape) - h = np.append(h, np.zeros((N - nh, 1))) - H = np.fft.fft(h, axis=0) - Y = np.fft.fft(y, axis=0) - [c ,l] = pyyawt.wavedec(abs(y), 1, 'db2') - sigma = pyyawt.wnoisest(c, l, 1) - Phh = np.square(abs(H)) - sqrdtempnorm = ((((np.linalg.norm(y-np.mean(y), 2)**2) - (N-1)*(sigma**2)))/(np.linalg.norm(h,1))**2) - Nf = (sigma**2)*N - tempreg = Nf/sqrdtempnorm - Pxx0 = np.square(abs(np.multiply(Y, (np.divide(np.conj(H), (Phh + N*tempreg)))))) - Pxx = Pxx0 - Sf = Pxx.reshape(-1, 1) + N = y.shape[0] + nh = max(h.shape) + h = np.append(h, np.zeros((N - nh, 1))) + H = np.fft.fft(h, axis=0) + Y = np.fft.fft(y, axis=0) + coeffs = pywt.wavedec(abs(y), 'db2', level=1) + sigma = np.median(np.abs(coeffs[-1])) / 0.6745 + Phh = np.square(abs(H)) + sqrdtempnorm = ((((np.linalg.norm(y-np.mean(y), 2)**2) - (N-1)*(sigma**2))) / (np.linalg.norm(h,1))**2) + Nf = (sigma**2)*N + tempreg = Nf/sqrdtempnorm + Pxx0 = np.square(abs(np.multiply(Y, (np.divide(np.conj(H), (Phh + N*tempreg)))))) + Pxx = Pxx0 + Sf = Pxx.reshape(-1, 1) for i in range(0, Iterations): - M = np.divide(np.multiply(np.multiply(np.conjugate(H), Pxx), Y), np.add(np.multiply(np.square(abs(H)), Pxx), Nf)) - PxxY = np.divide(np.multiply(Pxx, Nf), np.add(np.multiply(np.square(abs(H)), Pxx), Nf)) - Pxx = np.add(PxxY, np.square(abs(M))) - Sf = np.concatenate((Sf, Pxx.reshape(-1, 1)), axis = 1) - dSf = np.diff(Sf, 1, 1) - dSfmse = np.mean(np.square(dSf), axis=1) - _, idx = knee.knee_pt(dSfmse) - idm = np.argmin(dSfmse) - ratio = np.abs(dSfmse[idx] - dSfmse[idm])/(np.abs(np.max(dSfmse) - np.min(dSfmse))) - if ratio > 0.5: + M = np.divide(np.multiply(np.multiply(np.conjugate(H), Pxx), Y), np.add(np.multiply(np.square(abs(H)), Pxx), Nf)) + PxxY = np.divide(np.multiply(Pxx, Nf), np.add(np.multiply(np.square(abs(H)), Pxx), Nf)) + Pxx = np.add(PxxY, np.square(abs(M))) + Sf = np.concatenate((Sf, Pxx.reshape(-1, 1)), axis=1) + dSf = np.diff(Sf, 1, 1) + dSfmse = np.mean(np.square(dSf), axis=1) + _, idx = knee.knee_pt(dSfmse) + idm = np.argmin(dSfmse) + ratio = np.abs(dSfmse[idx] - dSfmse[idm])/(np.abs(np.max(dSfmse) - np.min(dSfmse))) + if ratio > 0.5: id0 = idm else: id0 = idx - Pxx = Sf[:,id0+1] + # Safe indexing to avoid IndexError + safe_index = min(id0 + 1, Sf.shape[1] - 1) + Pxx = Sf[:, safe_index] WienerFilterEst = np.divide(np.multiply(np.conj(H), Pxx), np.add(np.multiply(np.square(abs(H)), Pxx), Nf)) return np.real(np.fft.ifft(np.multiply(WienerFilterEst, Y))) \ No newline at end of file diff --git a/rsHRF/processing/knee.py b/rsHRF/processing/knee.py index eeac2cb..d1cf95d 100755 --- a/rsHRF/processing/knee.py +++ b/rsHRF/processing/knee.py @@ -19,7 +19,7 @@ def knee_pt_helper(y, x=None): res_x = np.nan idx_of_result = np.nan - if type(y) is not np.ndarray: + if not isinstance(y, np.ndarray): print('knee_pt: y must be a numpy 1D vector') return res_x, idx_of_result else: @@ -32,7 +32,7 @@ def knee_pt_helper(y, x=None): else: if x is None: x_was_none = True - x = np.arange(1, np.amax(y.shape) + 1, dtype=np.int) + x = np.arange(1, np.amax(y.shape) + 1, dtype=int) if x.shape != y.shape: print('knee_pt: y and x must have the same dimensions') return res_x, idx_of_result diff --git a/rsHRF/processing/rest_filter.py b/rsHRF/processing/rest_filter.py index 3d4f177..9cb1e1b 100755 --- a/rsHRF/processing/rest_filter.py +++ b/rsHRF/processing/rest_filter.py @@ -13,7 +13,7 @@ def rest_IdealFilter(x, TR, Bands, m=5000): else: ind_X = [j for j in range((i-1)*m, nvar)] x1 = x[:, ind_X] - x1 = conn_filter(TR,Bands,x1) + np.matlib.repmat(np.mean(x1), x1.shape[0], 1) + x1 = conn_filter(TR,Bands,x1) + np.mean(x1, keepdims=True).repeat(x1.shape[0], axis=0) x[:,ind_X] = x1 return x diff --git a/rsHRF/rsHRF_GUI/core/core.py b/rsHRF/rsHRF_GUI/core/core.py index 3912042..ad48722 100755 --- a/rsHRF/rsHRF_GUI/core/core.py +++ b/rsHRF/rsHRF_GUI/core/core.py @@ -197,7 +197,7 @@ def makeInputFile(self, input_file, mask_file, file_type, mode): # if the time-series is not present if not subject.is_present("BOLD", input_file): if file_type == ".nii" or file_type == ".nii.gz": - data = v1.get_data() + data = v1.get_fdata() nobs = data.shape[3] data1 = np.reshape(data, (-1, nobs), order='F').T elif file_type == ".gii" or file_type == ".gii.gz": diff --git a/rsHRF/spm_dep/spm.py b/rsHRF/spm_dep/spm.py index d2ee2b0..c693867 100755 --- a/rsHRF/spm_dep/spm.py +++ b/rsHRF/spm_dep/spm.py @@ -18,7 +18,7 @@ def spm_read_vols(mapped_image_volume): """ Read in entire image volumes """ - data = mapped_image_volume.get_data() + data = mapped_image_volume.get_fdata() data = data.flatten(order='F') return data diff --git a/rsHRF/unit_tests/test_knee.py b/rsHRF/unit_tests/test_knee.py index 60093af..df3d5be 100644 --- a/rsHRF/unit_tests/test_knee.py +++ b/rsHRF/unit_tests/test_knee.py @@ -6,8 +6,10 @@ def test_knee_pt(): size = np.random.randint(10, 99) y = np.random.random(size) + np.random.random(size) * 1j out1, out2 = knee.knee_pt(y) - assert type(out1) == np.int_ - assert type(out2) == np.int_ + # assert type(out1) == np.int_ + # assert type(out2) == np.int_ + assert isinstance(out1, (int, np.integer)) + assert isinstance(out2, (int, np.integer)) y = np.zeros(size) + np.zeros(size) * 1j out1, out2 = knee.knee_pt(y) assert out1 == 2 diff --git a/rsHRF/utils/hrf_estimation.py b/rsHRF/utils/hrf_estimation.py index e857bf0..6b05587 100755 --- a/rsHRF/utils/hrf_estimation.py +++ b/rsHRF/utils/hrf_estimation.py @@ -30,7 +30,11 @@ def compute_hrf(bold_sig, para, temporal_mask, p_jobs, bf = None): shutil.rmtree(folder) except: print("Failed to delete: " + folder) - return np.array(beta_hrf).T, np.array(event_bold) + # return np.array(beta_hrf).T, np.array(event_bold) + return np.array(beta_hrf).T, event_bold + # events_bold can have shapes + # (6,), (12,), (2,), (5,), (11,), (8,), (14,), (4,), (7,), (10,), (13,), (3,), (9,) + def estimate_hrf(bold_sig, i, para, N, bf = None): """ diff --git a/setup.py b/setup.py index 8673a3b..a742815 100644 --- a/setup.py +++ b/setup.py @@ -45,8 +45,9 @@ def run(self): include_package_data=True, zip_safe=False, python_requires=">=3.6", - install_requires=["numpy", "nibabel", "matplotlib", "scipy", "pybids==0.11.1", "pandas", "patsy", "mpld3", "duecredit", "joblib", "pyyawt"], + install_requires=["numpy", "nibabel", "matplotlib", "scipy", "pybids==0.11.1", "pandas", "patsy", "mpld3", "duecredit", "joblib", "PyWavelets"], cmdclass={ 'verify': VerifyVersionCommand, }, ) +