%% Function to create physiological nuisance regressors % For full documentation of the parameters, see also tapas_physio_new %(e.g., via edit tapas_physio_new) %% possible to do: fix artifacts in respiratory phase caused by wrong sign % of derivative oversample = false; %% Create default parameter structure with all fields physio = tapas_physio_new(); %% Individual Parameter settings. Modify to your need and remove default settings physio.save_dir = {'regress'}; physio.log_files.vendor = 'Siemens_Tics'; physio.log_files.cardiac = {'Physio_example_PULS.log'}; physio.log_files.respiration = {'Physio_example_RESP.log'}; physio.log_files.scan_timing = {'Physio_example_Info.log'}; physio.log_files.relative_start_acquisition = 0; physio.log_files.align_scan = 'last'; physio.scan_timing.sqpar.Nslices = 72; physio.scan_timing.sqpar.TR = 0.8; physio.scan_timing.sqpar.Ndummies = 0; physio.scan_timing.sqpar.Nscans = 420; physio.scan_timing.sqpar.onset_slice = 2;% for single slice % physio.scan_timing.sqpar.onset_slice = [1:9];% for slicewise physio.scan_timing.sync.method = 'scan_timing_log'; physio.preproc.cardiac.modality = 'PPU'; physio.preproc.cardiac.filter.include = false; physio.preproc.cardiac.filter.type = 'butter'; physio.preproc.cardiac.filter.passband = [0.3 9]; physio.preproc.cardiac.initial_cpulse_select.method = 'auto_matched'; physio.preproc.cardiac.initial_cpulse_select.max_heart_rate_bpm = 90; physio.preproc.cardiac.initial_cpulse_select.file = 'initial_cpulse_kRpeakfile.mat'; physio.preproc.cardiac.initial_cpulse_select.min = 0.4; physio.preproc.cardiac.posthoc_cpulse_select.method = 'off'; physio.preproc.cardiac.posthoc_cpulse_select.percentile = 80; physio.preproc.cardiac.posthoc_cpulse_select.upper_thresh = 60; physio.preproc.cardiac.posthoc_cpulse_select.lower_thresh = 60; physio.preproc.respiratory.filter.passband = [0.01 2]; physio.preproc.respiratory.despike = false; physio.model.orthogonalise = 'none'; physio.model.censor_unreliable_recording_intervals = false; physio.model.output_multiple_regressors = 'example.txt'; physio.model.output_physio = 'example.mat'; physio.model.retroicor.include = true; physio.model.retroicor.order.c = 3;%3 physio.model.retroicor.order.r = 4;%4 physio.model.retroicor.order.cr = 1;%1 %% modified: added new values for retroicor to produce 1 volume with voxel- % wise (slicewise) delays and to save the raw phases % physio.model.retroicor.voxelwise = true; physio.model.retroicor.voxelwise_volume = 'saufunctional.nii'; %mask for slicewise regressors physio.model.retroicor.voxelwise_volume_mask = 'saufunctional_brain_mask.nii.gz'; physio.model.retroicor.save_phase = 0; % physio.model.rvt.include = false; physio.model.rvt.method = 'hilbert'; physio.model.rvt.delays = 0; physio.model.hrv.include = false; physio.model.hrv.delays = 0; physio.model.noise_rois.include = false; physio.model.noise_rois.fmri_files = {'swaufunctional.nii'}; physio.model.noise_rois.roi_files = {'white_matter.nii', 'csf.nii'}; physio.model.noise_rois.thresholds = 0.9; physio.model.noise_rois.n_voxel_crop = 0; physio.model.noise_rois.n_components = 1; physio.model.noise_rois.force_coregister = 0; physio.model.movement.include = false; physio.model.movement.order = 6; physio.model.movement.censoring_threshold = 0.5; physio.model.movement.censoring_method = 'FD'; physio.model.other.include = false; physio.verbose.level = 0; physio.verbose.process_log = cell(0, 1); physio.verbose.fig_handles = zeros(1, 0); physio.verbose.use_tabs = true; physio.verbose.show_figs = true; physio.verbose.save_figs = false; physio.verbose.close_figs = false; physio.ons_secs.c_scaling = 1; physio.ons_secs.r_scaling = 1; physio.version = 'R2022a-v8.1.0'; %% Taken and modified from tapas_physio_main_create_regressors % Include subfolders of code to path as well pathThis = fileparts(mfilename('fullpath')); addpath(genpath(pathThis)); % These parameters could become toolbox inputs... minConstantIntervalAlertSeconds = 0.2; % fill up empty parameters physio = tapas_physio_fill_empty_parameters(physio); % replace cellstrings physio = tapas_physio_cell2char(physio); % prepend absolute directories - save_dir physio = tapas_physio_prepend_absolute_paths(physio); % set sub-structures for readability; NOTE: copy by value, physio-structure % not updated! ons_secs = physio.ons_secs; save_dir = physio.save_dir; log_files = physio.log_files; preproc = physio.preproc; scan_timing = physio.scan_timing; model = physio.model; verbose = physio.verbose; hasPhaseLogfile = strcmpi(log_files.vendor, 'CustomPhase'); doesNeedPhyslogFiles = model.retroicor.include || model.rvt.include || model.hrv.include; hasPhyslogFiles = ~isempty(log_files.cardiac) || ~isempty(log_files.respiration); if ~hasPhaseLogfile % read and preprocess logfiles only, if model-based physiological noise correction is needed if doesNeedPhyslogFiles if ~hasPhyslogFiles verbose = tapas_physio_log(['No physlog files specified, but models relying on ' ... 'physiological recordings selected. I will skip those.'], ... verbose, 1); sqpar = scan_timing.sqpar; else %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% 1. Read in vendor-specific physiological log-files %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [ons_secs.c, ons_secs.r, ons_secs.t, ons_secs.cpulse, ons_secs.acq_codes, ... verbose] = tapas_physio_read_physlogfiles(... log_files, preproc.cardiac.modality, verbose, scan_timing.sqpar); % also: normalize cardiac/respiratory data, if wanted doNormalize = true; % Normalize and pad time series after read-In ons_secs = tapas_physio_preprocess_phys_timeseries(ons_secs, ... scan_timing.sqpar, doNormalize); hasCardiacData = ~isempty(ons_secs.c); hasRespData = ~isempty(ons_secs.r); verbose = tapas_physio_plot_raw_physdata(ons_secs, verbose); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% 2. Create scan timing nominally or from logfile % nominal: using entered sequence parameters (nSlices, nScans etc) % Philips: via gradient time-course or existing acq_codes in logfile % GE: nominal % Siemens: from tics (Release VD/E), from .resp/.ecg files (Release VB) % Biopac: using triggers from Digital input (mat file) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [ons_secs, VOLLOCS, LOCS, verbose] = tapas_physio_create_scan_timing(... log_files, scan_timing, ons_secs, verbose); minConstantIntervalAlertSamples = ceil(minConstantIntervalAlertSeconds/ons_secs.dt); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% 3. Extract and preprocess physiological data, crop to scan aquisition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if hasCardiacData % preproc.cardiac.modality = 'OXY'; % 'ECG' or 'OXY' (for pulse oximetry) %% initial pulse select via load from logfile or autocorrelation with 1 %% cardiac pulse [ons_secs.c, verbose] = tapas_physio_filter_cardiac(... ons_secs.t, ons_secs.c, preproc.cardiac.filter, verbose); switch preproc.cardiac.initial_cpulse_select.method case {'load_from_logfile', ''} % do nothing otherwise % run one of the various cardiac pulse detection algorithms [ons_secs.cpulse, verbose] = ... tapas_physio_get_cardiac_pulses(ons_secs.t, ons_secs.c, ... preproc.cardiac.initial_cpulse_select, ... preproc.cardiac.modality, verbose); end %% post-hoc: hand pick additional cardiac pulses or load from previous %% time switch preproc.cardiac.posthoc_cpulse_select.method case {'manual'} % additional manual fill-in of more missed pulses [ons_secs, outliersHigh, outliersLow, verbose] = ... tapas_physio_correct_cardiac_pulses_manually(ons_secs, ... preproc.cardiac.posthoc_cpulse_select, verbose); case {'load'} hasPosthocLogFile = exist(preproc.cardiac.posthoc_cpulse_select.file, 'file') || ... exist([preproc.cardiac.posthoc_cpulse_select.file '.mat'], 'file'); if hasPosthocLogFile % load or set selection to manual, if no file exists osload = load(preproc.cardiac.posthoc_cpulse_select.file, 'ons_secs'); ons_secs = osload.ons_secs; else [ons_secs, outliersHigh, outliersLow, verbose] = ... tapas_physio_correct_cardiac_pulses_manually(ons_secs,... preproc.cardiac.posthoc_cpulse_select, verbose); end case {'off', ''} end % label constant samples as unreliable (clipping/detachment) [ons_secs.c_is_reliable, ~, verbose] = tapas_physio_detect_constants(ons_secs.c, ... minConstantIntervalAlertSamples, [], verbose); ons_secs.c_is_reliable = 1 - ons_secs.c_is_reliable; end if hasRespData % filter respiratory signal [ons_secs.fr, verbose] = tapas_physio_filter_respiratory( ... ons_secs.r, ons_secs.dt, ... preproc.respiratory.filter.passband, preproc.respiratory.despike, ... doNormalize, verbose); % label constant samples as unreliable (clipping/detachment) [ons_secs.r_is_reliable, ~, verbose] = tapas_physio_detect_constants(ons_secs.fr, ... minConstantIntervalAlertSamples, [], verbose); ons_secs.r_is_reliable = 1 - ons_secs.r_is_reliable; end %% modified if oversample % remove cropping to preserve initial size and resolution ons_secs.raw = ons_secs; %required for verbosity else %original [ons_secs, scan_timing.sqpar, verbose] = tapas_physio_crop_scanphysevents_to_acq_window(... ons_secs, scan_timing.sqpar, verbose); end sqpar = scan_timing.sqpar; if verbose.level >= 2 verbose.fig_handles(end+1) = ... tapas_physio_plot_cropped_phys_to_acqwindow(ons_secs, sqpar); end [verbose, ons_secs.c_outliers_low, ons_secs.c_outliers_high, ... ons_secs.r_hist] = ... tapas_physio_plot_raw_physdata_diagnostics(ons_secs.cpulse, ... ons_secs.r, preproc.cardiac.posthoc_cpulse_select, verbose, ... ons_secs.t, ons_secs.c); end else % does NOT NeedPhyslogFiles sqpar = scan_timing.sqpar; end % doesNeedPhyslogFiles else % Phase data saved in log-file already % Read logged phases into object directly load(log_files.cardiac, 'c_phase_probe_regressors') load(log_files.respiratory, 'r_phase_probe_regressors'); ons_secs.r_sample_phase = r_phase_probe_regressors; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% 4. Create physiological noise model regressors for GLM for all specified % slices %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% onset_slices = reshape(sqpar.onset_slice, 1, []); nOnsetSlices = numel(onset_slices); if nOnsetSlices < 1 error('Please specify an onset slice.'); end %% modified to save retroicor regressors to voxelwise regressor; happens % automatically, if multiple onset_slices are specified: % for each onset_slice a delayed regressor with the respective onset is % calculated and then each point in the slice is set equal to the regressor % at this point a reference volume is needed to extract the correct % dimensions and crucially the number of repeats if nOnsetSlices > 1 if model.retroicor.include % read in example volume to get correct dimensions regressor = niftiread(model.retroicor.voxelwise_volume); % % change data type to double (default is uint16) % regressor = cast(regressor, 'double'); % get number of repeats of slice delays; equal to multiband factor % in multiband acquisitions; total number of slices should be % divisible by number of onset_slices slice_dimensions = size(regressor, 1, 2); slice_repeats = idivide(size(regressor, 3), uint16(nOnsetSlices)); % total number of retroicor regressors retroicor_order = 2*(model.retroicor.order.c + ... model.retroicor.order.r) + 4*model.retroicor.order.cr; % increase order by 2 if phases are included if model.retroicor.save_phase retroicor_order = retroicor_order + 2; end % predefine array for every regressor per slice retroicor_regressors = zeros(size(regressor, 4), retroicor_order,... nOnsetSlices); % save some memory clear regressor end end % for onset_slice = onset_slices model.R = []; model.R_column_names = {}; %% 4.1. Slice specific parameter adaptation sqpar.onset_slice = onset_slice; if hasPhaseLogfile % explicit down-sampling of pre-existing phases % for Field Probe + Physiological Noise Paper ons_secs.c_sample_phase = ... c_phase_probe_regressors((140+sqpar.onset_slice):(sqpar.Nslices):end); ons_secs.r_sample_phase = ... r_phase_probe_regressors((sqpar.onset_slice):(sqpar.Nslices):end); else % otherwise reset, since phases will be estimated from raw % c(ardiac) and r(espiratory) time courses ons_secs.c_sample_phase = []; ons_secs.r_sample_phase = []; end %% Physiological measures if hasPhyslogFiles %% 4.2. Create RETROICOR regressors (Fourier expansion of cardiac/respiratory phase) if model.retroicor.include %% modified for oversampling if oversample [cardiac_sess, respire_sess, mult_sess, ons_secs, ... model.retroicor.order, verbose] = ... tapas_physio_create_retroicor_regressors_mod(ons_secs, sqpar, ... model.retroicor.order, verbose); else [cardiac_sess, respire_sess, mult_sess, ons_secs, ... model.retroicor.order, verbose] = ... tapas_physio_create_retroicor_regressors(ons_secs, sqpar, ... model.retroicor.order, verbose); end % if model.censor_unreliable_recording_intervals [ons_secs, cardiac_sess, respire_sess, mult_sess, verbose] = ... tapas_physio_censor_unreliable_regressor_parts_retroicor(... ons_secs, sqpar, cardiac_sess, respire_sess, mult_sess, verbose); end %% modified to include raw phase if required if model.retroicor.save_phase [model.R, model.R_column_names] = append_regressors(model.R, model.R_column_names, ... ons_secs.c_sample_phase, 'cardiac phase'); [model.R, model.R_column_names] = append_regressors(model.R, model.R_column_names, ... ons_secs.r_sample_phase, 'respiratory phase'); end % [model.R, model.R_column_names] = append_regressors(model.R, model.R_column_names, ... cardiac_sess, 'RETROICOR (cardiac)'); [model.R, model.R_column_names] = append_regressors(model.R, model.R_column_names, ... respire_sess, 'RETROICOR (respiratory)'); [model.R, model.R_column_names] = append_regressors(model.R, model.R_column_names, ... mult_sess, 'RETROICOR (multiplicative)'); %% modified to save voxelwise regressors if nOnsetSlices > 1 retroicor_regressors(:,:,onset_slice) = model.R; end % end %% 4.3. Create a heart-rate variability regressor using the cardiac response % function if model.hrv.include %% modified if oversample [convHRV, ons_secs.hrv, verbose] = tapas_physio_create_hrv_regressors_mod(... ons_secs, sqpar, model.hrv, verbose); else %standard path, also modified! [convHRV, ons_secs.hrv, verbose] = tapas_physio_create_hrv_regressors(... ons_secs, sqpar, model.hrv, verbose); %% modified -> also normalize standard regressors to (-1 ,1) convHRV = convHRV/max(abs(convHRV)); end [model.R, model.R_column_names] = append_regressors(model.R, model.R_column_names, ... convHRV, 'HR * CRF'); end %% 4.4. Create a respiratory volume/time regressor using the respiratory response % function if model.rvt.include %% modified if oversample [convRVT, ons_secs.rvt, verbose] = tapas_physio_create_rvt_regressors_mod(... ons_secs, sqpar, model.rvt, verbose); else %standard path, also modified! [convRVT, ons_secs.rvt, verbose] = tapas_physio_create_rvt_regressors(... ons_secs, sqpar, model.rvt, verbose); %% modified -> also normalize standard regressors to (-1 ,1) convRVT = convRVT/max(abs(convRVT)); % end [model.R, model.R_column_names] = append_regressors(model.R, model.R_column_names, ... convRVT, 'RVT * RRF'); end end % hasPhyslogFiles %% 4.5. Extract anatomical defined (ROI) principal component regressors if model.noise_rois.include [noise_rois_R, model.noise_rois, verbose] = tapas_physio_create_noise_rois_regressors(... model.noise_rois, verbose); [model.R, model.R_column_names] = append_regressors(model.R, model.R_column_names, ... noise_rois_R, 'Noise ROIs'); end %% 4.6. Load other (physiological) confound regressors if model.other.include && ~isempty(model.other.input_multiple_regressors) [other_R, verbose] = tapas_physio_load_other_multiple_regressors(... model.other.input_multiple_regressors, verbose); [model.R, model.R_column_names] = append_regressors(model.R, model.R_column_names, ... other_R, 'Other'); end %% 4.7. Load and manipulate movement parameters as confound regressors if model.movement.include && ~isempty(model.movement.file_realignment_parameters) [movement_R, model.movement, verbose] = ... tapas_physio_create_movement_regressors(model.movement, verbose); [model.R, model.R_column_names] = append_regressors(model.R, model.R_column_names, ... movement_R(:, 1:model.movement.order), 'Movement'); [model.R, model.R_column_names] = append_regressors(model.R, model.R_column_names, ... movement_R(:, model.movement.order+1:end), 'Motion outliers'); end %% 4.8. Orthogonalisation of regressors ensures numerical stability for % otherwise correlated cardiac regressors [model.R, verbose] = tapas_physio_orthogonalise_physiological_regressors(... model.R, model.R_column_names, model.orthogonalise, ... verbose); %% 4.9 Save Multiple Regressors file for SPM physio.save_dir = save_dir; physio.log_files = log_files; physio.preproc = preproc; physio.scan_timing = scan_timing; physio.model = model; physio.verbose = verbose; physio.ons_secs = ons_secs; %% modified to include onset_slice in the final physio structure physio.scan_timing.sqpar = sqpar; % % determine file names for output, append slice index, if multiple slices % chosen if nOnsetSlices > 1 % save final physio-structure in .mat-file [fp, fn, ext] = fileparts(model.output_physio); file_output_physio = ... fullfile(fp, [fn, sprintf('_slice%03d', onset_slice), ext]); [fp, fn, ext] = fileparts(model.output_multiple_regressors); file_output_multiple_regressors = ... fullfile(fp, [fn, sprintf('_slice%03d', onset_slice), ext]); else file_output_physio = model.output_physio; file_output_multiple_regressors = model.output_multiple_regressors; end if ~isempty(model.output_physio) save(file_output_physio, 'physio'); end if isempty(model.R) disp(['No model estimated. Only saving read-in log-files data into physio ' ... 'mat-output-file instead: Check variable physio.ons_secs']); else [fpfx, fn, fsfx] = fileparts(file_output_multiple_regressors); R = model.R; names = model.R_column_names; switch fsfx case '.mat' save(file_output_multiple_regressors, 'R', 'names'); % SPM understands `names` otherwise save(file_output_multiple_regressors, 'R', '-ascii', '-double', '-tabs'); end end end % onset_slices %% modified to save retroicor regressors with multiple onset slices to % one nii volume with slicewise regressors if nOnsetSlices > 1 % get proper header of fmri file nii_info = niftiinfo(model.retroicor.voxelwise_volume); % change data type in header to match the regressor nii_info.Datatype = 'double'; [fp, fn, ~] = fileparts(model.output_multiple_regressors); ext = '.nii'; % save one volume for each retroicor order for order = [1: retroicor_order] % get regressor with wrong dimensions (time as first dim, slices as % third dim) regressor_slice_by_time = retroicor_regressors(:, order, :); % resize regressor such that time is the 4th axis and slices % are the 3rd (matching the output volume) reg_size = size(regressor_slice_by_time); reg_size = [1, 1, reg_size(3), reg_size(1)]; regressor = zeros(reg_size); % need to do resizing through loop because reshape() gives wrong % result (counts timepoints through slices first instead of keeping % axes consistent) for slice = onset_slices regressor(1, 1, slice, :) = regressor_slice_by_time(:, 1, slice); end % fill the entire volume with repeats of the regressor regressor = repmat(regressor, [slice_dimensions, slice_repeats, 1]); % get name of regressor reg_name = model.R_column_names{order}; filename = append(fn, '_'); % append abbreviations to filename to indicate regressor type if contains(reg_name, 'cardiac') filename = append(filename, 'c'); elseif contains(reg_name, 'respiratory') filename = append(filename, 'r'); else filename = append(filename, 'cr'); end % append regressor order as number to filename (0 for phase, % positive integers otherwise) if contains(reg_name, 'phase') filename = append(filename, '0'); else % find the first instance of the current regressor name % (all cardiac/resp regressors get the same name) lowest_order = find(ismember(model.R_column_names, ... reg_name), 1); % subtract first instance from current order and add 1 to get % the proper order of the current regressor type (c/r/cr) filename = ... append(filename, num2str(order - lowest_order + 1)); end % save regressor volume as nii output_filename = fullfile(fp, [filename, ext]); niftiwrite(regressor, output_filename, nii_info); % mask output volume masked_filename = strcat(filename, '_masked'); masked_filename = fullfile(fp, masked_filename); system(strcat("fslmaths ", output_filename, " -mul ", ... model.retroicor.voxelwise_volume_mask, " ", masked_filename)); end % loop over order end % nOnsetslices %end of mod %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% 5. Save output figures to files %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [physio.verbose] = tapas_physio_print_figs_to_file(physio.verbose); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Helper functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [R, names] = append_regressors(R, names, regressors, name) R = [R, regressors]; names = [names, repmat({name}, 1, size(regressors, 2))]; end