From 2bae0af06b287b83498d56e27082118e908a592b Mon Sep 17 00:00:00 2001 From: Lars Kasper Date: Thu, 3 Aug 2017 23:28:18 +0200 Subject: [PATCH 01/14] Added PhysIO as submodule, pointing to tag R2017.2 --- .gitmodules | 3 + PhysIO | 1 + PhysIO/README.txt | 230 --- PhysIO/code/COPYING | 674 -------- .../tapas_physio_cardiac_detect_outliers.m | 120 -- PhysIO/code/tapas_physio_cell2char.m | 38 - PhysIO/code/tapas_physio_cfg_matlabbatch.m | 1482 ----------------- PhysIO/code/tapas_physio_check_efficacy.m | 179 -- ...tapas_physio_check_get_regressor_columns.m | 109 -- .../code/tapas_physio_check_get_xcon_index.m | 37 - .../code/tapas_physio_check_job_contrasts.mat | Bin 40848 -> 0 bytes PhysIO/code/tapas_physio_check_job_report.mat | Bin 1376 -> 0 bytes ...tapas_physio_check_prepare_job_contrasts.m | 65 - PhysIO/code/tapas_physio_compute_tsnr_gains.m | 135 -- PhysIO/code/tapas_physio_compute_tsnr_spm.m | 217 --- PhysIO/code/tapas_physio_corrcoef12.m | 80 - ...s_physio_correct_cardiac_pulses_manually.m | 127 -- .../tapas_physio_count_physio_regressors.m | 65 - .../tapas_physio_create_LOCS_from_VOLLOCS.m | 82 - .../code/tapas_physio_create_hrv_regressors.m | 136 -- ...s_physio_create_missing_physio_contrasts.m | 79 - .../tapas_physio_create_movement_regressors.m | 117 -- ...apas_physio_create_noise_rois_regressors.m | 213 --- ...tapas_physio_create_retroicor_regressors.m | 173 -- .../code/tapas_physio_create_rvt_regressors.m | 140 -- PhysIO/code/tapas_physio_create_scan_timing.m | 115 -- ...physio_create_scan_timing_from_acq_codes.m | 104 -- ...eate_scan_timing_from_end_marker_philips.m | 172 -- ..._scan_timing_from_gradients_auto_philips.m | 479 ------ ...reate_scan_timing_from_gradients_philips.m | 316 ---- ...sio_create_scan_timing_from_tics_siemens.m | 83 - .../tapas_physio_create_scan_timing_nominal.m | 98 -- PhysIO/code/tapas_physio_crf.m | 37 - ...physio_crop_scanphysevents_to_acq_window.m | 134 -- PhysIO/code/tapas_physio_detect_constants.m | 67 - PhysIO/code/tapas_physio_downsample_phase.m | 18 - .../tapas_physio_example_ge_ppu3t_spm_job.m | 39 - .../tapas_physio_example_multi_subjects.m | 83 - ...pas_physio_example_philips_ecg3t_spm_job.m | 46 - ...pas_physio_example_philips_ecg7t_spm_job.m | 42 - ...pas_physio_example_siemens_ecg3t_spm_job.m | 39 - PhysIO/code/tapas_physio_fieldnamesr.m | 177 -- PhysIO/code/tapas_physio_filename2path.m | 41 - .../code/tapas_physio_fill_empty_parameters.m | 63 - PhysIO/code/tapas_physio_filter_respiratory.m | 67 - PhysIO/code/tapas_physio_find_ecg_r_peaks.m | 143 -- ...physio_find_matching_scanphyslog_philips.m | 63 - PhysIO/code/tapas_physio_find_string.m | 66 - PhysIO/code/tapas_physio_findpeaks.m | 206 --- ...as_physio_findpeaks_template_correlation.m | 270 --- .../tapas_physio_findpeaks_template_xcorr.m | 209 --- PhysIO/code/tapas_physio_get_cardiac_phase.m | 104 -- .../tapas_physio_get_cardiac_pulse_template.m | 210 --- PhysIO/code/tapas_physio_get_cardiac_pulses.m | 91 - ...s_physio_get_cardiac_pulses_auto_matched.m | 106 -- ...hysio_get_cardiac_pulses_manual_template.m | 88 - .../tapas_physio_get_contrast_names_default.m | 37 - .../tapas_physio_get_default_fig_params.m | 83 - ...apas_physio_get_filename_from_id_philips.m | 26 - .../code/tapas_physio_get_fourier_expansion.m | 38 - .../code/tapas_physio_get_onsets_from_locs.m | 134 -- .../tapas_physio_get_oxy_pulses_filtered.m | 76 - .../code/tapas_physio_get_respiratory_phase.m | 164 -- PhysIO/code/tapas_physio_get_sample_points.m | 53 - PhysIO/code/tapas_physio_hr.m | 49 - PhysIO/code/tapas_physio_init_philips.m | 71 - PhysIO/code/tapas_physio_job2physio.m | 66 - ...as_physio_load_other_multiple_regressors.m | 42 - PhysIO/code/tapas_physio_log.m | 57 - .../tapas_physio_main_create_regressors.m | 396 ----- PhysIO/code/tapas_physio_maxfilter.m | 50 - PhysIO/code/tapas_physio_new.m | 631 ------- ...o_orthogonalise_physiological_regressors.m | 97 -- PhysIO/code/tapas_physio_overlay_contrasts.m | 197 --- ...as_physio_plot_cropped_phys_to_acqwindow.m | 142 -- PhysIO/code/tapas_physio_plot_gradient.m | 64 - PhysIO/code/tapas_physio_plot_raw_physdata.m | 90 - ...pas_physio_plot_raw_physdata_diagnostics.m | 84 - .../tapas_physio_plot_raw_physdata_siemens.m | 55 - .../tapas_physio_plot_retroicor_regressors.m | 106 -- PhysIO/code/tapas_physio_prctile.m | 43 - .../tapas_physio_prepend_absolute_paths.m | 51 - .../tapas_physio_preprocess_phys_timeseries.m | 99 -- PhysIO/code/tapas_physio_print_figs_to_file.m | 73 - PhysIO/code/tapas_physio_propval.m | 210 --- .../code/tapas_physio_read_physlog_PHILIPS.m | 200 --- PhysIO/code/tapas_physio_read_physlogfiles.m | 102 -- .../code/tapas_physio_read_physlogfiles_GE.m | 138 -- ...apas_physio_read_physlogfiles_biopac_mat.m | 190 --- .../tapas_physio_read_physlogfiles_custom.m | 131 -- .../tapas_physio_read_physlogfiles_philips.m | 146 -- ..._physio_read_physlogfiles_philips_matrix.m | 46 - .../tapas_physio_read_physlogfiles_siemens.m | 248 --- ...pas_physio_read_physlogfiles_siemens_raw.m | 66 - ...as_physio_read_physlogfiles_siemens_resp.m | 460 ----- ...as_physio_read_physlogfiles_siemens_tics.m | 241 --- .../tapas_physio_repair_scan_events_PHILIPS.m | 184 -- .../tapas_physio_replace_absolute_paths.m | 58 - PhysIO/code/tapas_physio_report_contrasts.m | 179 -- ...hysio_rescale_gradient_gain_fluctuations.m | 170 -- PhysIO/code/tapas_physio_review.m | 146 -- PhysIO/code/tapas_physio_rrf.m | 39 - PhysIO/code/tapas_physio_rvt.m | 170 -- .../code/tapas_physio_save_batch_mat_script.m | 115 -- .../tapas_physio_scaleorthmean_regressors.m | 39 - PhysIO/code/tapas_physio_siemens_line2table.m | 118 -- .../code/tapas_physio_siemens_table2cardiac.m | 102 -- ...apas_physio_sort_images_by_cardiac_phase.m | 342 ---- ..._physio_sort_images_by_respiratory_phase.m | 312 ---- .../tapas_physio_split_regressor_slices.m | 43 - PhysIO/code/tapas_physio_strip_fields.m | 44 - PhysIO/code/tapas_physio_uddpvparse.m | 135 -- PhysIO/code/tapas_physio_update_from_job.m | 116 -- PhysIO/manual/QuickStart_PhysIO_Toolbox.pdf | Bin 4730854 -> 0 bytes 114 files changed, 4 insertions(+), 15838 deletions(-) create mode 100644 .gitmodules create mode 160000 PhysIO delete mode 100644 PhysIO/README.txt delete mode 100644 PhysIO/code/COPYING delete mode 100644 PhysIO/code/tapas_physio_cardiac_detect_outliers.m delete mode 100644 PhysIO/code/tapas_physio_cell2char.m delete mode 100644 PhysIO/code/tapas_physio_cfg_matlabbatch.m delete mode 100644 PhysIO/code/tapas_physio_check_efficacy.m delete mode 100644 PhysIO/code/tapas_physio_check_get_regressor_columns.m delete mode 100644 PhysIO/code/tapas_physio_check_get_xcon_index.m delete mode 100644 PhysIO/code/tapas_physio_check_job_contrasts.mat delete mode 100644 PhysIO/code/tapas_physio_check_job_report.mat delete mode 100644 PhysIO/code/tapas_physio_check_prepare_job_contrasts.m delete mode 100644 PhysIO/code/tapas_physio_compute_tsnr_gains.m delete mode 100644 PhysIO/code/tapas_physio_compute_tsnr_spm.m delete mode 100644 PhysIO/code/tapas_physio_corrcoef12.m delete mode 100644 PhysIO/code/tapas_physio_correct_cardiac_pulses_manually.m delete mode 100644 PhysIO/code/tapas_physio_count_physio_regressors.m delete mode 100644 PhysIO/code/tapas_physio_create_LOCS_from_VOLLOCS.m delete mode 100644 PhysIO/code/tapas_physio_create_hrv_regressors.m delete mode 100644 PhysIO/code/tapas_physio_create_missing_physio_contrasts.m delete mode 100644 PhysIO/code/tapas_physio_create_movement_regressors.m delete mode 100644 PhysIO/code/tapas_physio_create_noise_rois_regressors.m delete mode 100644 PhysIO/code/tapas_physio_create_retroicor_regressors.m delete mode 100644 PhysIO/code/tapas_physio_create_rvt_regressors.m delete mode 100644 PhysIO/code/tapas_physio_create_scan_timing.m delete mode 100644 PhysIO/code/tapas_physio_create_scan_timing_from_acq_codes.m delete mode 100644 PhysIO/code/tapas_physio_create_scan_timing_from_end_marker_philips.m delete mode 100644 PhysIO/code/tapas_physio_create_scan_timing_from_gradients_auto_philips.m delete mode 100644 PhysIO/code/tapas_physio_create_scan_timing_from_gradients_philips.m delete mode 100644 PhysIO/code/tapas_physio_create_scan_timing_from_tics_siemens.m delete mode 100644 PhysIO/code/tapas_physio_create_scan_timing_nominal.m delete mode 100644 PhysIO/code/tapas_physio_crf.m delete mode 100644 PhysIO/code/tapas_physio_crop_scanphysevents_to_acq_window.m delete mode 100644 PhysIO/code/tapas_physio_detect_constants.m delete mode 100644 PhysIO/code/tapas_physio_downsample_phase.m delete mode 100644 PhysIO/code/tapas_physio_example_ge_ppu3t_spm_job.m delete mode 100644 PhysIO/code/tapas_physio_example_multi_subjects.m delete mode 100644 PhysIO/code/tapas_physio_example_philips_ecg3t_spm_job.m delete mode 100644 PhysIO/code/tapas_physio_example_philips_ecg7t_spm_job.m delete mode 100644 PhysIO/code/tapas_physio_example_siemens_ecg3t_spm_job.m delete mode 100644 PhysIO/code/tapas_physio_fieldnamesr.m delete mode 100644 PhysIO/code/tapas_physio_filename2path.m delete mode 100644 PhysIO/code/tapas_physio_fill_empty_parameters.m delete mode 100644 PhysIO/code/tapas_physio_filter_respiratory.m delete mode 100644 PhysIO/code/tapas_physio_find_ecg_r_peaks.m delete mode 100644 PhysIO/code/tapas_physio_find_matching_scanphyslog_philips.m delete mode 100644 PhysIO/code/tapas_physio_find_string.m delete mode 100644 PhysIO/code/tapas_physio_findpeaks.m delete mode 100644 PhysIO/code/tapas_physio_findpeaks_template_correlation.m delete mode 100644 PhysIO/code/tapas_physio_findpeaks_template_xcorr.m delete mode 100644 PhysIO/code/tapas_physio_get_cardiac_phase.m delete mode 100644 PhysIO/code/tapas_physio_get_cardiac_pulse_template.m delete mode 100644 PhysIO/code/tapas_physio_get_cardiac_pulses.m delete mode 100644 PhysIO/code/tapas_physio_get_cardiac_pulses_auto_matched.m delete mode 100644 PhysIO/code/tapas_physio_get_cardiac_pulses_manual_template.m delete mode 100644 PhysIO/code/tapas_physio_get_contrast_names_default.m delete mode 100644 PhysIO/code/tapas_physio_get_default_fig_params.m delete mode 100644 PhysIO/code/tapas_physio_get_filename_from_id_philips.m delete mode 100644 PhysIO/code/tapas_physio_get_fourier_expansion.m delete mode 100644 PhysIO/code/tapas_physio_get_onsets_from_locs.m delete mode 100644 PhysIO/code/tapas_physio_get_oxy_pulses_filtered.m delete mode 100644 PhysIO/code/tapas_physio_get_respiratory_phase.m delete mode 100644 PhysIO/code/tapas_physio_get_sample_points.m delete mode 100644 PhysIO/code/tapas_physio_hr.m delete mode 100644 PhysIO/code/tapas_physio_init_philips.m delete mode 100644 PhysIO/code/tapas_physio_job2physio.m delete mode 100644 PhysIO/code/tapas_physio_load_other_multiple_regressors.m delete mode 100644 PhysIO/code/tapas_physio_log.m delete mode 100644 PhysIO/code/tapas_physio_main_create_regressors.m delete mode 100644 PhysIO/code/tapas_physio_maxfilter.m delete mode 100644 PhysIO/code/tapas_physio_new.m delete mode 100644 PhysIO/code/tapas_physio_orthogonalise_physiological_regressors.m delete mode 100644 PhysIO/code/tapas_physio_overlay_contrasts.m delete mode 100644 PhysIO/code/tapas_physio_plot_cropped_phys_to_acqwindow.m delete mode 100644 PhysIO/code/tapas_physio_plot_gradient.m delete mode 100644 PhysIO/code/tapas_physio_plot_raw_physdata.m delete mode 100644 PhysIO/code/tapas_physio_plot_raw_physdata_diagnostics.m delete mode 100644 PhysIO/code/tapas_physio_plot_raw_physdata_siemens.m delete mode 100644 PhysIO/code/tapas_physio_plot_retroicor_regressors.m delete mode 100644 PhysIO/code/tapas_physio_prctile.m delete mode 100644 PhysIO/code/tapas_physio_prepend_absolute_paths.m delete mode 100644 PhysIO/code/tapas_physio_preprocess_phys_timeseries.m delete mode 100644 PhysIO/code/tapas_physio_print_figs_to_file.m delete mode 100644 PhysIO/code/tapas_physio_propval.m delete mode 100644 PhysIO/code/tapas_physio_read_physlog_PHILIPS.m delete mode 100644 PhysIO/code/tapas_physio_read_physlogfiles.m delete mode 100644 PhysIO/code/tapas_physio_read_physlogfiles_GE.m delete mode 100644 PhysIO/code/tapas_physio_read_physlogfiles_biopac_mat.m delete mode 100644 PhysIO/code/tapas_physio_read_physlogfiles_custom.m delete mode 100644 PhysIO/code/tapas_physio_read_physlogfiles_philips.m delete mode 100644 PhysIO/code/tapas_physio_read_physlogfiles_philips_matrix.m delete mode 100644 PhysIO/code/tapas_physio_read_physlogfiles_siemens.m delete mode 100644 PhysIO/code/tapas_physio_read_physlogfiles_siemens_raw.m delete mode 100644 PhysIO/code/tapas_physio_read_physlogfiles_siemens_resp.m delete mode 100644 PhysIO/code/tapas_physio_read_physlogfiles_siemens_tics.m delete mode 100644 PhysIO/code/tapas_physio_repair_scan_events_PHILIPS.m delete mode 100644 PhysIO/code/tapas_physio_replace_absolute_paths.m delete mode 100644 PhysIO/code/tapas_physio_report_contrasts.m delete mode 100644 PhysIO/code/tapas_physio_rescale_gradient_gain_fluctuations.m delete mode 100644 PhysIO/code/tapas_physio_review.m delete mode 100644 PhysIO/code/tapas_physio_rrf.m delete mode 100644 PhysIO/code/tapas_physio_rvt.m delete mode 100644 PhysIO/code/tapas_physio_save_batch_mat_script.m delete mode 100644 PhysIO/code/tapas_physio_scaleorthmean_regressors.m delete mode 100644 PhysIO/code/tapas_physio_siemens_line2table.m delete mode 100644 PhysIO/code/tapas_physio_siemens_table2cardiac.m delete mode 100644 PhysIO/code/tapas_physio_sort_images_by_cardiac_phase.m delete mode 100644 PhysIO/code/tapas_physio_sort_images_by_respiratory_phase.m delete mode 100644 PhysIO/code/tapas_physio_split_regressor_slices.m delete mode 100644 PhysIO/code/tapas_physio_strip_fields.m delete mode 100644 PhysIO/code/tapas_physio_uddpvparse.m delete mode 100644 PhysIO/code/tapas_physio_update_from_job.m delete mode 100644 PhysIO/manual/QuickStart_PhysIO_Toolbox.pdf diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 00000000..aaec2302 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "PhysIO"] + path = PhysIO + url = git@tnurepository.ethz.ch:physio/physio-public.git diff --git a/PhysIO b/PhysIO new file mode 160000 index 00000000..aa1e0da2 --- /dev/null +++ b/PhysIO @@ -0,0 +1 @@ +Subproject commit aa1e0da2889adeb61175280619515e7be144bec9 diff --git a/PhysIO/README.txt b/PhysIO/README.txt deleted file mode 100644 index 901015e1..00000000 --- a/PhysIO/README.txt +++ /dev/null @@ -1,230 +0,0 @@ -TAPAS PhysIO Toolbox Version 2016 - -************************************************************************ -Copyright (C) 2012-2016 Lars Kasper -Translational Neuromodeling Unit (TNU) -Institute for Biomedical Engineering -University of Zurich and ETH Zurich ------------------------------------------------------------------------- - -The PhysIO Toolbox is free software: you can redistribute it and/or -modify it under the terms of the GNU General Public License as -published by the Free Software Foundation, either version 3 of the -License, or (at your option) any later version. - -This program is distributed in the hope that it will be useful, but -WITHOUT ANY WARRANTY; without even the implied warranty of -MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -General Public License for more details. - -You should have received a copy of the GNU General Public License -along with this program (see the file COPYING). If not, see -. -************************************************************************ - -PURPOSE - -The PhysIO Toolbox provides model-based physiological noise correction of -fMRI data using peripheral measures of respiration and cardiac pulsation. -It incorporates noise models of cardiac/respiratory phase (RETROICOR, -Glover et al. 2000), as well as heart rate variability and respiratory -volume per time (cardiac response function, Chang et. al, 2009, respiratory -response function, Birn et al. 2006). The toolbox is usable via the SPM -batch editor, performs automatic pre-processing of noisy peripheral data -and outputs nuisance regressor files directly suitable for SPM -(multiple_regressors.txt). - -BACKGROUND - -The PhysIO Toolbox provides physiological noise correction for fMRI-data -from peripheral measures (ECG/pulse oximetry, breathing belt). It is -model-based, i.e. creates nuisance regressors from the physiological -monitoring that can enter a General Linear Model (GLM) analysis, e.g. -SPM8/12. Furthermore, for scanner vendor logfiles (PHILIPS, GE, Siemens), -it provides means to statistically assess peripheral data (e.g. heart rate variability) -and recover imperfect measures (e.g. distorted R-peaks of the ECG). - -Facts about physiological noise in fMRI: -- Physiological noise can explain 20-60 % of variance in fMRI voxel time - series (Birn2006, Hutton2011, Harvey2008. - - Physiological noise affects a lot of brain regions (s. figure, e.g. - brainstem or OFC), especially next to CSF, arteries (Hutton2011). - - If not accounted for, this is a key factor limiting sensitivity for - effects of interest. -- Physiological noise contributions increase with field strength; they - become a particular concern at and above 3 Tesla (Kasper2009, Hutton2011). -- In resting state fMRI, disregarding physiological noise leads to wrong - connectivity results (Birn2006). - -=> Some kind of physiological noise correction is highly recommended for - every statistical fMRI analysis. - -Model-based correction of physiological noise: -- Physiological noise can be decomposed into periodic time series following - heart rate and breathing cycle. -- The Fourier expansion of cardiac and respiratory phases was introduced as - RETROICOR (RETROspective Image CORrection, Glover2000, - see also Josephs1997). -- These Fourier Terms can enter a General Linear Model (GLM) as nuisance - regressors, analogous to movement parameters. -- As the physiological noise regressors augment the GLM and explain - variance in the time series, they increase sensitivity in all contrasts - of interest. - - -FEATURES OF THIS TOOLBOX - -Physiological Noise Modeling : -- Modeling physiological noise regressors from peripheral data - (breathing belt, ECG, pulse oximeter) -- State of the art RETROICOR cardiac and respiratory phase expansion -- Cardiac response function (Chang et al, 2009) and respiratory response - function (Birn et al. 2006) modelling of heart-rate variability and - respiratory volume per time influence on physiological noise -- Flexible expansion orders to model different contributions of cardiac, - respiratory and interaction terms (see Harvey2008, Hutton2011) -- Automatic creation of nuisance regressors, full integration into standard - GLMs, tested for SPM8/12 ("multiple_regressors.mat") - -Flexible Read-in: -The toolbox is dedicated to seamless integration into a clinical research s -etting and therefore offers correction methods to recover physiological -data from imperfect peripheral measures. -- General Electric logfiles -- Siemens logfiles (synchronization via DICOM time stamp or tics time scale) -- Philips SCANPHYSLOG-file handling, including automatic alignment of scan volume timing and physiological time series - from logged gradient timecourses -- BioPac .mat-files - - -COMPATIBILITY & SUPPORT - -- Matlab Toolbox -- Input: - - Fully integrated to work with physiological logfiles for Philips MR systems (SCANPHYSLOG) - - tested for General Electric (GE) log-files - - implementation for Siemens log-files - - also: interface for 'Custom', i.e. general heart-beat time stamps - & breathing volume time courses from other log formats -- Output: - - Nuisance regressors for mass-univariate statistical analysis with SPM5,8,12 - or as text file for export to any other package - - raw and processed physiological logfile data -- Part of the TNU Software Edition: long term support and ongoing development - - -DOWNLOADS & RELEASE INFORMATION - -- Current Release: - -PhysIO_Toolbox_R2017.1 - -February 19, 2017 - -Minor Release Notes (R2017.1) -* Substantially improved Siemens interface, both for VB/VD and 3T/7T releases -** several bugfixes -** based on extensive user feedback from Berlin and Brisbane -* New functionality tapas_physio_overlay_contrasts.m to display non-physio - contrasts automatically as well - -Major Release Notes (r904 / R2016.1): -- Software version for accepted PhysIO Toolbox Paper: doi:10.1016/j.jneumeth.2016.10.019 -- Tested and expanded versions of examples -- Improved stability by bugfixes and compatibility to Matlab R2016 -- Slice-wise regressor creation -- Detection of constant physiological time series (detachment, clipping) -- Refactoring of report_contrasts and compute_tsnr_gains as standalone functionality -- Improved Read-in capabilities (Siemens respiration data, BioPac .mat) -- Migration from svn (r904) to git (tnurepository) for version control - -Major Release Notes (r835): -- Software version for Toolbox Paper submission -- Noise ROIs modeling -- Extended motion models (24 parameters, Volterra expansion) -- HRV/RVT models with optional multiple delay regressors -- Report_contrasts with automatic contrast generation for all regressor groups -- compute_tsnr_gains for individual physiological regressor groups -- consistent module naming (scan_timing, preproc) -- Visualisation improvement (color schemes, legends) - -Minor Release Notes (r666): -- Compatibility tested for SPM12, small bugfixes Batch Dependencies -- Cleaner Batch Interface with grouped sub-menus (cfg_choice) -- new model: 'none' to just read out physiological raw data and preprocess, - without noise modelling -- Philips: Scan-timing via gradient log now automatized (gradient_log_auto) -- Siemens: Tics-Logfile read-in (proprietary, needs Siemens-agreement) -- All peak detections (cardiac/respiratory) now via auto_matched algorithm -- Adapt plots/saving for Matlab R2014b - -Major Release Notes (r534): -- Read-in of Siemens plain text log files; new example dataset for Siemens -- Speed up and debugging of auto-detection method for noisy cardiac data => new method thresh.cardiac.initial_cpulse_select.method = ???auto_matched??? -- Error handling for temporary breathing belt failures (Eduardo Aponte, TNU Zurich) -- slice-wise regressors can be created by setting sqpar.onset_slice to a index vector of slices - -Major Release Notes (r497): -- SPM matlabbatch GUI implemented (Call via Batch -> SPM -> Tools -> TAPAS PhysIO Toolbox) -- improved, automatic heartbeat detection for noisy ECG now standard for ECG and Pulse oximetry (courtesy of Steffen Bollmann) -- QuickStart-Manual and PhysIO-Background presentation expanded/updated -- job .m/.mat-files created for all example datasets -- bugfixes cpulse-initial-select method-handling (auto/manual/load) - -Major Release Notes (r429): -- Cardiac and Respiratory response function regressors integrated in workflow (heart rate and breathing volume computation) -- Handling of Cardiac and Respiratory Logfiles only -- expanded documentation (Quickstart.pdf and Handbook.pdf) -- read-in of custom log files, e.g. for BrainVoyager peripheral data -- more informative plots and commenting (especially in tapas_physio_new). - -Minor Release Notes (r354): -- computation of heart and breathing rate in Philips/PPU/main_PPU.m -- prefix of functions with tapas_* - -Major Release Notes (r241): -- complete modularization of reading/preprocessing/regressor creation for peripheral physiological data -- manual selection of missed heartbeats in ECG/pulse oximetry (courtesy of Jakob Heinzle) -- support for logfiles from GE scanners (courtesy of Steffen Bollmann, KiSpi Zuerich) -- improved detection of pulse oximetry peaks (courtesy of Steffen Bollmann) -- improved documentation -- consistent function names (prefixed by "physio_") - -NOTE: Your main_ECG/PPU.m etc. scripts from previous versions (<=r159) will not work with this one any more. Please adapt one of the example scripts for your needs (~5 min of work). The main benefit of this version is a complete new variable structure that is more sustainable and makes the code more readable. - - -Lead Programmer: Lars Kasper, TNU & MR-Technology Group, IBT, University & ETH Zurich - -Contributors: -Steffen Bollmann, Children's Hospital Zurich & ETH Zurich -Jakob Heinzle, TNU Zurich -Eduardo Aponte, TNU Zurich - -Send bug reports and suggestions to: kasper@biomed.ee.ethz.ch - - -TUTORIAL - -Run main_ECG3T.m in subdirectory "examples" of the toolbox -See subdirectory "manual" - - -REFERENCES - -Kasper, L., Bollmann, S., Diaconescu, A.O., Hutton, C., Heinzle, J., Iglesias, S., Hauser, T.U., Sebold, M., Manjaly, Z.-M., Pruessmann, K.P., Stephan, K.E., 2016. The PhysIO Toolbox for Modeling Physiological Noise in fMRI Data. Journal of Neuroscience Methods accepted. doi:10.1016/j.jneumeth.2016.10.019 - -Birn, Rasmus M., Jason B. Diamond, Monica A. Smith, and Peter A. Bandettini. 2006. Separating Respiratory-variation-related Fluctuations from Neuronal-activity-related Fluctuations in fMRI. NeuroImage 31 (4) (July 15): 1536?1548. doi:10.1016/j.neuroimage.2006.02.048. - -Glover, G H, T Q Li, and D Ress. 2000. Image-based Method for Retrospective Correction of Physiological Motion Effects in fMRI: RETROICOR. Magnetic Resonance in Medicine: Official Journal of the Society of Magnetic Resonance in Medicine 44 (1) (July): 162(7). doi:10893535. - -Harvey, Ann K., Kyle T.S. Pattinson, Jonathan C.W. Brooks, Stephen D. Mayhew, Mark Jenkinson, and Richard G. Wise. 2008. Brainstem Functional Magnetic Resonance Imaging: Disentangling Signal from Physiological Noise. Journal of Magnetic Resonance Imaging 28 (6): 1337?1344. doi:10.1002/jmri.21623. - -Hutton, C., O. Josephs, J. Stadler, E. Featherstone, A. Reid, O. Speck, J. Bernarding, and N. Weiskopf. 2011. The Impact of Physiological Noise Correction on fMRI at 7 T. NeuroImage 57 (1) (July 1): 101?112. doi:10.1016/j.neuroimage.2011.04.018. - -Josephs, O., Howseman, A.M., Friston, K., Turner, R., 1997. Physiological noise modelling for multi-slice EPI fMRI using SPM. Proceedings of the 5th Annual Meeting of ISMRM, Vancouver, Canada, p. 1682 - -Kasper, Lars, Sarah Marti, S. Johanna Vannesjo, Chloe Hutton, Ray Dolan, Nikolaus Weiskopf, Klaas Enno Stephan, and Klaas Paul Pruessmann. 2009. Cardiac Artefact Correction for Human Brainstem fMRI at 7 Tesla. In Proc. Org. Hum. Brain Mapping 15, 395. San Francisco. - - -VERSION OF THIS FILE -$Id$ diff --git a/PhysIO/code/COPYING b/PhysIO/code/COPYING deleted file mode 100644 index 94a9ed02..00000000 --- a/PhysIO/code/COPYING +++ /dev/null @@ -1,674 +0,0 @@ - GNU GENERAL PUBLIC LICENSE - Version 3, 29 June 2007 - - Copyright (C) 2007 Free Software Foundation, Inc. - Everyone is permitted to copy and distribute verbatim copies - of this license document, but changing it is not allowed. - - Preamble - - The GNU General Public License is a free, copyleft license for -software and other kinds of works. - - The licenses for most software and other practical works are designed -to take away your freedom to share and change the works. By contrast, -the GNU General Public License is intended to guarantee your freedom to -share and change all versions of a program--to make sure it remains free -software for all its users. We, the Free Software Foundation, use the -GNU General Public License for most of our software; it applies also to -any other work released this way by its authors. You can apply it to -your programs, too. - - When we speak of free software, we are referring to freedom, not -price. Our General Public Licenses are designed to make sure that you -have the freedom to distribute copies of free software (and charge for -them if you wish), that you receive source code or can get it if you -want it, that you can change the software or use pieces of it in new -free programs, and that you know you can do these things. - - To protect your rights, we need to prevent others from denying you -these rights or asking you to surrender the rights. Therefore, you have -certain responsibilities if you distribute copies of the software, or if -you modify it: responsibilities to respect the freedom of others. - - For example, if you distribute copies of such a program, whether -gratis or for a fee, you must pass on to the recipients the same -freedoms that you received. You must make sure that they, too, receive -or can get the source code. And you must show them these terms so they -know their rights. - - Developers that use the GNU GPL protect your rights with two steps: -(1) assert copyright on the software, and (2) offer you this License -giving you legal permission to copy, distribute and/or modify it. - - For the developers' and authors' protection, the GPL clearly explains -that there is no warranty for this free software. For both users' and -authors' sake, the GPL requires that modified versions be marked as -changed, so that their problems will not be attributed erroneously to -authors of previous versions. - - Some devices are designed to deny users access to install or run -modified versions of the software inside them, although the manufacturer -can do so. This is fundamentally incompatible with the aim of -protecting users' freedom to change the software. The systematic -pattern of such abuse occurs in the area of products for individuals to -use, which is precisely where it is most unacceptable. Therefore, we -have designed this version of the GPL to prohibit the practice for those -products. If such problems arise substantially in other domains, we -stand ready to extend this provision to those domains in future versions -of the GPL, as needed to protect the freedom of users. - - Finally, every program is threatened constantly by software patents. -States should not allow patents to restrict development and use of -software on general-purpose computers, but in those that do, we wish to -avoid the special danger that patents applied to a free program could -make it effectively proprietary. To prevent this, the GPL assures that -patents cannot be used to render the program non-free. - - The precise terms and conditions for copying, distribution and -modification follow. - - TERMS AND CONDITIONS - - 0. Definitions. - - "This License" refers to version 3 of the GNU General Public License. - - "Copyright" also means copyright-like laws that apply to other kinds of -works, such as semiconductor masks. - - "The Program" refers to any copyrightable work licensed under this -License. Each licensee is addressed as "you". "Licensees" and -"recipients" may be individuals or organizations. - - To "modify" a work means to copy from or adapt all or part of the work -in a fashion requiring copyright permission, other than the making of an -exact copy. The resulting work is called a "modified version" of the -earlier work or a work "based on" the earlier work. - - A "covered work" means either the unmodified Program or a work based -on the Program. - - To "propagate" a work means to do anything with it that, without -permission, would make you directly or secondarily liable for -infringement under applicable copyright law, except executing it on a -computer or modifying a private copy. Propagation includes copying, -distribution (with or without modification), making available to the -public, and in some countries other activities as well. - - To "convey" a work means any kind of propagation that enables other -parties to make or receive copies. Mere interaction with a user through -a computer network, with no transfer of a copy, is not conveying. - - An interactive user interface displays "Appropriate Legal Notices" -to the extent that it includes a convenient and prominently visible -feature that (1) displays an appropriate copyright notice, and (2) -tells the user that there is no warranty for the work (except to the -extent that warranties are provided), that licensees may convey the -work under this License, and how to view a copy of this License. If -the interface presents a list of user commands or options, such as a -menu, a prominent item in the list meets this criterion. - - 1. Source Code. - - The "source code" for a work means the preferred form of the work -for making modifications to it. "Object code" means any non-source -form of a work. - - A "Standard Interface" means an interface that either is an official -standard defined by a recognized standards body, or, in the case of -interfaces specified for a particular programming language, one that -is widely used among developers working in that language. - - The "System Libraries" of an executable work include anything, other -than the work as a whole, that (a) is included in the normal form of -packaging a Major Component, but which is not part of that Major -Component, and (b) serves only to enable use of the work with that -Major Component, or to implement a Standard Interface for which an -implementation is available to the public in source code form. A -"Major Component", in this context, means a major essential component -(kernel, window system, and so on) of the specific operating system -(if any) on which the executable work runs, or a compiler used to -produce the work, or an object code interpreter used to run it. - - The "Corresponding Source" for a work in object code form means all -the source code needed to generate, install, and (for an executable -work) run the object code and to modify the work, including scripts to -control those activities. However, it does not include the work's -System Libraries, or general-purpose tools or generally available free -programs which are used unmodified in performing those activities but -which are not part of the work. For example, Corresponding Source -includes interface definition files associated with source files for -the work, and the source code for shared libraries and dynamically -linked subprograms that the work is specifically designed to require, -such as by intimate data communication or control flow between those -subprograms and other parts of the work. - - The Corresponding Source need not include anything that users -can regenerate automatically from other parts of the Corresponding -Source. - - The Corresponding Source for a work in source code form is that -same work. - - 2. Basic Permissions. - - All rights granted under this License are granted for the term of -copyright on the Program, and are irrevocable provided the stated -conditions are met. This License explicitly affirms your unlimited -permission to run the unmodified Program. The output from running a -covered work is covered by this License only if the output, given its -content, constitutes a covered work. This License acknowledges your -rights of fair use or other equivalent, as provided by copyright law. - - You may make, run and propagate covered works that you do not -convey, without conditions so long as your license otherwise remains -in force. You may convey covered works to others for the sole purpose -of having them make modifications exclusively for you, or provide you -with facilities for running those works, provided that you comply with -the terms of this License in conveying all material for which you do -not control copyright. Those thus making or running the covered works -for you must do so exclusively on your behalf, under your direction -and control, on terms that prohibit them from making any copies of -your copyrighted material outside their relationship with you. - - Conveying under any other circumstances is permitted solely under -the conditions stated below. Sublicensing is not allowed; section 10 -makes it unnecessary. - - 3. Protecting Users' Legal Rights From Anti-Circumvention Law. - - No covered work shall be deemed part of an effective technological -measure under any applicable law fulfilling obligations under article -11 of the WIPO copyright treaty adopted on 20 December 1996, or -similar laws prohibiting or restricting circumvention of such -measures. - - When you convey a covered work, you waive any legal power to forbid -circumvention of technological measures to the extent such circumvention -is effected by exercising rights under this License with respect to -the covered work, and you disclaim any intention to limit operation or -modification of the work as a means of enforcing, against the work's -users, your or third parties' legal rights to forbid circumvention of -technological measures. - - 4. Conveying Verbatim Copies. - - You may convey verbatim copies of the Program's source code as you -receive it, in any medium, provided that you conspicuously and -appropriately publish on each copy an appropriate copyright notice; -keep intact all notices stating that this License and any -non-permissive terms added in accord with section 7 apply to the code; -keep intact all notices of the absence of any warranty; and give all -recipients a copy of this License along with the Program. - - You may charge any price or no price for each copy that you convey, -and you may offer support or warranty protection for a fee. - - 5. Conveying Modified Source Versions. - - You may convey a work based on the Program, or the modifications to -produce it from the Program, in the form of source code under the -terms of section 4, provided that you also meet all of these conditions: - - a) The work must carry prominent notices stating that you modified - it, and giving a relevant date. - - b) The work must carry prominent notices stating that it is - released under this License and any conditions added under section - 7. This requirement modifies the requirement in section 4 to - "keep intact all notices". - - c) You must license the entire work, as a whole, under this - License to anyone who comes into possession of a copy. This - License will therefore apply, along with any applicable section 7 - additional terms, to the whole of the work, and all its parts, - regardless of how they are packaged. This License gives no - permission to license the work in any other way, but it does not - invalidate such permission if you have separately received it. - - d) If the work has interactive user interfaces, each must display - Appropriate Legal Notices; however, if the Program has interactive - interfaces that do not display Appropriate Legal Notices, your - work need not make them do so. - - A compilation of a covered work with other separate and independent -works, which are not by their nature extensions of the covered work, -and which are not combined with it such as to form a larger program, -in or on a volume of a storage or distribution medium, is called an -"aggregate" if the compilation and its resulting copyright are not -used to limit the access or legal rights of the compilation's users -beyond what the individual works permit. Inclusion of a covered work -in an aggregate does not cause this License to apply to the other -parts of the aggregate. - - 6. Conveying Non-Source Forms. - - You may convey a covered work in object code form under the terms -of sections 4 and 5, provided that you also convey the -machine-readable Corresponding Source under the terms of this License, -in one of these ways: - - a) Convey the object code in, or embodied in, a physical product - (including a physical distribution medium), accompanied by the - Corresponding Source fixed on a durable physical medium - customarily used for software interchange. - - b) Convey the object code in, or embodied in, a physical product - (including a physical distribution medium), accompanied by a - written offer, valid for at least three years and valid for as - long as you offer spare parts or customer support for that product - model, to give anyone who possesses the object code either (1) a - copy of the Corresponding Source for all the software in the - product that is covered by this License, on a durable physical - medium customarily used for software interchange, for a price no - more than your reasonable cost of physically performing this - conveying of source, or (2) access to copy the - Corresponding Source from a network server at no charge. - - c) Convey individual copies of the object code with a copy of the - written offer to provide the Corresponding Source. This - alternative is allowed only occasionally and noncommercially, and - only if you received the object code with such an offer, in accord - with subsection 6b. - - d) Convey the object code by offering access from a designated - place (gratis or for a charge), and offer equivalent access to the - Corresponding Source in the same way through the same place at no - further charge. You need not require recipients to copy the - Corresponding Source along with the object code. If the place to - copy the object code is a network server, the Corresponding Source - may be on a different server (operated by you or a third party) - that supports equivalent copying facilities, provided you maintain - clear directions next to the object code saying where to find the - Corresponding Source. Regardless of what server hosts the - Corresponding Source, you remain obligated to ensure that it is - available for as long as needed to satisfy these requirements. - - e) Convey the object code using peer-to-peer transmission, provided - you inform other peers where the object code and Corresponding - Source of the work are being offered to the general public at no - charge under subsection 6d. - - A separable portion of the object code, whose source code is excluded -from the Corresponding Source as a System Library, need not be -included in conveying the object code work. - - A "User Product" is either (1) a "consumer product", which means any -tangible personal property which is normally used for personal, family, -or household purposes, or (2) anything designed or sold for incorporation -into a dwelling. In determining whether a product is a consumer product, -doubtful cases shall be resolved in favor of coverage. For a particular -product received by a particular user, "normally used" refers to a -typical or common use of that class of product, regardless of the status -of the particular user or of the way in which the particular user -actually uses, or expects or is expected to use, the product. A product -is a consumer product regardless of whether the product has substantial -commercial, industrial or non-consumer uses, unless such uses represent -the only significant mode of use of the product. - - "Installation Information" for a User Product means any methods, -procedures, authorization keys, or other information required to install -and execute modified versions of a covered work in that User Product from -a modified version of its Corresponding Source. The information must -suffice to ensure that the continued functioning of the modified object -code is in no case prevented or interfered with solely because -modification has been made. - - If you convey an object code work under this section in, or with, or -specifically for use in, a User Product, and the conveying occurs as -part of a transaction in which the right of possession and use of the -User Product is transferred to the recipient in perpetuity or for a -fixed term (regardless of how the transaction is characterized), the -Corresponding Source conveyed under this section must be accompanied -by the Installation Information. But this requirement does not apply -if neither you nor any third party retains the ability to install -modified object code on the User Product (for example, the work has -been installed in ROM). - - The requirement to provide Installation Information does not include a -requirement to continue to provide support service, warranty, or updates -for a work that has been modified or installed by the recipient, or for -the User Product in which it has been modified or installed. Access to a -network may be denied when the modification itself materially and -adversely affects the operation of the network or violates the rules and -protocols for communication across the network. - - Corresponding Source conveyed, and Installation Information provided, -in accord with this section must be in a format that is publicly -documented (and with an implementation available to the public in -source code form), and must require no special password or key for -unpacking, reading or copying. - - 7. Additional Terms. - - "Additional permissions" are terms that supplement the terms of this -License by making exceptions from one or more of its conditions. -Additional permissions that are applicable to the entire Program shall -be treated as though they were included in this License, to the extent -that they are valid under applicable law. If additional permissions -apply only to part of the Program, that part may be used separately -under those permissions, but the entire Program remains governed by -this License without regard to the additional permissions. - - When you convey a copy of a covered work, you may at your option -remove any additional permissions from that copy, or from any part of -it. (Additional permissions may be written to require their own -removal in certain cases when you modify the work.) You may place -additional permissions on material, added by you to a covered work, -for which you have or can give appropriate copyright permission. - - Notwithstanding any other provision of this License, for material you -add to a covered work, you may (if authorized by the copyright holders of -that material) supplement the terms of this License with terms: - - a) Disclaiming warranty or limiting liability differently from the - terms of sections 15 and 16 of this License; or - - b) Requiring preservation of specified reasonable legal notices or - author attributions in that material or in the Appropriate Legal - Notices displayed by works containing it; or - - c) Prohibiting misrepresentation of the origin of that material, or - requiring that modified versions of such material be marked in - reasonable ways as different from the original version; or - - d) Limiting the use for publicity purposes of names of licensors or - authors of the material; or - - e) Declining to grant rights under trademark law for use of some - trade names, trademarks, or service marks; or - - f) Requiring indemnification of licensors and authors of that - material by anyone who conveys the material (or modified versions of - it) with contractual assumptions of liability to the recipient, for - any liability that these contractual assumptions directly impose on - those licensors and authors. - - All other non-permissive additional terms are considered "further -restrictions" within the meaning of section 10. If the Program as you -received it, or any part of it, contains a notice stating that it is -governed by this License along with a term that is a further -restriction, you may remove that term. If a license document contains -a further restriction but permits relicensing or conveying under this -License, you may add to a covered work material governed by the terms -of that license document, provided that the further restriction does -not survive such relicensing or conveying. - - If you add terms to a covered work in accord with this section, you -must place, in the relevant source files, a statement of the -additional terms that apply to those files, or a notice indicating -where to find the applicable terms. - - Additional terms, permissive or non-permissive, may be stated in the -form of a separately written license, or stated as exceptions; -the above requirements apply either way. - - 8. Termination. - - You may not propagate or modify a covered work except as expressly -provided under this License. Any attempt otherwise to propagate or -modify it is void, and will automatically terminate your rights under -this License (including any patent licenses granted under the third -paragraph of section 11). - - However, if you cease all violation of this License, then your -license from a particular copyright holder is reinstated (a) -provisionally, unless and until the copyright holder explicitly and -finally terminates your license, and (b) permanently, if the copyright -holder fails to notify you of the violation by some reasonable means -prior to 60 days after the cessation. - - Moreover, your license from a particular copyright holder is -reinstated permanently if the copyright holder notifies you of the -violation by some reasonable means, this is the first time you have -received notice of violation of this License (for any work) from that -copyright holder, and you cure the violation prior to 30 days after -your receipt of the notice. - - Termination of your rights under this section does not terminate the -licenses of parties who have received copies or rights from you under -this License. If your rights have been terminated and not permanently -reinstated, you do not qualify to receive new licenses for the same -material under section 10. - - 9. Acceptance Not Required for Having Copies. - - You are not required to accept this License in order to receive or -run a copy of the Program. Ancillary propagation of a covered work -occurring solely as a consequence of using peer-to-peer transmission -to receive a copy likewise does not require acceptance. However, -nothing other than this License grants you permission to propagate or -modify any covered work. These actions infringe copyright if you do -not accept this License. Therefore, by modifying or propagating a -covered work, you indicate your acceptance of this License to do so. - - 10. Automatic Licensing of Downstream Recipients. - - Each time you convey a covered work, the recipient automatically -receives a license from the original licensors, to run, modify and -propagate that work, subject to this License. You are not responsible -for enforcing compliance by third parties with this License. - - An "entity transaction" is a transaction transferring control of an -organization, or substantially all assets of one, or subdividing an -organization, or merging organizations. If propagation of a covered -work results from an entity transaction, each party to that -transaction who receives a copy of the work also receives whatever -licenses to the work the party's predecessor in interest had or could -give under the previous paragraph, plus a right to possession of the -Corresponding Source of the work from the predecessor in interest, if -the predecessor has it or can get it with reasonable efforts. - - You may not impose any further restrictions on the exercise of the -rights granted or affirmed under this License. For example, you may -not impose a license fee, royalty, or other charge for exercise of -rights granted under this License, and you may not initiate litigation -(including a cross-claim or counterclaim in a lawsuit) alleging that -any patent claim is infringed by making, using, selling, offering for -sale, or importing the Program or any portion of it. - - 11. Patents. - - A "contributor" is a copyright holder who authorizes use under this -License of the Program or a work on which the Program is based. The -work thus licensed is called the contributor's "contributor version". - - A contributor's "essential patent claims" are all patent claims -owned or controlled by the contributor, whether already acquired or -hereafter acquired, that would be infringed by some manner, permitted -by this License, of making, using, or selling its contributor version, -but do not include claims that would be infringed only as a -consequence of further modification of the contributor version. For -purposes of this definition, "control" includes the right to grant -patent sublicenses in a manner consistent with the requirements of -this License. - - Each contributor grants you a non-exclusive, worldwide, royalty-free -patent license under the contributor's essential patent claims, to -make, use, sell, offer for sale, import and otherwise run, modify and -propagate the contents of its contributor version. - - In the following three paragraphs, a "patent license" is any express -agreement or commitment, however denominated, not to enforce a patent -(such as an express permission to practice a patent or covenant not to -sue for patent infringement). To "grant" such a patent license to a -party means to make such an agreement or commitment not to enforce a -patent against the party. - - If you convey a covered work, knowingly relying on a patent license, -and the Corresponding Source of the work is not available for anyone -to copy, free of charge and under the terms of this License, through a -publicly available network server or other readily accessible means, -then you must either (1) cause the Corresponding Source to be so -available, or (2) arrange to deprive yourself of the benefit of the -patent license for this particular work, or (3) arrange, in a manner -consistent with the requirements of this License, to extend the patent -license to downstream recipients. "Knowingly relying" means you have -actual knowledge that, but for the patent license, your conveying the -covered work in a country, or your recipient's use of the covered work -in a country, would infringe one or more identifiable patents in that -country that you have reason to believe are valid. - - If, pursuant to or in connection with a single transaction or -arrangement, you convey, or propagate by procuring conveyance of, a -covered work, and grant a patent license to some of the parties -receiving the covered work authorizing them to use, propagate, modify -or convey a specific copy of the covered work, then the patent license -you grant is automatically extended to all recipients of the covered -work and works based on it. - - A patent license is "discriminatory" if it does not include within -the scope of its coverage, prohibits the exercise of, or is -conditioned on the non-exercise of one or more of the rights that are -specifically granted under this License. You may not convey a covered -work if you are a party to an arrangement with a third party that is -in the business of distributing software, under which you make payment -to the third party based on the extent of your activity of conveying -the work, and under which the third party grants, to any of the -parties who would receive the covered work from you, a discriminatory -patent license (a) in connection with copies of the covered work -conveyed by you (or copies made from those copies), or (b) primarily -for and in connection with specific products or compilations that -contain the covered work, unless you entered into that arrangement, -or that patent license was granted, prior to 28 March 2007. - - Nothing in this License shall be construed as excluding or limiting -any implied license or other defenses to infringement that may -otherwise be available to you under applicable patent law. - - 12. No Surrender of Others' Freedom. - - If conditions are imposed on you (whether by court order, agreement or -otherwise) that contradict the conditions of this License, they do not -excuse you from the conditions of this License. If you cannot convey a -covered work so as to satisfy simultaneously your obligations under this -License and any other pertinent obligations, then as a consequence you may -not convey it at all. For example, if you agree to terms that obligate you -to collect a royalty for further conveying from those to whom you convey -the Program, the only way you could satisfy both those terms and this -License would be to refrain entirely from conveying the Program. - - 13. Use with the GNU Affero General Public License. - - Notwithstanding any other provision of this License, you have -permission to link or combine any covered work with a work licensed -under version 3 of the GNU Affero General Public License into a single -combined work, and to convey the resulting work. The terms of this -License will continue to apply to the part which is the covered work, -but the special requirements of the GNU Affero General Public License, -section 13, concerning interaction through a network will apply to the -combination as such. - - 14. Revised Versions of this License. - - The Free Software Foundation may publish revised and/or new versions of -the GNU General Public License from time to time. Such new versions will -be similar in spirit to the present version, but may differ in detail to -address new problems or concerns. - - Each version is given a distinguishing version number. If the -Program specifies that a certain numbered version of the GNU General -Public License "or any later version" applies to it, you have the -option of following the terms and conditions either of that numbered -version or of any later version published by the Free Software -Foundation. If the Program does not specify a version number of the -GNU General Public License, you may choose any version ever published -by the Free Software Foundation. - - If the Program specifies that a proxy can decide which future -versions of the GNU General Public License can be used, that proxy's -public statement of acceptance of a version permanently authorizes you -to choose that version for the Program. - - Later license versions may give you additional or different -permissions. However, no additional obligations are imposed on any -author or copyright holder as a result of your choosing to follow a -later version. - - 15. Disclaimer of Warranty. - - THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY -APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT -HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY -OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, -THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM -IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF -ALL NECESSARY SERVICING, REPAIR OR CORRECTION. - - 16. Limitation of Liability. - - IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING -WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS -THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY -GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE -USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF -DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD -PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), -EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF -SUCH DAMAGES. - - 17. Interpretation of Sections 15 and 16. - - If the disclaimer of warranty and limitation of liability provided -above cannot be given local legal effect according to their terms, -reviewing courts shall apply local law that most closely approximates -an absolute waiver of all civil liability in connection with the -Program, unless a warranty or assumption of liability accompanies a -copy of the Program in return for a fee. - - END OF TERMS AND CONDITIONS - - How to Apply These Terms to Your New Programs - - If you develop a new program, and you want it to be of the greatest -possible use to the public, the best way to achieve this is to make it -free software which everyone can redistribute and change under these terms. - - To do so, attach the following notices to the program. It is safest -to attach them to the start of each source file to most effectively -state the exclusion of warranty; and each file should have at least -the "copyright" line and a pointer to where the full notice is found. - - - Copyright (C) - - This program is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with this program. If not, see . - -Also add information on how to contact you by electronic and paper mail. - - If the program does terminal interaction, make it output a short -notice like this when it starts in an interactive mode: - - Copyright (C) - This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'. - This is free software, and you are welcome to redistribute it - under certain conditions; type `show c' for details. - -The hypothetical commands `show w' and `show c' should show the appropriate -parts of the General Public License. Of course, your program's commands -might be different; for a GUI interface, you would use an "about box". - - You should also get your employer (if you work as a programmer) or school, -if any, to sign a "copyright disclaimer" for the program, if necessary. -For more information on this, and how to apply and follow the GNU GPL, see -. - - The GNU General Public License does not permit incorporating your program -into proprietary programs. If your program is a subroutine library, you -may consider it more useful to permit linking proprietary applications with -the library. If this is what you want to do, use the GNU Lesser General -Public License instead of this License. But first, please read -. diff --git a/PhysIO/code/tapas_physio_cardiac_detect_outliers.m b/PhysIO/code/tapas_physio_cardiac_detect_outliers.m deleted file mode 100644 index 0b1e64a0..00000000 --- a/PhysIO/code/tapas_physio_cardiac_detect_outliers.m +++ /dev/null @@ -1,120 +0,0 @@ -function [outliersHigh, outliersLow, verbose] = ... - tapas_physio_cardiac_detect_outliers(tCardiac,percentile,... - deviationPercentUp,deviationPercentDown, verbose, ah) -% detects outliers (missed or erroneous pulses) in heartrate given a sequence of heartbeat pulses -% -% output = tapas_physio_cardiac_detect_outliers(input) -% -% IN -% tCardiac [nPulses,1] onset time of cardiac pulses -% percentile -% deviationPercentUp -% deviationPercentDown -% isVerbose if false, only warnings are output, no figures -% ah axes handle (optional)...specifies where plot is provided -% -% OUT -% outliersHigh -% outliersLow -% fh -% -% EXAMPLE -% tapas_physio_cardiac_detect_outliers -% -% See also tapas_physio_correct_cardiac_pulses_manually -% -% Author: Lars Kasper -% Created: 2013-04-25 -% Copyright (C) 2013 TNU, Institute for Biomedical Engineering, University of Zurich and ETH Zurich. -% -% This file is part of the TNU CheckPhysRETROICOR toolbox, which is released under the terms of the GNU General Public -% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL -% (either version 3 or, at your option, any later version). For further details, see the file -% COPYING or . -% -% $Id$ - -dt = diff(tCardiac); - -if nargin < 5 - isVerbose = 1; -else - isVerbose = verbose.level >=1; -end - -if isVerbose - if nargin < 6 - fh = tapas_physio_get_default_fig_params(); - verbose.fig_handles(end+1,1) = fh; - set(fh, 'Name','Diagnostics raw phys time series'); - else - fh = get(ah, 'Parent'); - figure(fh); - set(fh, 'CurrentAxes', ah); - end - - - - hp(1) = plot(tCardiac(2:end), dt, 'Color', [0 0.5 0]); - xlabel('t (seconds)'); - ylabel('lag between heartbeats (seconds)'); - title('Temporal lag between heartbeats'); - legend('Temporal lag between subsequent heartbeats'); - -else - fh = []; -end - -if ~isempty(percentile) && ~isempty(deviationPercentUp) && ~isempty(deviationPercentDown) - - nBins = length(dt)/10; - [dtSort,dtInd]=sort(dt); - percentile=percentile/100; - upperThresh=(1+deviationPercentUp/100)*dtSort(ceil(percentile*length(dtSort))); - lowerThresh=(1-deviationPercentDown/100)*dtSort(ceil((1-percentile)*length(dtSort))); - outliersHigh=dtInd(find(dtSort>upperThresh)); - outliersLow=dtInd(find(dtSort 5 - warningMessage = sprintf(['%4.1f %% of all %d heartbeat durations are below %3.1f s ' ... - 'or above %3.1f s \n - consider refining the pulse detection algorithm!' ... - 'Alternatively, do not model cardiac and interaction noise terms'], ... - percentageOutliers, nHeartBeatDurations, ... - lowerThresh, upperThresh); - verbose = tapas_physio_log(warningMessage, verbose, 1); - end - -end - -end diff --git a/PhysIO/code/tapas_physio_cell2char.m b/PhysIO/code/tapas_physio_cell2char.m deleted file mode 100644 index 9176a792..00000000 --- a/PhysIO/code/tapas_physio_cell2char.m +++ /dev/null @@ -1,38 +0,0 @@ -function physio = tapas_physio_cell2char(physio) -%convert all cellstrings (e.g. created via matlabbatch) back into regular -%strings for all relevant parameters of physio-structure -% -% physio = tapas_physio_cell2char(physio) -% -% IN -% -% OUT -% -% EXAMPLE -% tapas_physio_cell2char -% -% See also -% -% Author: Lars Kasper -% Created: 2014-04-28 -% Copyright (C) 2014 TNU, Institute for Biomedical Engineering, University of Zurich and ETH Zurich. -% -% This file is part of the TAPAS PhysIO Toolbox, which is released under the terms of the GNU General Public -% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL -% (either version 3 or, at your option, any later version). For further details, see the file -% COPYING or . -% -% $Id: teditRETRO.m 464 2014-04-27 11:58:09Z kasperla $ -physio.save_dir = char(physio.save_dir); -physio.log_files.cardiac = char(physio.log_files.cardiac); -physio.log_files.respiration = char(physio.log_files.respiration); -physio.log_files.scan_timing = char(physio.log_files.scan_timing); -physio.model.movement.file_realignment_parameters = char(... - physio.model.movement.file_realignment_parameters); -physio.model.other.input_multiple_regressors = char(... - physio.model.other.input_multiple_regressors); -physio.model.output_multiple_regressors = char(... - physio.model.output_multiple_regressors); -physio.model.output_physio = char(... - physio.model.output_physio); -physio.verbose.fig_output_file = char(physio.verbose.fig_output_file); \ No newline at end of file diff --git a/PhysIO/code/tapas_physio_cfg_matlabbatch.m b/PhysIO/code/tapas_physio_cfg_matlabbatch.m deleted file mode 100644 index 69e618aa..00000000 --- a/PhysIO/code/tapas_physio_cfg_matlabbatch.m +++ /dev/null @@ -1,1482 +0,0 @@ -function physio = tapas_physio_cfg_matlabbatch -% Lars Kasper, March 2013 -% -% Copyright (C) 2013, Institute for Biomedical Engineering, ETH/Uni Zurich. -% -% This file is part of the TAPAS PhysIO Toolbox, which is released under the terms of the GNU General Public -% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL -% (either version 3 or, at your option, any later version). For further details, see the file -% COPYING or . -% -% $Id$ - - -pathThis = fileparts(mfilename('fullpath')); % TODO: more elegant via SPM! -addpath(pathThis); - - -%-------------------------------------------------------------------------- -%% Save_dir (directory) - where all data is saved -%-------------------------------------------------------------------------- -save_dir = cfg_files; -save_dir.tag = 'save_dir'; -save_dir.name = 'save_dir'; -save_dir.val = {{''}}; -save_dir.help = {'Specify a directory where output of modelling and figures shall be saved.' - 'Default: current working directory'}; -save_dir.filter = 'dir'; -save_dir.ufilter = '.*'; -save_dir.num = [0 1]; - -%========================================================================== -%% Sub-structure log_files -%========================================================================== - -%-------------------------------------------------------------------------- -% vendor -%-------------------------------------------------------------------------- -vendor = cfg_menu; -vendor.tag = 'vendor'; -vendor.name = 'vendor'; -vendor.help = {' vendor Name depending on your MR Scanner system' - ' ''Philips''' - ' ''GE''' - ' ''Siemens''' - ' ''Siemens_Tics'' - new Siemens physiological' - ' logging with time stamps in tics' - ' (= steps of 2.5 ms since midnight) and' - ' extra acquisition (scan_timing) logfile with' - ' time stamps of all volumes and slices' - ' ''Biopac_Mat'' - exported mat files from Biopac system' - ' ' - ' or' - ' ''Custom''' - ' ' - ' ''Custom'' expects the logfiles (separate files for cardiac and respiratory)' - ' to be plain text, with one cardiac (or' - ' respiratory) sample per row;' - ' If heartbeat (R-wave peak) events are' - ' recorded as well, they have to be put' - ' as a 2nd column in the cardiac logfile' - ' by specifying a 1; 0 in all other rows' - ' e.g.:' - ' 0.2 0' - ' 0.4 1 <- cardiac pulse event' - ' 0.2 0' - ' -0.3 0' - ' ' - ' ' - ' NOTE: the sampling interval has to be specified for these files as' - ' well (s.b.)' - }; -vendor.labels = {'Philips', 'GE', 'Siemens (VB, *.puls/*.ecg/*.resp)', 'Siemens_Tics (VD: *_PULS.log/*_ECG1.log/*_RESP.log/*_AcquisitionInfo*.log)', 'Biopac_Mat', 'Custom'}; -vendor.values = {'Philips', 'GE', 'Siemens', 'Siemens_Tics', 'Biopac_Mat', 'Custom'}; -vendor.val = {'Philips'}; - -%-------------------------------------------------------------------------- -% cardiac -%-------------------------------------------------------------------------- -cardiac = cfg_files; -cardiac.tag = 'cardiac'; -cardiac.name = 'log_cardiac'; -%cardiac.val = {{'/Users/kasperla/Documents/code/matlab/smoothing_trunk/tSNR_fMRI_SPM/CheckPhysRETROICOR/PhysIOToolbox/examples/Philips/ECG3T/SCANPHYSLOG.log'}}; -cardiac.help = {'logfile with cardiac, i.e. ECG/PPU (pulse oximetry) data' - 'Select 0 files, if only respiratory data is available' - 'For Philips, same as respiratory logfile.' - }; -cardiac.filter = 'any'; -cardiac.ufilter = '.*'; -cardiac.num = [0 1]; - -%-------------------------------------------------------------------------- -% respiration (filename) -%-------------------------------------------------------------------------- -respiration = cfg_files; -respiration.tag = 'respiration'; -respiration.name = 'log_respiration'; -% respiration.val = {{'/Users/kasperla/Documents/code/matlab/smoothing_trunk/tSNR_fMRI_SPM/CheckPhysRETROICOR/PhysIOToolbox/examples/Philips/ECG3T/SCANPHYSLOG.log'}}; -respiration.help = {'logfile with respiratory, i.e. breathing belt amplitude data' - 'Select 0 files, if only cardiac data available' - 'For Philips, same as cardiac logfile.' - }; -respiration.filter = 'any'; -respiration.ufilter = '.*'; -respiration.num = [0 1]; - -%-------------------------------------------------------------------------- -% respiration (filename) -%-------------------------------------------------------------------------- -log_scan_timing = cfg_files; -log_scan_timing.tag = 'scan_timing'; -log_scan_timing.name = 'log_scan_timing'; -log_scan_timing.help = { - 'additional file for relative timing information between logfiles and' - ' MRI scans.' - '' - ' Currently implemented for 2 cases:' - ' Siemens: Enter the first or last Dicom volume of your session here,' - ' The time stamp in the dicom header is on the same time' - ' axis as the time stamp in the physiological log file' - ' Siemens_Tics: log-file which holds table conversion for tics axis to' - ' time conversion' - }; -log_scan_timing.filter = 'any'; -log_scan_timing.ufilter = '.*'; -log_scan_timing.val = {{''}}; -log_scan_timing.num = [0 1]; - -%-------------------------------------------------------------------------- -% sampling_interval -%-------------------------------------------------------------------------- -sampling_interval = cfg_entry; -sampling_interval.tag = 'sampling_interval'; -sampling_interval.name = 'sampling_interval'; -sampling_interval.help = { - 'sampling interval of phys log files (in seconds)' - ' If empty, default values are used: 2 ms for Philips, 25 ms for GE and others' - ' For Biopac, sampling rate is read directly from logfile' - ' If cardiac and respiratory sampling rate differ, enter them as vector' - ' [sampling_interval_cardiac, sampling_interval_respiratory]' - ' ' - ' If cardiac, respiratory and acquisition timing (tics) sampling rate differ,' - 'enter them as a vector:' - ' [sampling_interval_cardiac, sampling_interval_respiratory sampling_interval_tics_acquisition_timing]' - '' - ' Note: If you use a WiFi-Philips device for peripheral monitoring' - ' (Ingenia system), please change this value to 1/496, ' - ' i.e. a sampling rate of 496 Hz)' - }; -sampling_interval.strtype = 'e'; -sampling_interval.num = [Inf Inf]; -sampling_interval.val = {[]}; - -%-------------------------------------------------------------------------- -% relative_start_acquisition -%-------------------------------------------------------------------------- -relative_start_acquisition = cfg_entry; -relative_start_acquisition.tag = 'relative_start_acquisition'; -relative_start_acquisition.name = 'relative_start_acquisition'; -relative_start_acquisition.help = { - 'start time (ins seconds) of 1st scan (or dummy)' - 'relative to start of physiological logfile'}; -relative_start_acquisition.strtype = 'e'; -relative_start_acquisition.num = [Inf Inf]; -relative_start_acquisition.val = {0}; - - -%-------------------------------------------------------------------------- -% align_scan -%-------------------------------------------------------------------------- -align_scan = cfg_menu; -align_scan.tag = 'align_scan'; -align_scan.name = 'align_scan'; -align_scan.help = { - ' Determines which scan shall be aligned to which part of the logfile' - ' Typically, aligning the last scan to the end of the logfile is' - ' beneficial, since start of logfile and scans might be shifted due' - ' to pre-scans' - '' - ' NOTE: In all cases, log_files.relative_start_acquisition is' - ' added to timing after the initial alignmnent to first/last scan' - '' - ' ''first'' start of logfile will be aligned to first scan volume' - ' ''last'' end of logfile will be aligned to last scan volume' - - }; -align_scan.labels = {'first', 'last'}; -align_scan.values = {'first', 'last'}; -align_scan.val = {'last'}; - - -%-------------------------------------------------------------------------- -% files -%-------------------------------------------------------------------------- -files = cfg_branch; -files.tag = 'log_files'; -files.name = 'log_files'; -files.val = {vendor cardiac respiration log_scan_timing, ... - sampling_interval, relative_start_acquisition, align_scan}; -files.help = {'Specify log files where peripheral data was stored, and their properties.'}; - - - - -%========================================================================== -%% Subsub-structure sqpar -%========================================================================== - - - -%-------------------------------------------------------------------------- -% Nscans -%-------------------------------------------------------------------------- -Nscans = cfg_entry; -Nscans.tag = 'Nscans'; -Nscans.name = 'Nscans'; -Nscans.help = { - 'Number of scans (volumes) in design matrix.' - 'Put exactly the same number as you have image volumes in your SPM GLM' - 'design specification.' - }; -Nscans.strtype = 'e'; -Nscans.num = [Inf Inf]; -%Nscans.val = {495}; - -%-------------------------------------------------------------------------- -% Ndummies -%-------------------------------------------------------------------------- -Ndummies = cfg_entry; -Ndummies.tag = 'Ndummies'; -Ndummies.name = 'Ndummies'; -Ndummies.help = { - 'Number of dummies that were acquired (but will not show up in design matrix' - '(also enter correct number, if dummies are not saved in imaging file)' - }; -Ndummies.strtype = 'e'; -Ndummies.num = [Inf Inf]; -%Ndummies.val = {3}; - -%-------------------------------------------------------------------------- -% TR -%-------------------------------------------------------------------------- -TR = cfg_entry; -TR.tag = 'TR'; -TR.name = 'TR'; -TR.help = {'Repetition time (in seconds) between consecutive image volumes'}; -TR.strtype = 'e'; -TR.num = [Inf Inf]; -%TR.val = {2.5}; - -%-------------------------------------------------------------------------- -% NslicesPerBeat -%-------------------------------------------------------------------------- -NslicesPerBeat = cfg_entry; -NslicesPerBeat.tag = 'NslicesPerBeat'; -NslicesPerBeat.name = 'NslicesPerBeat'; -NslicesPerBeat.help = {'Only for triggered (gated) sequences: ' - 'Number of slices acquired per heartbeat'}; -NslicesPerBeat.strtype = 'e'; -NslicesPerBeat.num = [Inf Inf]; -NslicesPerBeat.val = {[]}; - - -%-------------------------------------------------------------------------- -% Nslices -%-------------------------------------------------------------------------- -Nslices = cfg_entry; -Nslices.tag = 'Nslices'; -Nslices.name = 'Nslices'; -Nslices.help = {'Number of slices in one volume'}; -Nslices.strtype = 'e'; -Nslices.num = [Inf Inf]; -%Nslices.val = {37}; - - - -%-------------------------------------------------------------------------- -% onset_slice -%-------------------------------------------------------------------------- -onset_slice = cfg_entry; -onset_slice.tag = 'onset_slice'; -onset_slice.name = 'onset_slice'; -onset_slice.help = { - 'Slice to which regressors are temporally aligned.' - 'Typically the slice where your most important activation is expected.'}; -onset_slice.strtype = 'e'; -onset_slice.num = [Inf Inf]; -%onset_slice.val = {19}; - -%-------------------------------------------------------------------------- -% Nprep -%-------------------------------------------------------------------------- -Nprep = cfg_entry; -Nprep.tag = 'Nprep'; -Nprep.name = 'Nprep'; -Nprep.help = { - ' Count of preparation pulses BEFORE 1st dummy scan.' - ' Only important, if log_files.scan_align = ''first'', since then' - ' preparation pulses and dummiy triggers are counted and discarded ' - ' as first scan onset' - }; -Nprep.strtype = 'e'; -Nprep.num = [Inf Inf]; -Nprep.val = {[]}; - -%-------------------------------------------------------------------------- -% time_slice_to_slice -%-------------------------------------------------------------------------- -time_slice_to_slice = cfg_entry; -time_slice_to_slice.tag = 'time_slice_to_slice'; -time_slice_to_slice.name = 'time_slice_to_slice'; -time_slice_to_slice.help = { - 'Duration between acquisition of two different slices' - 'if empty, set to default value TR/Nslices' - 'differs e.g. if slice timing was minimal and TR was bigger than needed' - 'to acquire Nslices' - }; -time_slice_to_slice.strtype = 'e'; -time_slice_to_slice.num = [Inf Inf]; -time_slice_to_slice.val = {[]}; - -%-------------------------------------------------------------------------- -% sqpar -%-------------------------------------------------------------------------- -sqpar = cfg_branch; -sqpar.tag = 'sqpar'; -sqpar.name = 'sqpar (Sequence timing parameters)'; -sqpar.val = {Nslices NslicesPerBeat TR Ndummies Nscans onset_slice time_slice_to_slice Nprep}; -sqpar.help = {'Sequence timing parameters, (number of slices, volumes, dummies, volume TR, slice TR ...)'}; - - -% ========================================================================== -%% Subsub-structure sync -%========================================================================== - - -%-------------------------------------------------------------------------- -% grad_direction -%-------------------------------------------------------------------------- -grad_direction = cfg_menu; -grad_direction.tag = 'grad_direction'; -grad_direction.name = 'grad_direction'; -grad_direction.help = {'...'}; -grad_direction.labels = {'x' 'y' 'z'}; -grad_direction.values = {'x' 'y' 'z'}; -grad_direction.val = {'y'}; - -%-------------------------------------------------------------------------- -% vol_spacing -%-------------------------------------------------------------------------- -vol_spacing = cfg_entry; -vol_spacing.tag = 'vol_spacing'; -vol_spacing.name = 'vol_spacing'; -vol_spacing.help = {'time (in seconds) between last slice of n-th volume' - 'and 1st slice of n+1-th volume(overrides .vol-threshold)' - 'NOTE: Leave empty if .vol shall be used'}; -vol_spacing.strtype = 'e'; -vol_spacing.num = [Inf Inf]; -vol_spacing.val = {[]}; - -%-------------------------------------------------------------------------- -% vol -%-------------------------------------------------------------------------- -vol = cfg_entry; -vol.tag = 'vol'; -vol.name = 'vol'; -vol.help = {'Gradient Amplitude Threshold for Start of new Volume'}; -vol.strtype = 'e'; -vol.num = [Inf Inf]; -vol.val = {[]}; - -%-------------------------------------------------------------------------- -% slice -%-------------------------------------------------------------------------- -slice = cfg_entry; -slice.tag = 'slice'; -slice.name = 'slice'; -slice.help = {'Gradient Amplitude Threshold for Start of new slice'}; -slice.strtype = 'e'; -slice.num = [Inf Inf]; -slice.val = {1800}; - -%-------------------------------------------------------------------------- -% zero -%-------------------------------------------------------------------------- -zero = cfg_entry; -zero.tag = 'zero'; -zero.name = 'zero'; -zero.help = {'Gradient Amplitude Threshold below which values will be set to 0.'}; -zero.strtype = 'e'; -zero.num = [Inf Inf]; -zero.val = {1700}; - - -%-------------------------------------------------------------------------- -% sync_method_gradient_log -%-------------------------------------------------------------------------- - -sync_method_gradient_log = cfg_branch; -sync_method_gradient_log.tag = 'gradient_log'; -sync_method_gradient_log.name = 'gradient_log'; -sync_method_gradient_log.val = { - grad_direction - zero - slice - vol - vol_spacing -}; -sync_method_gradient_log.help = { ... - ' Derive scan-timing from logged gradient time courses' - ' in SCANPHYSLOG-files (Philips only)'}; - - -%-------------------------------------------------------------------------- -% sync_method_gradient_log_auto -%-------------------------------------------------------------------------- - -sync_method_gradient_log_auto = cfg_branch; -sync_method_gradient_log_auto.tag = 'gradient_log_auto'; -sync_method_gradient_log_auto.name = 'gradient_log_auto'; -sync_method_gradient_log_auto.val = {}; -sync_method_gradient_log_auto.help = { ... - ' Derive scan-timing from logged gradient time courses' - ' in SCANPHYSLOG-files automatically (Philips only), ' - ' using prior information on TR and number of slices, ' - 'i.e. without manual threshold settings.' -}; - - -%-------------------------------------------------------------------------- -% sync_method_nominal -%-------------------------------------------------------------------------- - -sync_method_nominal = cfg_branch; -sync_method_nominal.tag = 'nominal'; -sync_method_nominal.name = 'nominal'; -sync_method_nominal.val = {}; -sync_method_nominal.help = { ... - ' Derive scan-timing for sqpar (nominal scan timing parameters)'}; - - -%-------------------------------------------------------------------------- -% sync_method_sync_log -%-------------------------------------------------------------------------- - -sync_method_scan_timing_log = cfg_branch; -sync_method_scan_timing_log.tag = 'scan_timing_log'; -sync_method_scan_timing_log.name = 'scan_timing_log'; -sync_method_scan_timing_log.val = {}; -sync_method_scan_timing_log.help = { ... - ' Derive scan-timing from individual scan timing logfile with time ' - ' stamps ("tics") for each slice and volume (e.g. Siemens_Cologne)'}; - - -%-------------------------------------------------------------------------- -% sync -%-------------------------------------------------------------------------- - -sync = cfg_choice; -sync.tag = 'sync'; -sync.name = 'Scan/Physlog Time Synchronization'; -sync.values = {sync_method_nominal, ... - sync_method_gradient_log, ... - sync_method_gradient_log_auto, ... - sync_method_scan_timing_log}; -sync.val = {sync_method_nominal}; -sync.help = {'Determines scan timing from nominal scan parameters or logged gradient time courses' - '' -' Available methods to determine slice onset times' -' ''nominal'' - to derive slice acquisition timing from sqpar directly' -' ''gradient_log'' - derive from logged gradient time courses' -' in SCANPHYSLOG-files (Philips only)' -' ''gradient_log_auto'' - as gradient_log, but thresholds are determined' -' automatically from TR and number of slices expected' -' ''scan_timing_log'' - individual scan timing logfile with time stamps ("tics") for each slice and volume (e.g. Siemens_Cologne)' -}; - - -%-------------------------------------------------------------------------- -% scan_timing -%-------------------------------------------------------------------------- -scan_timing = cfg_branch; -scan_timing.tag = 'scan_timing'; -scan_timing.name = 'scan_timing (Parameters for sequence timing & synchronization)'; -scan_timing.val = {sqpar sync}; -scan_timing.help = {'Parameters for sequence timing & synchronization, i.e.' - 'scan_tming.sqpar = slice and volume acquisition starts, TR,' - ' number of scans etc.' - 'scan_timing.sync = synchronize phys logfile to scan acquisition via logged MR gradient time courses/time stamps' - }; - - -% ========================================================================== -%% Sub-structure preproc -%========================================================================== - -% ========================================================================== -%% Subsub-structure cardiac -%========================================================================== - - -%-------------------------------------------------------------------------- -% modality -%-------------------------------------------------------------------------- -modality = cfg_menu; -modality.tag = 'modality'; -modality.name = 'modality'; -modality.help = {'Shall ECG or PPU data be read from logfiles?'}; -modality.labels = {'ECG', 'OXY/PPU', 'ECG_WiFi', 'PPU_WiFi'}; -modality.values = {'ECG', 'PPU', 'ECG_WiFi', 'PPU_Wifi'}; -modality.val = {'ECG'}; - - - -%-------------------------------------------------------------------------- -% initial_cpulse_select_file -%-------------------------------------------------------------------------- -initial_cpulse_select_file = cfg_entry; -initial_cpulse_select_file.tag = 'file'; -initial_cpulse_select_file.name = 'file'; -initial_cpulse_select_file.help = {'...'}; -initial_cpulse_select_file.strtype = 's'; -initial_cpulse_select_file.num = [0 Inf]; -initial_cpulse_select_file.val = {'initial_cpulse_kRpeakfile.mat'}; - - - -%-------------------------------------------------------------------------- -% min -%-------------------------------------------------------------------------- -min = cfg_entry; -min.tag = 'min'; -min.name = 'min'; -min.help = {'Minimum correlation value considered a peak (for auto, manual, load-methods).'}; -min.strtype = 'e'; -min.num = [Inf Inf]; -min.val = {0.4}; - - -%-------------------------------------------------------------------------- -% initial_cpulse_select_method_auto_template -%-------------------------------------------------------------------------- - -initial_cpulse_select_method_auto_template = cfg_branch; -initial_cpulse_select_method_auto_template.tag = 'auto_template'; -initial_cpulse_select_method_auto_template.name = 'auto_template'; -initial_cpulse_select_method_auto_template.val = { - min - initial_cpulse_select_file -}; -initial_cpulse_select_method_auto_template.help = { ... - ' Auto generation of representative QRS-wave; detection via' - ' maximising auto-correlation with it (default)'}; - - -%-------------------------------------------------------------------------- -% initial_cpulse_select_method_auto_matched -%-------------------------------------------------------------------------- - -initial_cpulse_select_method_auto_matched = cfg_branch; -initial_cpulse_select_method_auto_matched.tag = 'auto_matched'; -initial_cpulse_select_method_auto_matched.name = 'auto_matched'; -initial_cpulse_select_method_auto_matched.val = { - min - initial_cpulse_select_file -}; -initial_cpulse_select_method_auto_matched.help = { ... - 'Auto generation of template QRS wave, ' - ' matched-filter/autocorrelation detection of heartbeats' - }; - - -%-------------------------------------------------------------------------- -% initial_cpulse_select_method_manual_template -%-------------------------------------------------------------------------- - -initial_cpulse_select_method_manual_template = cfg_branch; -initial_cpulse_select_method_manual_template.tag = 'manual_template'; -initial_cpulse_select_method_manual_template.name = 'manual_template'; -initial_cpulse_select_method_manual_template.val = { - min - initial_cpulse_select_file -}; -initial_cpulse_select_method_manual_template.help = { ... - 'Manually select QRS-wave for autocorrelation detection' - }; - - -%-------------------------------------------------------------------------- -% initial_cpulse_select_method_load_template -%-------------------------------------------------------------------------- - -initial_cpulse_select_method_load_template = cfg_branch; -initial_cpulse_select_method_load_template.tag = 'load_template'; -initial_cpulse_select_method_load_template.name = 'load_template'; -initial_cpulse_select_method_load_template.val = { - min - initial_cpulse_select_file -}; -initial_cpulse_select_method_load_template.help = { ... - 'Load template from previous manual/auto run to perform autocorrelation detection of hearbeats' - }; - - - -%-------------------------------------------------------------------------- -% initial_cpulse_select_method_load_from_logfile -%-------------------------------------------------------------------------- - -initial_cpulse_select_method_load_from_logfile = cfg_branch; -initial_cpulse_select_method_load_from_logfile.tag = 'load_from_logfile'; -initial_cpulse_select_method_load_from_logfile.name = 'load_from_logfile'; -initial_cpulse_select_method_load_from_logfile.val = {}; -initial_cpulse_select_method_load_from_logfile.help = { ... - 'Load heartbeat data from Phys-logfile, detected R-peaks of scanner'}; - - -%-------------------------------------------------------------------------- -% initial_cpulse_select -%-------------------------------------------------------------------------- -initial_cpulse_select = cfg_choice; -initial_cpulse_select.tag = 'initial_cpulse_select'; -initial_cpulse_select.name = 'Initial Detection of Heartbeats'; -initial_cpulse_select.val = {initial_cpulse_select_method_load_from_logfile}; -initial_cpulse_select.values = { - initial_cpulse_select_method_auto_matched, ... - initial_cpulse_select_method_auto_template, ... - initial_cpulse_select_method_load_from_logfile, ... - initial_cpulse_select_method_manual_template, ... - initial_cpulse_select_method_load_template, ... - }; - - -initial_cpulse_select.help = { - 'The initial cardiac pulse selection structure: Determines how the' - 'majority of cardiac pulses is detected in a first pass.' - ' ''auto_matched'' - auto generation of template QRS wave, ' - ' matched-filter/autocorrelation detection of heartbeats' - ' ''auto_template'' - auto generation of representative QRS-wave; detection via' - ' maximising auto-correlation with it (default)' - ' ''load_from_logfile'' - from phys logfile, detected R-peaks of scanner' - ' ''manual_template'' - via manually selected QRS-wave for autocorrelations' - ' ''load_template'' - from previous manual/auto run' - }; - - -%-------------------------------------------------------------------------- -% posthoc_cpulse_select_file -%-------------------------------------------------------------------------- -posthoc_cpulse_select_file = cfg_entry; -posthoc_cpulse_select_file.tag = 'file'; -posthoc_cpulse_select_file.name = 'file'; -posthoc_cpulse_select_file.help = {'...'}; -posthoc_cpulse_select_file.strtype = 's'; -posthoc_cpulse_select_file.num = [0 Inf]; -posthoc_cpulse_select_file.val = {'posthoc_cpulse.mat'}; - - -%-------------------------------------------------------------------------- -% posthoc_cpulse_select_percentile -%-------------------------------------------------------------------------- - -posthoc_cpulse_select_percentile = cfg_entry; -posthoc_cpulse_select_percentile.tag = 'percentile'; -posthoc_cpulse_select_percentile.name = 'percentile'; -posthoc_cpulse_select_percentile.help = { - 'percentile of beat-2-beat interval histogram that constitutes the' - 'average heart beat duration in the session'}; -posthoc_cpulse_select_percentile.strtype = 'e'; -posthoc_cpulse_select_percentile.num = [Inf Inf]; -posthoc_cpulse_select_percentile.val = {80}; - - -%-------------------------------------------------------------------------- -% posthoc_cpulse_select_upper_thresh -%-------------------------------------------------------------------------- -posthoc_cpulse_select_upper_thresh = cfg_entry; -posthoc_cpulse_select_upper_thresh.tag = 'upper_thresh'; -posthoc_cpulse_select_upper_thresh.name = 'upper_thresh'; -posthoc_cpulse_select_upper_thresh.help = { - 'minimum exceedance (in %) from average heartbeat duration ' - 'to be classified as missing heartbeat'}; -posthoc_cpulse_select_upper_thresh.strtype = 'e'; -posthoc_cpulse_select_upper_thresh.num = [Inf Inf]; -posthoc_cpulse_select_upper_thresh.val = {60}; - - -%-------------------------------------------------------------------------- -% posthoc_cpulse_select_lower_thresh -%-------------------------------------------------------------------------- - -posthoc_cpulse_select_lower_thresh = cfg_entry; -posthoc_cpulse_select_lower_thresh.tag = 'lower_thresh'; -posthoc_cpulse_select_lower_thresh.name = 'lower_thresh'; -posthoc_cpulse_select_lower_thresh.help = { - 'minimum reduction (in %) from average heartbeat duration' - 'to be classified an abundant heartbeat'}; -posthoc_cpulse_select_lower_thresh.strtype = 'e'; -posthoc_cpulse_select_lower_thresh.num = [Inf Inf]; -posthoc_cpulse_select_lower_thresh.val = {60}; - - -%-------------------------------------------------------------------------- -% posthoc_cpulse_select_method_off -%-------------------------------------------------------------------------- - -posthoc_cpulse_select_method_off = cfg_branch; -posthoc_cpulse_select_method_off.tag = 'off'; -posthoc_cpulse_select_method_off.name = 'Off'; -posthoc_cpulse_select_method_off.val = {}; -posthoc_cpulse_select_method_off.help = {'No manual post-hoc pulse selection'}; - - -%-------------------------------------------------------------------------- -% posthoc_cpulse_select_method_manual -%-------------------------------------------------------------------------- - -posthoc_cpulse_select_method_manual = cfg_branch; -posthoc_cpulse_select_method_manual.tag = 'manual'; -posthoc_cpulse_select_method_manual.name = 'Manual'; -posthoc_cpulse_select_method_manual.help = {'Manual post-hoc cardiac pulse selection by clicking'}; - -posthoc_cpulse_select_method_manual.val = {... - posthoc_cpulse_select_file ... - posthoc_cpulse_select_percentile ... - posthoc_cpulse_select_upper_thresh ... - posthoc_cpulse_select_lower_thresh}; - - -%-------------------------------------------------------------------------- -% posthoc_cpulse_select_method_load -%-------------------------------------------------------------------------- - -posthoc_cpulse_select_method_load = cfg_branch; -posthoc_cpulse_select_method_load.tag = 'load'; -posthoc_cpulse_select_method_load.name = 'Load'; -posthoc_cpulse_select_method_load.help = {'Loads manually selected cardiac pulses from file'}; -posthoc_cpulse_select_method_load.val = { - posthoc_cpulse_select_file}; - - -%-------------------------------------------------------------------------- -% posthoc_cpulse_select -%-------------------------------------------------------------------------- - -posthoc_cpulse_select = cfg_choice; -posthoc_cpulse_select.tag = 'posthoc_cpulse_select'; -posthoc_cpulse_select.name = 'Post-hoc Selection of Cardiac Pulses'; -posthoc_cpulse_select.val = {posthoc_cpulse_select_method_off}; -posthoc_cpulse_select.values = {posthoc_cpulse_select_method_off, ... - posthoc_cpulse_select_method_manual, ... - posthoc_cpulse_select_method_load}; - - -posthoc_cpulse_select.help = { - 'The posthoc cardiac pulse selection structure: If only few (<20)' - 'cardiac pulses are missing in a session due to bad signal quality, a' - 'manual selection after visual inspection is possible using the' - 'following parameters. The results are saved for reproducibility.' - '' - 'Refers to physio.preproc.cardiac.posthoc_cpulse_select.method in physio-structure' - }; - - - -%-------------------------------------------------------------------------- -% cardiac -%-------------------------------------------------------------------------- -cardiac = cfg_branch; -cardiac.tag = 'cardiac'; -cardiac.name = 'cardiac'; -cardiac.val = {modality initial_cpulse_select posthoc_cpulse_select}; -cardiac.help = {'...'}; - - -%-------------------------------------------------------------------------- -% preproc -%-------------------------------------------------------------------------- -preproc = cfg_branch; -preproc.tag = 'preproc'; -preproc.name = 'preproc (Thresholding parameters for de-noising and timing)'; -preproc.val = {cardiac}; -preproc.help = {'Thresholding parameters for de-noising of raw peripheral data' - 'and determination of sequence timing from logged MR gradient time courses'}; - - - -%========================================================================== -%% Sub-structure model -%========================================================================== - - -%-------------------------------------------------------------------------- -% orthog -%-------------------------------------------------------------------------- - -orthog = cfg_menu; -orthog.tag = 'orthogonalise'; -orthog.name = 'orthogonalise'; -orthog.help = { - 'Orthogonalize physiological regressors with respect to each other.' - 'Note: This is only recommended for triggered/gated acquisition sequences.' - }; -orthog.labels = {'none' 'cardiac' 'resp' 'mult' 'RETROCOR', 'HRV', 'RVT', 'Noise_ROIs'}; -orthog.values = {'none' 'cardiac' 'resp' 'mult' 'retroicor', 'hrv', 'rvt', 'noise_rois'}; -orthog.val = {'none'}; - - -%-------------------------------------------------------------------------- -% output_multiple_regressors -%-------------------------------------------------------------------------- - -output_multiple_regressors = cfg_entry; -output_multiple_regressors.tag = 'output_multiple_regressors'; -output_multiple_regressors.name = 'output_multiple_regressors'; -output_multiple_regressors.help = { - 'Output file for physiological regressors' - 'Choose file name with extension:' - '.txt for ASCII files with 1 regressor per column' - '.mat for matlab variable file' - }; -output_multiple_regressors.strtype = 's'; -output_multiple_regressors.num = [1 Inf]; -output_multiple_regressors.val = {'multiple_regressors.txt'}; - - -%-------------------------------------------------------------------------- -% output_physio -%-------------------------------------------------------------------------- - -output_physio = cfg_entry; -output_physio.tag = 'output_physio'; -output_physio.name = 'output_physio'; -output_physio.help = { - 'Output file for physio-structure with extracted physiological time' - 'series, detected peak and created regressors' - 'Choose mat-file name; structure will be saved as variable physio in there.' - }; -output_physio.strtype = 's'; -output_physio.num = [1 Inf]; -output_physio.val = {'physio.mat'}; - - -%-------------------------------------------------------------------------- -%% Sub-substructure retroicor -%-------------------------------------------------------------------------- - -%-------------------------------------------------------------------------- -% c -%-------------------------------------------------------------------------- -c = cfg_entry; -c.tag = 'c'; -c.name = 'cardiac'; -c.help = {'Order of Fourier expansion for cardiac phase' - ' - equals 1/2 number of cardiac regressors, since sine and cosine terms' - 'are computed, i.e. sin(phi), cos(phi), sin(2*phi), cos(2*phi), ..., sin(c*phi), cos(c*phi)' - }; -c.strtype = 'e'; -c.num = [1 1]; -c.val = {3}; - -%-------------------------------------------------------------------------- -% r -%-------------------------------------------------------------------------- -r = cfg_entry; -r.tag = 'r'; -r.name = 'respiratory'; -r.help = { - 'Order of Fourier expansion for respiratory phase' - ' - equals 1/2 number of respiratory regressors, since sine and cosine terms' - 'are computed, i.e. sin(phi), cos(phi), sin(2*phi), cos(2*phi), ..., sin(r*phi), cos(r*phi)' - }; -r.strtype = 'e'; -r.num = [1 1]; -r.val = {4}; - -%-------------------------------------------------------------------------- -% cr -%-------------------------------------------------------------------------- -cr = cfg_entry; -cr.tag = 'cr'; -cr.name = 'cardiac X respiratory'; -cr.help = { - 'Order of Fourier expansion for interaction of cardiac and respiratory phase' - ' - equals 1/4 number of interaction regressors, since sine and cosine terms' - 'are computed and multiplied, i.e. sin(phi_c)*cos(phi_r), sin(phi_r)*cos(phi_c)' - }; -cr.strtype = 'e'; -cr.num = [1 1]; -cr.val = {1}; - - -%-------------------------------------------------------------------------- -% order -%-------------------------------------------------------------------------- - -order = cfg_branch; -order.tag = 'order'; -order.name = 'order'; -order.val = {c r cr}; -order.help = {'...'}; - - -%-------------------------------------------------------------------------- -% retroicor_no -%-------------------------------------------------------------------------- - -retroicor_no = cfg_branch; -retroicor_no.tag = 'no'; -retroicor_no.name = 'No'; -retroicor_no.val = {}; -retroicor_no.help = {'RETROICOR not used'}; - - -%-------------------------------------------------------------------------- -% retroicor_yes -%-------------------------------------------------------------------------- - -retroicor_yes = cfg_branch; -retroicor_yes.tag = 'yes'; -retroicor_yes.name = 'Yes'; -retroicor_yes.val = {order}; -retroicor_yes.help = {'Include RETROICOR Model, as described in Glover et al., MRM 2000'}; - - - -%-------------------------------------------------------------------------- -% retroicor -%-------------------------------------------------------------------------- - -retroicor = cfg_choice; -retroicor.tag = 'retroicor'; -retroicor.name = 'RETROICOR'; -retroicor.val = {retroicor_yes}; -retroicor.values = {retroicor_no, retroicor_yes}; -retroicor.help = {'RETROICOR Model, as described in Glover et al., MRM 2000'}; - - -%-------------------------------------------------------------------------- -%% Sub-substructure RVT -%-------------------------------------------------------------------------- - -%-------------------------------------------------------------------------- -% rvt_delays -%-------------------------------------------------------------------------- -rvt_delays = cfg_entry; -rvt_delays.tag = 'delays'; -rvt_delays.name = 'Delays (seconds)'; -rvt_delays.help = { - 'Delays (in seconds) by which respiratory response function is ' - 'shifted with respect to RVT regressor before convolution' - }; -rvt_delays.num = [Inf Inf]; -rvt_delays.val = {0}; - -%-------------------------------------------------------------------------- -% rvt_no -%-------------------------------------------------------------------------- - -rvt_no = cfg_branch; -rvt_no.tag = 'no'; -rvt_no.name = 'No'; -rvt_no.val = {}; -rvt_no.help = {'Respiratory Volume per Time Model not used'}; - - -%-------------------------------------------------------------------------- -% rvt_yes -%-------------------------------------------------------------------------- - -rvt_yes = cfg_branch; -rvt_yes.tag = 'yes'; -rvt_yes.name = 'Yes'; -rvt_yes.val = {rvt_delays}; -rvt_yes.help = { - 'Include Respiratory Volume per Time (RVT) Model, ' - 'as described in Birn, R.M., et al. NeuroImage 40, 644?654. doi:10.1016/j.neuroimage.2007.11.059' - }; - - - -%-------------------------------------------------------------------------- -% rvt -%-------------------------------------------------------------------------- - -rvt = cfg_choice; -rvt.tag = 'rvt'; -rvt.name = 'Respiratory Volume per Time (RVT)'; -rvt.val = {rvt_no}; -rvt.values = {rvt_no, rvt_yes}; -rvt.help = { - 'Respiratory Volume per Time (RVT) Model, ' - 'as described in Birn, R.M., et al. NeuroImage 40, 644-654. doi:10.1016/j.neuroimage.2007.11.059' - }; - - -%-------------------------------------------------------------------------- -%% Sub-substructure HRV -%-------------------------------------------------------------------------- - -%-------------------------------------------------------------------------- -% hrv_delays -%-------------------------------------------------------------------------- -hrv_delays = cfg_entry; -hrv_delays.tag = 'delays'; -hrv_delays.name = 'Delays (seconds)'; -hrv_delays.help = { - 'Delays (in seconds) by which respiratory response function is ' - 'shifted with respect to HRV regressor before convolution' - }; -hrv_delays.num = [Inf Inf]; -hrv_delays.val = {0}; - -%-------------------------------------------------------------------------- -% hrv_no -%-------------------------------------------------------------------------- - -hrv_no = cfg_branch; -hrv_no.tag = 'no'; -hrv_no.name = 'No'; -hrv_no.val = {}; -hrv_no.help = {'Heart Rate Variability Model not used'}; - - -%-------------------------------------------------------------------------- -% hrv_yes -%-------------------------------------------------------------------------- - -hrv_yes = cfg_branch; -hrv_yes.tag = 'yes'; -hrv_yes.name = 'Yes'; -hrv_yes.val = {hrv_delays}; -hrv_yes.help = { - 'Include Heart Rate Variability (HRV) Model, ' - 'as described in Chang, C. et al., NeuroImage 44, 857-869. doi:10.1016/j.neuroimage.2008.09.029' - }; - - - -%-------------------------------------------------------------------------- -% hrv -%-------------------------------------------------------------------------- - -hrv = cfg_choice; -hrv.tag = 'hrv'; -hrv.name = 'Heart Rate Variability (HRV)'; -hrv.val = {hrv_no}; -hrv.values = {hrv_no, hrv_yes}; -hrv.help = { - 'Heart Rate Variability (HRV) Model, as described in ' - 'Chang, C. et al., NeuroImage 44, 857-869. doi:10.1016/j.neuroimage.2008.09.029' -}; - - -%-------------------------------------------------------------------------- -%% Sub-substructure Noise_Rois Model -%-------------------------------------------------------------------------- - -%-------------------------------------------------------------------------- -% fmri_files -%-------------------------------------------------------------------------- - -fmri_files = cfg_files; -fmri_files.tag = 'fmri_files'; -fmri_files.name = 'FMRI Time Series File(s)'; -fmri_files.val = {{''}}; -fmri_files.help = { - 'Preprocessed fmri nifti/analyze files, from which time series ' - 'shall be extracted'}; -fmri_files.filter = '.*'; -fmri_files.ufilter = '.nii$|.img$'; -fmri_files.num = [0 Inf]; - -%-------------------------------------------------------------------------- -% roi_files -%-------------------------------------------------------------------------- - -roi_files = cfg_files; -roi_files.tag = 'roi_files'; -roi_files.name = 'Noise ROI Image File(s)'; -roi_files.val = {{''}}; -roi_files.help = {'Masks/tissue probability maps characterizing where noise resides'}; -roi_files.filter = '.*'; -roi_files.ufilter = '.nii$|.img$'; -roi_files.num = [0 Inf]; - - -%-------------------------------------------------------------------------- -% roi_thresholds -%-------------------------------------------------------------------------- - -roi_thresholds = cfg_entry; -roi_thresholds.tag = 'thresholds'; -roi_thresholds.name = 'ROI thresholds'; -roi_thresholds.help = { - 'Single threshold or vector [1, nRois] of thresholds to be applied to' - 'mask files to decide which voxels to include ' - '(e.g. a probability like 0.99, if roi_files' - 'are tissue probability maps' - }; -roi_thresholds.num = [Inf Inf]; -roi_thresholds.val = {0.9}; - - -%-------------------------------------------------------------------------- -% n_voxel_crop -%-------------------------------------------------------------------------- - -n_voxel_crop = cfg_entry; -n_voxel_crop.tag = 'n_voxel_crop'; -n_voxel_crop.name = 'Number of ROI pixels to be cropped'; -n_voxel_crop.help = { - 'Single number or vector [1, nRois] of number of voxels to crop per ROI' - }; -n_voxel_crop.num = [Inf Inf]; -n_voxel_crop.val = {0}; - - -%-------------------------------------------------------------------------- -% n_components -%-------------------------------------------------------------------------- - -n_components = cfg_entry; -n_components.tag = 'n_components'; -n_components.name = 'Number of principal components'; -n_components.help = { - ' Single number or vector [1, nRois] of numbers' - ' integer >=1: number of principal components to be extracted' - ' from all voxel time series within each ROI' - ' float in [0,1[ choose as many components as needed to explain this' - ' relative share of total variance, e.g. 0.99 =' - ' add more components, until 99 % of variance explained' - ' NOTE: Additionally, the mean time series of the region is also' - ' extracted' }; -n_components.num = [Inf Inf]; -n_components.val = {1}; - - -%-------------------------------------------------------------------------- -% noise_rois_no -%-------------------------------------------------------------------------- - -noise_rois_no = cfg_branch; -noise_rois_no.tag = 'no'; -noise_rois_no.name = 'No'; -noise_rois_no.val = {}; -noise_rois_no.help = {'Noise ROIs not used'}; - - -%-------------------------------------------------------------------------- -% noise_rois_yes -%-------------------------------------------------------------------------- - -noise_rois_yes = cfg_branch; -noise_rois_yes.tag = 'yes'; -noise_rois_yes.name = 'Yes'; -noise_rois_yes.val = {fmri_files, roi_files, roi_thresholds, n_voxel_crop, ... - n_components}; -noise_rois_yes.help = { - 'Include Noise ROIs model' - '(Principal components of anatomical regions), similar to aCompCor, Behzadi et al. 2007' - }; - - - -%-------------------------------------------------------------------------- -% noise_rois -%-------------------------------------------------------------------------- - -noise_rois = cfg_choice; -noise_rois.tag = 'noise_rois'; -noise_rois.name = 'Noise ROIs model (Principal components of anatomical regions)'; -noise_rois.val = {noise_rois_no}; -noise_rois.values = {noise_rois_no, noise_rois_yes}; -noise_rois.help = {'Noise ROIs model (Principal components of anatomical regions), similar to aCompCor, Behzadi et al. 2007'}; - - - -%-------------------------------------------------------------------------- -%% Sub-substructure Movement Model -%-------------------------------------------------------------------------- - -%-------------------------------------------------------------------------- -% movement_file_realignment_parameters -%-------------------------------------------------------------------------- - -movement_file_realignment_parameters = cfg_files; -movement_file_realignment_parameters.tag = 'file_realignment_parameters'; -movement_file_realignment_parameters.name = 'Realignment Parameter File'; -movement_file_realignment_parameters.val = {{''}}; -movement_file_realignment_parameters.help = {'...'}; -movement_file_realignment_parameters.filter = '.*'; -movement_file_realignment_parameters.ufilter = '.mat$|.txt$'; -movement_file_realignment_parameters.num = [0 1]; - - -%-------------------------------------------------------------------------- -% movement_order -%-------------------------------------------------------------------------- -movement_order = cfg_menu; -movement_order.tag = 'order'; -movement_order.name = 'order'; -movement_order.help = {'Order of movement regressors 6/12/24, including derivatives and squared parameters/derivatives'}; -movement_order.labels = {'6' '12' '24'}; -movement_order.values = {6, 12 24}; -movement_order.val = {6}; - - -%-------------------------------------------------------------------------- -% movement_outlier_translation_mm -%-------------------------------------------------------------------------- - -movement_outlier_translation_mm = cfg_entry; -movement_outlier_translation_mm.tag = 'outlier_translation_mm'; -movement_outlier_translation_mm.name = 'Outlier Translation Threshold (mm)'; -movement_outlier_translation_mm.help = { - 'Threshold, above which a stick regressor is created for ' - 'corresponding volume of exceeding shift' - '' - 'Set to Inf to switch off regressor creation' - }; -movement_outlier_translation_mm.strtype = 'e'; -movement_outlier_translation_mm.num = [1 1]; -movement_outlier_translation_mm.val = {Inf}; - - -%-------------------------------------------------------------------------- -% movement_outlier_rotation_deg -%-------------------------------------------------------------------------- - -movement_outlier_rotation_deg = cfg_entry; -movement_outlier_rotation_deg.tag = 'outlier_rotation_deg'; -movement_outlier_rotation_deg.name = 'Outlier Rotation Threshold (degrees)'; -movement_outlier_rotation_deg.help = { - 'Threshold, above which a stick regressor is created for ' - 'corresponding volume of exceeding rotational movement' - '' - 'Set to Inf to switch off regressor creation' - }; -movement_outlier_rotation_deg.strtype = 'e'; -movement_outlier_rotation_deg.num = [1 1]; -movement_outlier_rotation_deg.val = {Inf}; - - -%-------------------------------------------------------------------------- -% movement_no -%-------------------------------------------------------------------------- - -movement_no = cfg_branch; -movement_no.tag = 'no'; -movement_no.name = 'No'; -movement_no.val = {}; -movement_no.help = {'Movement regressors not used.'}; - - -%-------------------------------------------------------------------------- -% movement_yes -%-------------------------------------------------------------------------- - -movement_yes = cfg_branch; -movement_yes.tag = 'yes'; -movement_yes.name = 'Yes'; -movement_yes.val = {movement_file_realignment_parameters, movement_order, ... - movement_outlier_translation_mm movement_outlier_rotation_deg}; -movement_yes.help = {'Include Movement Model, as described in Friston et al., 1996.'}; - - - -%-------------------------------------------------------------------------- -% movement -%-------------------------------------------------------------------------- - -movement = cfg_choice; -movement.tag = 'movement'; -movement.name = 'Movement'; -movement.val = {movement_yes}; -movement.values = {movement_no, movement_yes}; -movement.help = {'Movement Model, as described in Friston et al., 1996'}; - - -%-------------------------------------------------------------------------- -%% Sub-substructure Other (model) -%-------------------------------------------------------------------------- - -%-------------------------------------------------------------------------- -% input_other_multiple_regressors -%-------------------------------------------------------------------------- - -input_other_multiple_regressors = cfg_files; -input_other_multiple_regressors.tag = 'input_multiple_regressors'; -input_other_multiple_regressors.name = 'input_multiple_regressors'; -input_other_multiple_regressors.val = {{''}}; -input_other_multiple_regressors.help = {'...'}; -input_other_multiple_regressors.filter = '.*'; -input_other_multiple_regressors.ufilter = '.mat$|.txt$'; -input_other_multiple_regressors.num = [0 Inf]; - - -%-------------------------------------------------------------------------- -% other_model_no -%-------------------------------------------------------------------------- - -other_model_no = cfg_branch; -other_model_no.tag = 'no'; -other_model_no.name = 'No'; -other_model_no.val = {}; -other_model_no.help = {'Movement regressors not used'}; - - -%-------------------------------------------------------------------------- -% other_model_yes -%-------------------------------------------------------------------------- - -other_model_yes = cfg_branch; -other_model_yes.tag = 'yes'; -other_model_yes.name = 'Yes'; -other_model_yes.val = {input_other_multiple_regressors}; -other_model_yes.help = {'Include Other multiple regressor file(s)'}; - - - -%-------------------------------------------------------------------------- -% other_model -%-------------------------------------------------------------------------- - -other_model = cfg_choice; -other_model.tag = 'other'; -other_model.name = 'Other Multiple Regressors'; -other_model.val = {other_model_no}; -other_model.values = {other_model_no, other_model_yes}; -other_model.help = {'Other multiple regressor file(s)'}; - - -%-------------------------------------------------------------------------- -%% Sub-structure model -%-------------------------------------------------------------------------- -model = cfg_branch; -model.tag = 'model'; -model.name = 'model'; -model.val = {output_multiple_regressors, output_physio, orthog, retroicor, ... - rvt, hrv, noise_rois, movement, other_model}; -model.help = {['Physiological Model to be estimated and Included in GLM as ' ... - 'multiple_regressors.txt']}; - - - - -%========================================================================== -%% Sub-structure verbose -%========================================================================== - -%-------------------------------------------------------------------------- -% level -%-------------------------------------------------------------------------- -level = cfg_entry; -level.tag = 'level'; -level.name = 'level'; -level.help = {'...'}; -level.strtype = 'e'; -level.num = [Inf Inf]; -level.val = {2}; - -%-------------------------------------------------------------------------- -% fig_output_file -%-------------------------------------------------------------------------- -fig_output_file = cfg_entry; -fig_output_file.tag = 'fig_output_file'; -fig_output_file.name = 'fig_output_file'; -fig_output_file.help = {'file name where figures are saved to;' - 'supported figure formats(via filename-suffix): jpg, png, fig, ps' - 'leave empty to not save output figures'}; -fig_output_file.strtype = 's'; -fig_output_file.num = [0 Inf]; -fig_output_file.val = {''}; - - - - -%-------------------------------------------------------------------------- -% use_tabs -%-------------------------------------------------------------------------- -use_tabs = cfg_menu; -use_tabs.tag = 'use_tabs'; -use_tabs.name = 'use_tabs'; -use_tabs.help = {'use spm_tabs for plotting'}; -use_tabs.labels = {'true' 'false'}; -use_tabs.values = {true, false}; -use_tabs.val = {false}; - - -%-------------------------------------------------------------------------- -% verbose -%-------------------------------------------------------------------------- -verbose = cfg_branch; -verbose.tag = 'verbose'; -verbose.name = 'verbose'; -verbose.help = { - ' determines how many figures shall be generated to follow the workflow' - ' of the toolbox and whether the graphical output shall be saved (to a' - ' PostScript-file)' - ' 0 = no graphical output;' - ' 1 = (default) main plots : Fig 1: gradient scan timing (if selected) ;' - ' Fig 2: heart beat/breathing statistics & outlier;' - ' Fig 3: final multiple_regressors matrix' - ' 2 = debugging plots for setting up new study or if Fig 2 had' - ' outliers' - ' Fig 1: raw phys logfile data' - ' Fig 2: gradient scan timing (if selected)' - ' Fig 3: cutout interval of logfile for' - ' regressor creation (including scan timing' - ' and raw phys data)' - ' Fig 4: heart beat/breathing statistics & outlier;' - ' Fig 5: time course of all sampled RETROICOR' - ' regressors' - ' Fig 6: final multiple_regressors matrix' - '' - ' 3 = all plots' - ' Fig 1: raw phys logfile data' - ' Fig 2: gradient scan timing (if selected)' - ' Fig 3: Slice assignment to volumes' - ' Fig 4: cutout interval of logfile for' - ' regressor creation (including scan timing' - ' and raw phys data)' - ' Fig 5: heart beat/breathing statistics & outlier;' - ' Fig 6: cardiac phase data of all slices' - ' Fig 7: respiratory phase data and' - ' histogram transfer function' - ' Fig 8: time course of all sampled RETROICOR' - ' regressors' - ' Fig 9: final multiple_regressors matrix' - - }; -verbose.val = {level fig_output_file use_tabs}; - - - - -%========================================================================== -%% Structure physio Assemblance -%========================================================================== - - -%-------------------------------------------------------------------------- -% physio -%-------------------------------------------------------------------------- -physio = cfg_exbranch; -physio.tag = 'physio'; -physio.name = 'TAPAS PhysIO Toolbox'; -physio.val = {save_dir files scan_timing preproc model verbose}; -physio.help = {'...'}; -physio.prog = @run_physio; -physio.vout = @vout_physio; - - -%========================================================================== -% function out = run_physio(job) -%========================================================================== -function out = run_physio(job) - -%% Rename job fields to the ones actually used in physio-structure - -physio = tapas_physio_job2physio(job); - -[physio, R] = tapas_physio_main_create_regressors(physio); - -out.physnoisereg = cellstr(physio.model.output_multiple_regressors); -out.R = R; -out.physiofile = cellstr(physio.model.output_physio); - - -%========================================================================== -% function dep = vout_physio(job) -%========================================================================== -function dep = vout_physio(job) -dep(1) = cfg_dep; -dep(1).sname = 'physiological noise regressors file (Multiple Regresssors)'; -dep(1).src_output = substruct('.','physnoisereg'); -dep(1).tgt_spec = cfg_findspec({{'filter','mat','strtype','e'}}); - -dep(2) = cfg_dep; -dep(2).sname = 'PhysIO Structure File'; -dep(2).src_output = substruct('.','physiofile'); -dep(2).tgt_spec = cfg_findspec({{'filter','mat','strtype','e'}}); diff --git a/PhysIO/code/tapas_physio_check_efficacy.m b/PhysIO/code/tapas_physio_check_efficacy.m deleted file mode 100644 index 55434c98..00000000 --- a/PhysIO/code/tapas_physio_check_efficacy.m +++ /dev/null @@ -1,179 +0,0 @@ -% This script reports all relevant F-contrast-maps for physIO-created regressors -% for the specified subjects -% -% Author: Lars Kasper -% Created: 2014-01-21 -% Copyright (C) 2014 TNU, Institute for Biomedical Engineering, University of Zurich and ETH Zurich. -% -% This file is part of the TNU CheckPhysRETROICOR toolbox, which is released under the terms of the GNU General Public -% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL -% (either version 3 or, at your option, any later version). For further details, see the file -% COPYING or . -% -% $Id$ - -%% ======================================================================== -% START #MOD - -% general paths study -pathSPM = '~/Documents/code/matlab/spm12b'; -pathPhysIO = '~/Documents/code/matlab/spm12b/toolbox/PhysIO'; -fileReport = ['~/Dropbox/Andreiuta/physio_rest_ioio_pharm/' ... - 'physio_IOIO_pharm/results/PhysIOTest_cardiac_overview_inferior_sliceParallel.ps']; % where contrast maps are saved - -% logfile names sorted per session -nSess = 1; - -% subject directories to be included into analysis - -% root directory holding all subject directories -pathDataRoot = '/Users/kasperla/Dropbox/Andreiuta/physio_rest_ioio_pharm/physio_IOIO_pharm/glmAnalysis'; - -% prefix of all subject directories -maskSubjects = 'DMPAD_*'; - -% GLM analysis subdirectory of subject folder -dirGLM = ''; - -% includes subdirectory of subject folder and file name mask -maskStructural = '^meanfunct_rest.nii'; - -% names of physiological contrasts to be reported -% namesPhysContrasts = { -% 'All Phys Regressors' -% 'Cardiac Regressors' -% 'Respiratory Regressors' -% 'Cardiac X Respiratory Interaction' -% 'Movement Regressors' -% }; - -namesPhysContrasts = { - 'All Phys' - 'Cardiac' - 'Respiratory' - 'Card X Resp Interation' - 'Movement' - }; - -% selection of physiological contrasts to be reported, corresponding to -% namesPhysContrasts order -indReportPhysContrasts = 2; - - -reportContrastThreshold = 0.001; % 0.05; 0.001; -reportContrastCorrection = 'none'; % 'FWE'; 'none'; -%reportContrastPosition = [0 -15 -2*16]; 'max'; % 'max' to jump to max; or [x,y,z] in mm -%fovMillimeter = 50; %mm; choose 0 to plot whole FOV (bounding box) -reportContrastPosition = [0 0 -30]; 'max'; % 'max' to jump to max; or [x,y,z] in mm -fovMillimeter = 0; %mm; choose 0 to plot whole FOV (bounding box) - -% if true, voxel space (parallel to slices), not world space (with interpolation) is used -doPlotSliceParallel = true; - -physio = tapas_physio_new('RETROICOR'); -model = physio.model; % holding number of physiological regressors - -% END #MOD -%% ======================================================================== - -scans = dir(fullfile(pathDataRoot,maskSubjects)); -scans = {scans.name}; -subjectIndices = 1:length(scans); - -delete(fileReport); -addpath(pathPhysIO); -addpath(pathSPM); -spm('defaults', 'fMRI'); -spm_jobman('initcfg'); - -for s = subjectIndices - - try - dirSubject = scans{s}; - pathSubject = fullfile(pathDataRoot,dirSubject); %dirSubject = scan - pathAnalysis = fullfile(pathDataRoot,dirSubject,dirGLM); - fileSPM = fullfile(pathAnalysis, 'SPM.mat'); - fileStruct = spm_select('ExtFPList', pathSubject, maskStructural); - - - - %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % Create and plot phys regressors F-contrasts - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - if exist(fileSPM, 'file') - load(fileSPM); - nContrasts = numel(indReportPhysContrasts); - - - % Check whether contrasts already exist in SPM.mat - indContrasts = zeros(nContrasts,1); - for c = 1:nContrasts - iC = indReportPhysContrasts(c); - indContrasts(c) = tapas_physio_check_get_xcon_index(SPM, ... - namesPhysContrasts{iC}); - end - - % Generate contrasts only if not already existing - matlabbatch = tapas_physio_check_prepare_job_contrasts(fileSPM, ... - model, SPM, indReportPhysContrasts, pathPhysIO); - matlabbatch{1}.spm.stats.con.consess(find(indContrasts)) = []; - if ~isempty(matlabbatch{1}.spm.stats.con.consess) - spm_jobman('run', matlabbatch); - load(fileSPM); - end - - % report contrasts - for c = 1:nContrasts - iC = indReportPhysContrasts(c); - indContrasts(c) = tapas_physio_check_get_xcon_index(SPM, ... - namesPhysContrasts{iC}); - load(fullfile(pathPhysIO, 'tapas_physio_check_job_report')); - matlabbatch{1}.spm.stats.results.spmmat = cellstr(fileSPM); - matlabbatch{1}.spm.stats.results.conspec.titlestr = [dirSubject ' - ' namesPhysContrasts{iC}]; - matlabbatch{1}.spm.stats.results.conspec.contrasts = indContrasts(c); - - % contrast report correction - matlabbatch{1}.spm.stats.results.conspec.thresh = reportContrastThreshold; - matlabbatch{1}.spm.stats.results.conspec.threshdesc = reportContrastCorrection; - - spm_jobman('run', matlabbatch); % report result - % spm_print(fileReport) - spm_sections(xSPM,hReg, fileStruct); % overlay structural - - % voxel, not world space - if doPlotSliceParallel - spm_orthviews('Space',1) - end - - spm_orthviews('Zoom', fovMillimeter); % zoom to FOV*2 view - spm_orthviews('Interp', 0); % next neighbour interpolation plot - - - if isequal(reportContrastPosition, 'max'); - spm_mip_ui('Jump',spm_mip_ui('FindMIPax'),'glmax'); % goto global max - else - spm_mip_ui('SetCoords', reportContrastPosition, ... - spm_mip_ui('FindMIPax')); % goto global max - end - - - % spm_print always prepend current directory to print-file - % name :-( - [pathReport, filenameReport] = fileparts(fileReport); - pathTmp = pwd; - cd(pathReport); - spm_print(filenameReport); - cd(pathTmp); - end - - titstr = [dirSubject, ' - SPM.xX.X']; - title(regexprep(titstr,'_','\\_')); - set(gcf,'Name', titstr); - fprintf('good SPM: %s\n', dirSubject); - else % no file, report that - fprintf('no SPM.mat: %s\n', dirSubject); - end - catch - warning(sprintf('Subject ID %d: %s did not run through', s, dirSubject)); - end -end \ No newline at end of file diff --git a/PhysIO/code/tapas_physio_check_get_regressor_columns.m b/PhysIO/code/tapas_physio_check_get_regressor_columns.m deleted file mode 100644 index 4efe9dcf..00000000 --- a/PhysIO/code/tapas_physio_check_get_regressor_columns.m +++ /dev/null @@ -1,109 +0,0 @@ -function [colPhys, colCard, colResp, colMult, colHRV, colRVT, colRois, colMove, colAll] = ... - tapas_physio_check_get_regressor_columns(SPM, model) -% -% returns indices of physiological regressors in an SPM design-matrix, -% pooled over all sessions for later creation of an F-contrast -% -% INPUT: -% SPM SPM.mat -% model physIO.model-structure -% -% -% OUTPUT: -% colPhys - index vector of all physiological regressor columns in design matrix (from SPM.xX.names) -% colCard - index vector of cardiac regressor columns in design matrix (from SPM.xX.names) -% colResp - index vector of respiratory regressor columns in design matrix (from SPM.xX.names) -% colMult - index vector of interaction cardiac X respiration regressor columns in design matrix (from SPM.xX.names) -% colHRV - index vector of heart rate variability column in design matrix (from SPM.xX.names) -% colRVT - index vector of respiratory volume per time column in design matrix (from SPM.xX.names) -% colRois - index vector of noise rois regressors (from SPM.xX.names) -% colMove - index vector of movement regressor columns in design matrix (from SPM.xX.names) -% colAll - colPhys and colMove (all nuisance regressors!) -% -% Author: Lars Kasper -% Created: 2014-01-21 -% Copyright (C) 2014 TNU, Institute for Biomedical Engineering, University of Zurich and ETH Zurich. -% -% This file is part of the TNU CheckPhysRETROICOR toolbox, which is released under the terms of the GNU General Public -% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL -% (either version 3 or, at your option, any later version). For further details, see the file -% COPYING or . -% -% $Id$ - -nSess = length(SPM.Sess); - -iCard = 1; - -if nargin < 2 - nMove = 6; - nCard = 6; - nResp = 8; - nMult = 4; - nHRV = 0; - nRVT = 0; - nRois = 0; -else - - if model.retroicor.include - - nCard = ~isempty(model.retroicor.order.c); - if nCard - nCard = model.retroicor.order.c*2; - end - - nResp = ~isempty(model.retroicor.order.r); - if nResp - nResp = model.retroicor.order.r*2; - end - - nMult = ~isempty(model.retroicor.order.cr); - if nMult - nMult = model.retroicor.order.cr*4; - end - - - else - nCard = 0; - nResp = 0; - nMult = 0; - end - - - % only if following models were calculated - nRois = model.noise_rois.include.*sum(model.noise_rois.n_components); - nMove = model.movement.include.*model.movement.order; - nHRV = model.hrv.include.*numel(model.hrv.delays); - nRVT = model.rvt.include.*numel(model.rvt.delays); - -end - -cnames = SPM.xX.name'; - -colCard = []; -colResp = []; -colMult = []; -colRois = []; -colMove = []; -colHRV = []; -colRVT = []; -for s = 1:nSess - - cname = ['Sn(' int2str(s) ') R' int2str(iCard)]; - indC = find(strcmp(cnames, cname)); - - colCard = [colCard, indC:(indC+nCard - 1)]; - colResp = [colResp, (indC+nCard):(indC+nCard+nResp - 1)]; - colMult = [colMult, (indC+nCard+nResp):(indC+nCard+nResp+nMult - 1)]; - colHRV = [colHRV, (indC+nCard+nResp+nMult):... - (indC+nCard+nResp+nMult+nHRV-1)]; - colRVT = [colRVT, (indC+nCard+nResp+nMult+nHRV):... - (indC+nCard+nResp+nMult+nHRV+nRVT-1)]; - colRois = [colRois, (indC+nCard+nResp+nMult+nHRV+nRVT):... - (indC+nCard+nResp+nMult+nHRV+nRVT+nRois-1)]; - colMove = [colMove, (indC+nCard+nResp+nMult+nHRV+nRVT+nRois):... - (indC+nCard+nResp+nMult+nHRV+nRVT+nRois+nMove-1)]; -end - -colPhys = [colCard colResp colMult colHRV colRVT colRois]; -colAll = [colPhys colMove]; \ No newline at end of file diff --git a/PhysIO/code/tapas_physio_check_get_xcon_index.m b/PhysIO/code/tapas_physio_check_get_xcon_index.m deleted file mode 100644 index 56c4b376..00000000 --- a/PhysIO/code/tapas_physio_check_get_xcon_index.m +++ /dev/null @@ -1,37 +0,0 @@ -function indC = tapas_physio_check_get_xcon_index(SPM, cname) -% returns index of contrasts in SPM.xCon, whose name matches cname -% used to find existing contrasts in an SPM.mat -% -% INPUT: -% SPM - SPM.mat-file -% cname - name of contrast to be searched for (string) -% (works for incomplete names/patterns as well, e.g. to find all F-contrasts) -% -% OUTPUT: -% indC - vector of contrast indices in SPM.xCon with matching name -% returns 0, if none found -% -% Author: Lars Kasper -% Created: 2014-01-21 -% Copyright (C) 2014 TNU, Institute for Biomedical Engineering, University of Zurich and ETH Zurich. -% -% This file is part of the TNU CheckPhysRETROICOR toolbox, which is released under the terms of the GNU General Public -% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL -% (either version 3 or, at your option, any later version). For further details, see the file -% COPYING or . -% -% $Id$ -nContrasts = length(SPM.xCon); -cnames = cell(nContrasts,1); -for c = 1:nContrasts - cnames{c} = SPM.xCon(c).name; -end -indC = find(cell2mat(cellfun(@(x) ~isempty(x), strfind(cnames, cname), 'UniformOutput', false))); - -if isempty(indC) - indC = 0; -else - % if multiple regressors of same name found, take first one - indC = indC(1); -end -end \ No newline at end of file diff --git a/PhysIO/code/tapas_physio_check_job_contrasts.mat b/PhysIO/code/tapas_physio_check_job_contrasts.mat deleted file mode 100644 index 185f1638b70bdc575eb3b6b771b748d9e0caef68..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 40848 zcmeI*Pj4GV7{~D^Rf{TB>H#DOQfapmLL5Sqm;^;R;Z~GHs-p-hh`7MU*^sEVV|iDl zz3~!w0sOh}5_kz5cnQ1!4jeghbud{JF~OT?tW&Z*qPZ)r~U9j z`);y+XC>*hAAWFn+{>1d-CkN8o{ajN$)^u?_Szq@42lJpR^7sDIfDwM@3pSIYbO z`{LVwQ(CsdG8^e9)%tgKqP0hs+fezgwk&J?*uJ%%)_pD2>$2@@`?7WHGutj}bX7jk zmaog2eU*K#_4~G?*9}{@)w0&m`p;^#6~iOt`st}O*KvEwALz3{PIP0qY z)##&aSYFSb%BQv88`>VO)7$SorMA?1mdNSIo=P!Mn`CN*QQ2O3*`A@lUrnm`AFS`MS?z7-&YpKywRx|D(67{Vy_~BB>dtsreMIN7whR^7k~uWMl{BOqZ^` zXg0-o{>l5Az!Us}Q= zLq{Oo`5_y&JO7EWCjtl{fB*srAb( z#yRVF%QrF#;~_p?Ba0zuWj&qb`1LuitSACJ3yq3xbReRuuqhF(|@KmY**5I_I{1Q3`i5I#XptyF&_VPyuw>=5dA zN%lYh0R#|0009ILKmY**5NJ|>@4V6kVa5nej7uNguEt|8GA^CejC0m;=_fLrHZHBh zKI>>cE`3)iz7KHE#$(dY8>`;d9g*XwX!|p>3yo#_J(Ahf;E*84C!BH@Q!P*WDCl!2X;We*=y6>TS$gyYA%6j0 zqRaIo*#$1C6Gq6XDXxf5Jo6MI3~_+WM^_eid3V^ny3Td_sCTA1j;StjsjXbeEnDfz zxKbP4oO7C=+iJGj@))q6o%PMtGp#PYuU<6IfLiC#?q7{-_D{YxPkc$U@8Z. -% -% $Id$ - -if nargin < 6 - namesPhysContrasts = tapas_physio_get_contrast_names_default(); -end - -if ~exist('SPM', 'var'), load(fileSPM); end - -[colPhys, colCard, colResp, colMult, colHRV, colRVT, colRois, colMove, colAll] = ... - tapas_physio_check_get_regressor_columns(SPM, model); -con{1} = colPhys; -con{2} = colCard; -con{3} = colResp; -con{4} = colMult; -con{5} = colHRV; -con{6} = colRVT; -con{7} = colRois; -con{8} = colMove; -con{9} = colAll; - -pathCodePhysIO = fileparts(mfilename('fullpath')); -load(fullfile(pathCodePhysIO,'tapas_physio_check_job_contrasts.mat')); -matlabbatch{1}.spm.stats.con.spmmat{1} = fileSPM; - -% execute computation only for non-empty contrasts; -indValidContrasts = intersect(indReportPhysContrasts, ... - find(~cellfun(@isempty, con))); -nContrasts = numel(indValidContrasts); -for c = 1:nContrasts - iC = indValidContrasts(c); - F{c} = zeros(max(con{iC})); - F{c}(sub2ind(size(F{c}),con{iC}, con{iC})) = 1; - matlabbatch{1}.spm.stats.con.consess{c}.fcon.convec{1} = F{c}; - - matlabbatch{1}.spm.stats.con.consess{c}.fcon.name = ... - namesPhysContrasts{indValidContrasts(c)}; -end - -% remove additional contrasts in matlabbatch -matlabbatch{1}.spm.stats.con.consess(nContrasts+1:end) = []; \ No newline at end of file diff --git a/PhysIO/code/tapas_physio_compute_tsnr_gains.m b/PhysIO/code/tapas_physio_compute_tsnr_gains.m deleted file mode 100644 index 66199584..00000000 --- a/PhysIO/code/tapas_physio_compute_tsnr_gains.m +++ /dev/null @@ -1,135 +0,0 @@ -function [tSnrImageArray, fileTsnrArray, ... - tSnrRatioImageArray, fileTsnrRatioArray ] = ... - tapas_physio_compute_tsnr_gains(physio, SPM, indexContrastForSnrRatio, ... - namesPhysContrasts) -% Computes tSNR gains through physiological noise correction for all -% confound regressor sets modelled, after estimation of the SPM-GLM -% -% [fileTsnrGainArray, fileTsnrArray] = ... -% tapas_physio_compute_tsnr_gains(physio, SPM, doSave); -% -% -% IN -% physio physio-structure (or filename with structure), after estimating multiple_regressors.txt -% SPM SPM variable (or filename with saved structure) -% after estimation of general linear model (GLM) -% indexContrastForSnrRatio -% contrast id for comparison to get relative tSNR values -% default: 0 (raw time series after filtering/pre-whitening) -% OUT -% fileTsnrGainArray -% cell(nPhysioSets,1) of nii-filenames holding tSNR gain -% images for all physiological regressor sets (in percent?) -% -% fileTsnrArray -% cell(nPhysioSets+1,1) of nii-filenames for tSNR images of -% for all physiological regressor sets -% and, as last element, raw tSNR (after preprocessing) -% EXAMPLE -% tapas_physio_compute_tsnr_gains -% -% See also tapas_physio_compute_tsnr_spm spm_write_residuals -% -% Author: Lars Kasper -% Created: 2015-07-03 -% Copyright (C) 2015 TNU, Institute for Biomedical Engineering, -% University of Zurich and ETH Zurich. -% -% This file is part of the TAPAS PhysIO Toolbox, which is released under the terms of the GNU General Public -% License (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL -% (either version 3 or, at your option, any later version). For further details, see the file -% COPYING or . -% -% $Id$ - - - -%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% 0. Check variable structure, cast for convenience -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -if nargin < 3 - indexContrastForSnrRatio = 0; % contrast for comparison -end - -if nargin < 4 - namesPhysContrasts = tapas_physio_get_contrast_names_default(); -end - -% load physio-variable, if filename given -if ~isstruct(physio) - % load SPM variable from file - if iscell(physio) - filePhysio = physio{1}; - else - filePhysio = physio; - end - load(filePhysio, 'physio'); -end - -% load SPM-variable, if filename given -if ~isstruct(SPM) - % load SPM variable from file - if iscell(SPM) - fileSpm = SPM{1}; - else - fileSpm = SPM; - end - load(fileSpm, 'SPM'); -end - - -%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% 1. Determine which physiological contrasts could be created with existing -% model and create the missing ones -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -SPM = tapas_physio_create_missing_physio_contrasts(SPM, physio.model, ... - namesPhysContrasts); - -nContrasts = numel(namesPhysContrasts); - - -%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% 2. Compute tSNR maps via tapas_physio_compute_tsnr_spm, as follows -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -tSnrImageArray = cell(nContrasts,1); -fileTsnrArray = cell(nContrasts,1); -tSnrRatioImageArray = cell(nContrasts,1); -fileTsnrRatioArray = cell(nContrasts,1); - -for c = 1:nContrasts - % First try, whether inverted contrast already exists, then no - % additional computation has to be performed! - indContrast = tapas_physio_check_get_xcon_index(SPM, ... - ['All but: ' namesPhysContrasts{c}]); - doComputeInvertedContrast = indContrast == 0; % not found, do invert - doSaveNewContrasts = doComputeInvertedContrast; % and save contrasts for later - if doComputeInvertedContrast % find contrast before inversion - indContrast = tapas_physio_check_get_xcon_index(SPM, ... - namesPhysContrasts{c}); - end - - if indContrast > 0 - % if contrast exist, compute tSNR and save it! - [tSnrImageArray{c}, fileTsnrArray{c}, ... - tSnrRatioImageArray{c}, fileTsnrRatioArray{c}] = ... - tapas_physio_compute_tsnr_spm(SPM, indContrast, ... - indexContrastForSnrRatio, doComputeInvertedContrast, doSaveNewContrasts); - - % rename files for better human readability - filenameSfx = regexprep(namesPhysContrasts{c}, ' ', '_'); - fileTsnrOld = fileTsnrArray{c}; - fileTsnrArray{c} = regexprep(fileTsnrArray{c}, ... - sprintf('con%04d', indContrast), filenameSfx); - - fileTsnrRatioOld = fileTsnrRatioArray{c}; - fileTsnrRatioArray{c} = regexprep(fileTsnrRatioArray{c}, ... - sprintf('con%04dvs%04d', indContrast, indexContrastForSnrRatio), ... - [filenameSfx '_vs_Raw']); - - movefile(fileTsnrOld, fileTsnrArray{c}); - movefile(fileTsnrRatioOld, fileTsnrRatioArray{c}); - end -end \ No newline at end of file diff --git a/PhysIO/code/tapas_physio_compute_tsnr_spm.m b/PhysIO/code/tapas_physio_compute_tsnr_spm.m deleted file mode 100644 index 47382bc6..00000000 --- a/PhysIO/code/tapas_physio_compute_tsnr_spm.m +++ /dev/null @@ -1,217 +0,0 @@ -function [tSnrImage, fileTsnr, tSnrRatioImage, fileTsnrRatio] = ... - tapas_physio_compute_tsnr_spm(SPM, iC, iCForRatio, doInvert, doSaveNewContrasts) -% Computes temporal SNR image after correcting for a contrast -% from SPM general linear model -% -% tSnrImage = compute_tsnr_spm(input) -% -% This function executes the following steps: -% 1. Compute F-Contrasts of the kind "All-but-selected regressors" -% run spm_write_residuals to retrieve Y-Yc = Y0 + e -% (note that roles of Y0 and Yc are swapped compared to in spm_spm, -% since we are interested in Residuals that do *not* stem from our -% regressors, but the other ones! -% 2. Compute tSNR images using mean(Residuals)/std(Residuals) from -% F-contrasts from (1.), save as nii-files tsnr_con.nii -% 3. Compute tSNR for contrast of comparison, default: -% preprocessed time-series, after pre-whitening and -% highpass-filtering (i.e. K*W*Y), saved as nii-files -% using a recursive call to this function -% 4. Compute tSNR ratio images from (3.) and (2.), save as nii-files -% tsnr_convs.nii -% -% IN -% SPM SPM variable (in SPM.mat, or file name) after parameter and -% contrast estimation -% iC Contrast of interest for which tSNR is computed -% iC = 0 computes tSNR of the raw image time series after -% pre-whitening and filtering, but without any model -% confound regressors removed. (default) -% iC = NaN computes tSNR after removing the full model (including -% the mean), and is therefore equivalent to sqrt(1/ResMS) -% iCForRatio -% For computation of relative tSNR. -% Contrast index whose tSNR image should be used as a -% denominator for tSnrRatioImage = tSNRiC./tSNRiCForRatio; -% tSNR is compared to the tSNR image of the specified contrast -% (0 for raw tSNR). The corresponding tSNR-image will be created -% default = 0 (raw tSNR); Leave [] to not compute ratio of tSNRs -% doInvert true (default for iC > 0 and ~nan) or false -% typically, one is interested in the tSNR *after* correcting for -% the contrast in question. To compute this, one has to look at -% the residuals of the inverse F-contrast that excludes all but -% the contrast of interest, i.e. eye(nRegressors) - xcon -% Thus, this function computes this contrast per default and goes -% from there determining residuals etc. -% doSaveNewContrasts -% true or false (default) -% if true, the temporary contrasts created by inversion of -% selected columns will be saved into the SPM structure. -% Otherwise, all temporary data will be removed. -% -% OUT -% tSnrImage [nX,nY,nZ] image matrix holding tSNR when only including regressors in -% contrast iC in design matrix; -% i.e. gives the effect of -% mean(Xc*bc + e)/std(Xc*bc + e) -% if Xc hold only regressors of iC -% -% NOTE: If you want to estimate the effect of a noise -% correction using some regressors, the Contrast with index -% iC should actually contain an F-contrast *EXCLUDING* these -% regressors, because everything else constitutes the -% regressors of interest -% fileTsnr tsnr_con.nii -% path and file name of tSNR image, if doSave was true -% tSnrRatioImage -% tSNR_con./tSNR_con -% (if iCForRatio was specified) -% fileTsnrRatio -% tsnr_convs.nii -% file where tSnrRatio was saved -% -% EXAMPLE -% compute_tsnr_spm -% -% See also spm_write_residuals spm_FcUtil -% -% Author: Lars Kasper -% Created: 2014-12-10 -% Copyright (C) 2014 Institute for Biomedical Engineering, ETH/Uni Zurich. -% $Id$ - -if nargin < 2 - iC = 0; -end - -if nargin < 3 - iCForRatio = 0; -end - -if nargin < 4 - doInvert = true; -end - -if nargin < 5 - doSaveNewContrasts = false; -end - - -doComputeRatio = ~isempty(iCForRatio); - - -% load SPM-variable, if filename given -if ~isstruct(SPM) - % load SPM variable from file - if iscell(SPM) - fileSpm = SPM{1}; - else - fileSpm = SPM; - end - load(fileSpm, 'SPM'); -end - -oldDirSpm = SPM.swd; - -if ~doSaveNewContrasts - % temporary changes to SPM structure saved in sub-dir, removed later - newDirSpm = fullfile(SPM.swd, 'tmp'); - mkdir(newDirSpm); - copyfile(fullfile(SPM.swd, '*.nii'), newDirSpm); - copyfile(fullfile(SPM.swd, 'SPM.mat'), newDirSpm); - SPM.swd = newDirSpm; -end - -iCIn = iC; - -isInvertableContrast = iC > 0 && ~isnan(iC); - -if isInvertableContrast && doInvert - if ~isequal(SPM.xCon(iC).STAT, 'F') - error('Can only invert F-contrasts'); - end - - idxColumnsContrast = find(sum(SPM.xCon(iC).c)); - - Fc = spm_FcUtil('Set', ['All but: ' SPM.xCon(iC).name], 'F', ... - 'iX0', idxColumnsContrast, SPM.xX.xKXs); - SPM.xCon(end+1) = Fc; - SPM = spm_contrasts(SPM,numel(SPM.xCon)); - - % use this for computation - iC = numel(SPM.xCon); -end - - -%% Write residuals Y - Y0 = Yc + e; -VRes = spm_write_residuals(SPM, iC); -nVolumes = numel(VRes); -for iVol = 1:nVolumes - VRes(iVol).fname = fullfile(SPM.swd, VRes(iVol).fname); -end - -fileTsnr = fullfile(oldDirSpm, sprintf('tSNR_con%04d.nii', iCIn)); - -useMrImage = false; - -if useMrImage % use toolbox functionality - ResImage = MrImage(VRes(1).fname); - - % Create 4D image of contrast-specific "residuals", i.e. Xc*bc + e - for iVol = 1:nVolumes - ResImage.append(VRes(iVol).fname); - end - - % compute tSNR = mean(Xc*bc + e)/std(Xc*bc + e) - tSnrImage = ResImage.mean./ResImage.std; - if doSaveNewContrasts - tSnrImage.save(fileTsnr); - end -else - ResImage = spm_read_vols(VRes); - meanImage = mean(ResImage, 4); - stdImage = std(ResImage, 0, 4); - tSnrImage = meanImage./stdImage; - VTsnr = VRes(1); - VTsnr.fname = fileTsnr; - spm_write_vol(VTsnr, tSnrImage); -end - -%% -if doComputeRatio - fileTsnrCompare = fullfile(oldDirSpm, sprintf('tSNR_con%04d.nii', iCForRatio)); - - % Load or compute tSNR image for comparison contrast - if ~exist(fileTsnrCompare, 'file'); - % when computed, don't compute another ratio, and don't delete tmp - % here! (will be done at the end) - tSnrCompareImage = tapas_physio_compute_tsnr_spm(... - fullfile(oldDirSpm, 'SPM.mat'), ... - iCForRatio, doInvert, [], 1); - else - VCompare = spm_vol(fileTsnrCompare); - tSnrCompareImage = spm_read_vols(VCompare); - end - - % compute tSNR ratio and save - tSnrRatioImage = tSnrImage./tSnrCompareImage; - - fileTsnrRatio = fullfile(oldDirSpm, ... - sprintf('tSNRRatio_con%04dvs%04d.nii', iCIn, iCForRatio)); - VRatio = spm_vol(fileTsnrCompare); - VRatio.fname = fileTsnrRatio; - spm_write_vol(VRatio, tSnrRatioImage); -else - tSnrRatioImage = []; - fileTsnrRatio = []; -end - - -%% clean up all created residual files and temporary SPM folder -if ~doSaveNewContrasts - delete(fullfile(newDirSpm, '*')); - rmdir(newDirSpm); -else - % delete at least the Res-images - delete(fullfile(oldDirSpm, 'Res_*.nii')); -end \ No newline at end of file diff --git a/PhysIO/code/tapas_physio_corrcoef12.m b/PhysIO/code/tapas_physio_corrcoef12.m deleted file mode 100644 index 4ca9ec38..00000000 --- a/PhysIO/code/tapas_physio_corrcoef12.m +++ /dev/null @@ -1,80 +0,0 @@ -function [correlation,x,y] = tapas_physio_corrcoef12(x,y, isZtransformed) -% computes correlation coefficient (i.e. entry (1,2) of correlation matrix) -% quickly between two time series -% The speed-up is mainly due to in-line re-implementation of costly -% statistical functions, i.e. mean, std, cov, and usage of a pre-applied -% z-transformation of either of the inputs -% -% [correlation,x,y] = tapas_physio_corrcoef12(x,y, isZtransformed) -% -% -% IN -% x [nSamples,1] column vector of samples -% y [nSamples,1] column vector of samples -% isZtransformed [1,2] vector stating whether x,y or both are -% z-transformed, i.e. mean-corrected and divided by their -% standard deviation -% example: -% isZtransformed = [1,0] means that x is z-transformed, -% by y is not, i.e. (y-mean(y))/std(y) will be computed -% OUT -% correlation correlation(1,2) of correlation = corrcoef(x,y) -% x z-transformed x -% y z-transformed y -% EXAMPLE -% tapas_physio_corrcoef12 -% -% See also -% -% Author: Lars Kasper -% Created: 2014-08-14 -% Copyright (C) 2014 TNU, Institute for Biomedical Engineering, University of Zurich and ETH Zurich. -% -% This file is part of the TAPAS PhysIO Toolbox, which is released under the terms of the GNU General Public -% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL -% (either version 3 or, at your option, any later version). For further details, see the file -% COPYING or . -% -% $Id$ -if nargin < 3 - isZtransformed = [0 0]; -end - -% This is the old implementation; uncomment for comparison purposes -%doUseSlow = false; -%if doUseSlow -% correlation = corrcoef(x,y); -% correlation = correlation(1,2); -%else %fast, using shortcuts and in-line implementationf of mean/std/cov... - - -%C(i,j)/SQRT(C(i,i)*C(j,j)). - -correlation = 0; -if any(x) && any(y) % all-zero vectors result in zero correlation - - nSamples = numel(x); - normFactor = 1/(nSamples-1); - - % make column vectors - if size(x,1) ~= nSamples - x = x(:); - y = y(:); - end - - if ~isZtransformed(1) % perform z-transformation - x = x - sum(x)/nSamples; - x = x./sqrt(x'*x*normFactor); - end - if ~isZtransformed(2) % perform z-transformation - y = y - sum(y)/nSamples; - y = y./sqrt(y'*y*normFactor); - end - - if numel(x) == numel(y) - correlation = x'*y*normFactor; - % otherwise, correlation stays zero - end -end - -%end % else doUseSlow \ No newline at end of file diff --git a/PhysIO/code/tapas_physio_correct_cardiac_pulses_manually.m b/PhysIO/code/tapas_physio_correct_cardiac_pulses_manually.m deleted file mode 100644 index f77f2fe6..00000000 --- a/PhysIO/code/tapas_physio_correct_cardiac_pulses_manually.m +++ /dev/null @@ -1,127 +0,0 @@ -function [ons_secs, outliersHigh, outliersLow, verbose] = ... - tapas_physio_correct_cardiac_pulses_manually(ons_secs, thresh_cardiac, ... - verbose) -% this function takes the onsets from ECG measure and controls for -% outliers (more or less than a threshold given by a percentile increased -% or decreased by upperTresh or lowerThresh percent respectively. -% -% Author: Jakob Heinzle, TNU, -% adaptation: Lars Kasper -% -% Copyright (C) 2013, Institute for Biomedical Engineering, ETH/Uni Zurich. -% -% This file is part of the PhysIO toolbox, which is released under the terms of the GNU General Public -% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL -% (either version 3 or, at your option, any later version). For further details, see the file -% COPYING or . -% -% $Id$ - -if verbose.level < 1 - error('Manual picking of cardiac pulses requires verbose.level >= 1'); -end - -percentile = thresh_cardiac.percentile; -upperThresh = thresh_cardiac.upper_thresh; -lowerThresh = thresh_cardiac.lower_thresh; - -[outliersHigh,outliersLow, verbose] = tapas_physio_cardiac_detect_outliers(... - ons_secs.cpulse, percentile, upperThresh, lowerThresh, verbose); -if any(outliersHigh) - %disp('Press Enter to proceed to manual peak selection!'); - %pause; - additionalPulse=[]; - fh2=figure; - for outk=1:length(outliersHigh) - s=0; - while ~(s==1) - indStart = outliersHigh(outk)-1; indEnd = outliersHigh(outk)+2; - ind=find(ons_secs.t>=ons_secs.cpulse(indStart), 1, 'first')-100:find(ons_secs.t<=ons_secs.cpulse(indEnd), 1, 'last')+100; - figure(fh2); clf; - plot(ons_secs.t(ind),ons_secs.c(ind),'r') - hold on; - plot(ons_secs.cpulse(indStart:indEnd),ones(4,1)*max(ons_secs.c(ind)),'ok') - inpNum=input('How many triggers do you want to set? Enter a number between 0 and 10 : '); - I1=[]; - for ii=1:inpNum - figure(fh2); - [I1(ii), J1] = ginput(1); - plot(I1(ii),J1, 'b*', 'MarkerSize',10); - - end - - s=input('If you agree with the selected triggers, press 1 (then enter) : '); - if isempty(s) - s=0; - end - end - additionalPulse=[additionalPulse;I1']; - end - ons_secs.cpulse = sort([ons_secs.cpulse;additionalPulse]); - close(fh2); -end -close(verbose.fig_handles(end)); -verbose.fig_handles(end) = []; - -[outliersHigh,outliersLow, verbose] = tapas_physio_cardiac_detect_outliers(... - ons_secs.cpulse, percentile, upperThresh, lowerThresh, verbose); - -nPulses = length(ons_secs.cpulse); -finalIndex=1:nPulses; - -nWindow = 10; - -%show windows with suspicious outliers only -for indStart = round(nWindow/2):nWindow-2:nPulses - indEnd = indStart+nWindow-1; - if any(ismember(outliersLow, indStart:indEnd)) - %disp('Press Enter to proceed to manual peak deletion!'); - %pause; - fh3=figure('Position', [500 500 1000 500]); - ind=find(ons_secs.t>=ons_secs.cpulse(indStart), 1, 'first')-100:find(ons_secs.t<=ons_secs.cpulse(indEnd), 1, 'last')+100; - figure(fh3); clf; - plot(ons_secs.t(ind),ons_secs.c(ind),'r') - hold on; - plot(ons_secs.cpulse(indStart:indEnd),ones(nWindow,1)*max(ons_secs.c(ind)),'ko','MarkerFaceColor','r'); - alreadyDeleted=intersect(indStart:indEnd,setdiff(1:nPulses,finalIndex)); - plot(ons_secs.cpulse(alreadyDeleted),ones(size(alreadyDeleted))*max(ons_secs.c(ind)),'ro'); - for kk=indStart:indEnd - text(ons_secs.cpulse(kk),max(ons_secs.c(ind))*1.05,int2str(kk)); - end - - delInd= []; - - delInd=input('Enter the indices of pulses you want to delete (0 if none, put multiple in [], e.g. [19,20]): '); - - if ~isempty(delInd) && any(find(delInd~=0)) - plot(ons_secs.cpulse(delInd),max(ons_secs.c(ind))*ones(size(delInd)), 'rx', 'MarkerSize',20); - finalIndex=setdiff(finalIndex,delInd'); - end - close(fh3); - end - -end -ons_secs.cpulse = sort(ons_secs.cpulse(finalIndex)); - -close(verbose.fig_handles(end)); -verbose.fig_handles(end) = []; - -[outliersHigh,outliersLow, verbose] = tapas_physio_cardiac_detect_outliers(... - ons_secs.cpulse, percentile, upperThresh, lowerThresh, verbose); - -close(verbose.fig_handles(end)); -verbose.fig_handles(end) = []; - -% recursively determine outliers -if ~isempty(outliersHigh) || ~isempty(outliersLow) - doManualCorrectionAgain = input('More outliers detected after correction. Do you want to remove them? (1=yes, 0=no; ENTER)'); - if doManualCorrectionAgain - [ons_secs, outliersHigh, outliersLow] = ... - tapas_physio_correct_cardiac_pulses_manually(ons_secs, ... - thresh_cardiac, verbose); - end -end - -cpulse = ons_secs.cpulse; -save(thresh_cardiac.file, 'ons_secs', 'thresh_cardiac', 'cpulse'); -end diff --git a/PhysIO/code/tapas_physio_count_physio_regressors.m b/PhysIO/code/tapas_physio_count_physio_regressors.m deleted file mode 100644 index 15fd3a9c..00000000 --- a/PhysIO/code/tapas_physio_count_physio_regressors.m +++ /dev/null @@ -1,65 +0,0 @@ -function nPhysioRegressors = tapas_physio_count_physio_regressors(physio) -% Returns number of physiological regressors created, given model -% specification; -% NOTE: only reproducible numbers (data-independent) are -% returned, i.e. session-specific movement spikes and %-variance explained -% PCA-components are not included -% -% IN -% physio physio-structure, See also tapas_physio_new -% -% OUT -% nPhysioRegressors number of physiological regressors, e.g. motion, -% retroicor, noise_rois -% but: ignores -% -% EXAMPLE -% tapas_physio_report_contrasts -% -% See also -% -% Author: Lars Kasper -% Created: 2014-10-16 -% Copyright (C) 2014 TNU, Institute for Biomedical Engineering, University of Zurich and ETH Zurich. -% -% This file is part of the TAPAS PhysIO Toolbox, which is released under the terms of the GNU General Public -% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL -% (either version 3 or, at your option, any later version). For further details, see the file -% COPYING or . -% -% $Id: teditRETRO.m 464 2014-04-27 11:58:09Z kasperla $ - -model = physio.model; - -nPhysioRegressors = 0; - -if model.hrv.include - nPhysioRegressors = nPhysioRegressors + numel(model.hrv.delays); -end - -if model.rvt.include - nPhysioRegressors = nPhysioRegressors + numel(model.rvt.delays); -end - - -if model.retroicor.include - order = model.retroicor.order; - nPhysioRegressors = nPhysioRegressors + ... - 2*order.c + ... - 2*order.r + ... - 4* order.cr; -end - -if model.noise_rois.include - % TODO: what if number of components implicit?...shall we save this? - nPhysioRegressors = nPhysioRegressors + ... - numel(model.noise_rois.roi_files)*ceil(model.noise_rois.n_components+1); % + 1 for mean -end - -if model.movement.include - % TODO: what about variable regressors, that should not be - % concatenated, e.g. movement outlier censoring - nPhysioRegressors = nPhysioRegressors + model.movement.order; -end - -end diff --git a/PhysIO/code/tapas_physio_create_LOCS_from_VOLLOCS.m b/PhysIO/code/tapas_physio_create_LOCS_from_VOLLOCS.m deleted file mode 100644 index 531eee08..00000000 --- a/PhysIO/code/tapas_physio_create_LOCS_from_VOLLOCS.m +++ /dev/null @@ -1,82 +0,0 @@ -function LOCS = tapas_physio_create_LOCS_from_VOLLOCS(VOLLOCS, t, sqpar) -% Compute slice scan events as locations in time vector (LOCS) from volume -% scan locations and sequence parameters (TR, nSlices) -% -% LOCS = tapas_physio_create_LOCS_from_VOLLOCS(VOLLOCS, t, sqpar); -% -% -% IN -% NOTE: The detailed description of all input structures can be found as -% comments in tapas_physio_new -% -% VOLLOCS - index locations in time vector (of physiological recordings), -% when volume scan events started -% t - time vector of phys time course -% -% sqpar - sequence timing parameters, used for computation -% of scan events from 'nominal' timing -% .Nslices - number of slices per volume in fMRI scan -% .TR - repetition time in seconds -% .time_slice_to_slice - time between the acquisition of 2 subsequent -% slices; typically TR/Nslices or -% minTR/Nslices, if minimal temporal slice -% spacing was chosen -% -% OUT -% LOCS - locations in time vector, when slice or volume scan -% events started -% -% See also tapas_physio_new tapas_physio_main_create_regresssors -% -% Author: Lars Kasper -% Created: 2013-08-23 -% Copyright (C) 2016 Institute for Biomedical Engineering, ETH/Uni Zurich. -% -% This file is part of the PhysIO toolbox, which is released under the terms of the GNU General Public -% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL -% (either version 3 or, at your option, any later version). For further details, see the file -% COPYING or . -% -% $Id$ - -NallVols = sqpar.Ndummies+sqpar.Nscans; -Nslices = sqpar.Nslices; -TR = sqpar.TR; -LOCS = NaN(NallVols*Nslices,1); - -%default for equidistantly spaced slices -if isempty(sqpar.time_slice_to_slice) - sqpar.time_slice_to_slice = TR/Nslices; -end - -iFirstValidVolume = find(~isnan(VOLLOCS), 1); -tRef = t(VOLLOCS(iFirstValidVolume)); - - -for n = iFirstValidVolume:NallVols - - % first, restrict search interval to within TR - LOCVOLSTART = VOLLOCS(n); - - if n == NallVols - LOCVOLEND = numel(t); - else - LOCVOLEND = VOLLOCS(n+1); - end - - % Then, find slice-by-slice nearest neighbour of actual and - % reference time - tSearchInterval = t(LOCVOLSTART:LOCVOLEND) - ... - (tRef + sqpar.TR*(n-iFirstValidVolume)); - - % Note: Why not the same as other slices? Some jitter problem... - % if TR is not given with high enough precision, therefore reset for - % each volume! - LOCS((n-1)*Nslices + 1) = LOCVOLSTART; - for s = 2:sqpar.Nslices - [tmp, RELATIVELOC] = min(abs(tSearchInterval - ... - sqpar.time_slice_to_slice*(s-1)) ); - - LOCS((n-1)*Nslices + s) = LOCVOLSTART + RELATIVELOC - 1; - end -end \ No newline at end of file diff --git a/PhysIO/code/tapas_physio_create_hrv_regressors.m b/PhysIO/code/tapas_physio_create_hrv_regressors.m deleted file mode 100644 index 4efb2b8c..00000000 --- a/PhysIO/code/tapas_physio_create_hrv_regressors.m +++ /dev/null @@ -1,136 +0,0 @@ -function [convHRVOut, hrOut, verbose] = tapas_physio_create_hrv_regressors(... - ons_secs, sqpar, model_hrv, verbose) -% computes cardiac response function regressor and heart rate -% -% [convHRV, hr] = tapas_physio_create_hrv_regressors(ons_secs, sqpar ) -% -% Reference: -% Chang, Catie, John P. Cunningham, and Gary H. Glover. -% Influence of Heart Rate on the BOLD Signal: The Cardiac Response Function. -% NeuroImage 44, no. 3 (February 1, 2009): 857-869. -% doi:10.1016/j.neuroimage.2008.09.029. -% -% IN -% ons_secs. -% cpulse onset times (seconds) of heartbeat pulses (R-wave peak) -% spulse_per_vol See also tapas_physio_get_sample_points -% sqpar. -% onset_slice -% -% OUT -% convHRV cardiac response function regressor after convolution . See -% also -% EXAMPLE -% [convHRV, hr] = tapas_physio_create_hrv_regressors(physio_out.ons_secs, ... -% physio_out.sqpar); -% -% See also tapas_physio_hr tapas_physio_crf -% -% Author: Lars Kasper -% Created: 2013-07-26 -% Copyright (C) 2013 TNU, Institute for Biomedical Engineering, University of Zurich and ETH Zurich. -% -% This file is part of the physIO toolbox, which is released under the terms of the GNU General Public -% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL -% (either version 3 or, at your option, any later version). For further details, see the file -% COPYING or . -% -% $Id$ -if nargin < 3 - physio = tapas_physio_new; - model_hrv = physio.model.hrv; -end - -delays = model_hrv.delays; - -if nargin < 4 - verbose.level = []; - verbose.fig_handles = []; -end - -slicenum = 1:sqpar.Nslices; - -sample_points = tapas_physio_get_sample_points(ons_secs, sqpar, slicenum); -hr = tapas_physio_hr(ons_secs.cpulse, sample_points); - -if verbose.level>=2 - verbose.fig_handles(end+1) = tapas_physio_get_default_fig_params(); - set(gcf, 'Name', 'Regressors Heart Rate: HRV X CRF'); - subplot(2,2,1) - plot(sample_points,hr,'r');xlabel('time (seconds)'); - title('Heart Rate'); - ylabel('beats per min (bpm)'); -end - -% create convolution for whole time series first... -dt = sqpar.TR/sqpar.Nslices; -t = 0:dt:32; % 32 seconds regressor -crf = tapas_physio_crf(t); -crf = crf/max(abs(crf)); -% crf = spm_hrf(dt); -if verbose.level>=2 - subplot(2,2,2) - plot(t, crf,'r');xlabel('time (seconds)'); - title('Cardiac response function'); -end - -% NOTE: the removal of the mean was implemented to avoid over/undershoots -% at the 1st and last scans of the session due to convolution -convHRV = conv(hr-mean(hr), crf, 'same'); - -% rescaling to -1...1 for display purposes -convHRV = convHRV./max(abs(convHRV)); - -if verbose.level>=2 - subplot(2,2,3) - plot(sample_points, convHRV,'r');xlabel('time (seconds)'); - title('Heart rate X cardiac response function'); -end - - -% create shifted regressors convolved time series, which is equivalent to -% delayed response functions according to Wikipedia (convoution) -% -% "Translation invariance[edit] -% The convolution commutes with translations, meaning that -% -% \tau_x ({f}*g) = (\tau_x f)*g = {f}*(\tau_x g)\, -% where \tau_x fis the translation of the function f by x defined by -% (\tau_x f)(y) = f(y-x).\. - -% remove mean and linear trend to fulfill periodicity condition for -% shifting -convHRV = detrend(convHRV); - - -% TODO: what happens at the end/beginning of shifted convolutions? -nDelays = numel(delays); -nShiftSamples = ceil(delays/dt); - -% resample to slices needed -nSampleSlices = numel(sqpar.onset_slice); -nScans = numel(sample_points(sqpar.onset_slice:sqpar.Nslices:end)); - -hrOut = zeros(nScans,nSampleSlices); -convHRVOut = zeros(nScans,nDelays,nSampleSlices); -samplePointsOut = zeros(nScans,nSampleSlices); - -for iDelay = 1:nDelays - convHRVShifted = circshift(convHRV, nShiftSamples(iDelay)); - for iSlice = 1:nSampleSlices - onset_slice = sqpar.onset_slice(iSlice); - hrOut(:,iSlice) = hr(onset_slice:sqpar.Nslices:end)'; - convHRVOut(:,iDelay,iSlice) = convHRVShifted(onset_slice:sqpar.Nslices:end); - samplePointsOut(:,iSlice) = sample_points(onset_slice:sqpar.Nslices:end); - end -end - -if verbose.level>=2 - subplot(2,2,4) - [tmp, iShiftMin] = min(abs(delays)); - - hp{1} = plot(samplePointsOut, hrOut,'k--'); hold all; - hp{2} = plot(samplePointsOut, squeeze(convHRVOut(:,iShiftMin,:)),'r'); - xlabel('time (seconds)');ylabel('regessor'); - legend([hp{1}(1), hp{2}(1)], 'heart rate (bpm)', 'cardiac response regressor'); -end diff --git a/PhysIO/code/tapas_physio_create_missing_physio_contrasts.m b/PhysIO/code/tapas_physio_create_missing_physio_contrasts.m deleted file mode 100644 index 71e4609c..00000000 --- a/PhysIO/code/tapas_physio_create_missing_physio_contrasts.m +++ /dev/null @@ -1,79 +0,0 @@ -function [SPM, matlabbatch, indContrastsExisting, indContrastsCreate] = ... - tapas_physio_create_missing_physio_contrasts(SPM, model, namesPhysContrasts) -% Creates all valid physiological subset regressors for given PhysIO model -% that are not already existing in SPM, ignores -% -% [matlabbatch, indContrastsExisting, indContrastsCreate] = ... -% tapas_physio_create_missing_physio_contrasts(SPM, model, namesPhysContrasts) -% -% IN -% SPM SPM structured variable -% model physio.model -% namesPhysContrasts -% optional cell of strings determining which contrasts shall be -% listed for checking and, if needed, creation -% NOTE: invalid names will be ignored -% default: all phys contrasts will be queried for -% creation, as listed in -% -% OUT -% SPM SPM structured variable with updated contrasts -% matlabbatch -% that was used to create new valid contrasts -% indContrastsExisting -% indices of namesPhysContrasts that existed in design matrix -% indContrastsCreate -% indices of namesPhysContrasts that had to be created -% -% EXAMPLE -% tapas_physio_create_missing_physio_contrasts -% -% See also tapas_physio_check_prepare_job_contrasts -% See also tapas_physio_report_contrasts -% -% Author: Lars Kasper -% Created: 2016-10-03 -% Copyright (C) 2016 TNU, Institute for Biomedical Engineering, -% University of Zurich and ETH Zurich. -% -% This file is part of the TAPAS PhysIO Toolbox, which is released under the terms of the GNU General Public -% License (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL -% (either version 3 or, at your option, any later version). For further details, see the file -% COPYING or . -% -% $Id$ - -if nargin < 3 - namesPhysContrasts = tapas_physio_get_contrast_names_default(); -end - -nContrasts = numel(namesPhysContrasts); -fileSpm = fullfile(SPM.swd, 'SPM.mat'); - -%% Check whether contrasts already exist in SPM.mat -indContrasts = zeros(nContrasts,1); -for c = 1:nContrasts - indContrasts(c) = tapas_physio_check_get_xcon_index(SPM, ... - namesPhysContrasts{c}); -end - - -%% Generate contrasts only if not already existing - -indContrastsExisting = find(indContrasts); - -if ~isempty(model) - indContrastsCreate = find(~indContrasts); - namesContrastsCreate = namesPhysContrasts(indContrastsCreate); - matlabbatch = tapas_physio_check_prepare_job_contrasts(fileSpm, ... - model, SPM, indContrastsCreate, namesContrastsCreate); - if ~isempty(matlabbatch{1}.spm.stats.con.consess) - spm_jobman('run', matlabbatch); - end -else - error('No physio.model specified'); -end - -if nargout - load(fileSpm); -end \ No newline at end of file diff --git a/PhysIO/code/tapas_physio_create_movement_regressors.m b/PhysIO/code/tapas_physio_create_movement_regressors.m deleted file mode 100644 index bdd61244..00000000 --- a/PhysIO/code/tapas_physio_create_movement_regressors.m +++ /dev/null @@ -1,117 +0,0 @@ -function [R, verbose] = tapas_physio_create_movement_regressors(movement, verbose) -% Reads realignment parameters, creates derivative/squared & outlier regressor -% -% [R, verbose] = tapas_physio_create_movement_regressors(movement, verbose) -% -% The 6 realignment parameters can be augmented by their derivatives (in -% total 12 parameters) + the squares of parameters and derivatives (in -% total 24 parameters), leading to a Volterra expansion as in -% -% Friston KJ, Williams S, Howard R, Frackowiak -% RS, Turner R. Movement-related effects in fMRI -% time-series. Magn Reson Med. 1996;35:346?355.) -% -% IN -% movement physio.model.movement -% verbose physio.verbose -% OUT -% R [nScans, (6|12|24)+nOutliers] regressor matrix from movement -% -% EXAMPLE -% tapas_physio_create_movement_regressors -% -% See also tapas_physio_new tapas_physio_main_create_regressors -% -% Author: Lars Kasper -% Created: 2015-07-10 -% Copyright (C) 2015 TNU, Institute for Biomedical Engineering, -% University of Zurich and ETH Zurich. -% -% This file is part of the TAPAS PhysIO Toolbox, which is released under the terms of the GNU General Public -% License (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL -% (either version 3 or, at your option, any later version). For further details, see the file -% COPYING or . -% -% $Id$ - - -[fp, fn, fs] = fileparts(movement.file_realignment_parameters); - -if ~exist(movement.file_realignment_parameters, 'file') - verbose = tapas_physio_log('No input multiple regressors found', verbose, 1); - R = []; -else - rp = load(movement.file_realignment_parameters); - if strcmp('.txt', fs) % text-file - R = rp; - else % mat file - R = rp.R; - end - - % Include derivatives (needed later for outliers as well - - dR= diff(R); - dR = [zeros(1,size(R,2)); dR]; - - if movement.order > 6 - R = [R, dR]; - end - - - % Include squared regressors/derivatives - - if movement.order > 12 - R = [R , R.^2]; - end - - - % Include outlier movements exceeding thresholds as stick regressors - % Euclidean distance used! - - sumTrans = sqrt(sum(dR(:,1:3).^2,2)); - sumRot = sqrt(sum(dR(:,4:6).^2,2)); - iOutlierTrans = find(sumTrans > movement.outlier_translation_mm); - iOutlierRot = find( sumRot > movement.outlier_rotation_deg/180*pi); - - nOutlierTrans = numel(iOutlierTrans); - nOutlierRot = numel(iOutlierRot); - if verbose.level > 2 - verbose.fig_handles(end+1) = tapas_physio_get_default_fig_params; - subplot(2,1,1); - plot(R(:,1:3),'-'); hold on - plot(sumTrans, '-'); - if nOutlierTrans - stem( iOutlierTrans,... - movement.outlier_translation_mm*ones(1,nOutlierTrans)) - end - legend('shift x (mm)', 'shift y (mm)', 'shift z (mm)', 'RMS diff translation (mm)', ... - 'excess translation volumes'); - xlabel('scans'); - ylabel('translation (mm)'); - - subplot(2,1,2); - plot(R(:,4:6)*180/pi,'-'); hold on - plot(sumRot*180/pi, '-'); - if nOutlierRot - stem(iOutlierRot,... - movement.outlier_rotation_deg*ones(1,nOutlierRot)) - end - legend('pitch x (deg)', 'roll y (deg)', 'yaw z (deg)', 'RMS diff rotation (deg)', ... - 'excess rotation volumes'); - - xlabel('scans'); - ylabel('rotation (deg)'); - end - - iOutlierArray = unique([iOutlierTrans; iOutlierRot]); - nOutliers = numel(iOutlierArray); - - nScans = size(R,1); - R_outlier = zeros(nScans, nOutliers); - for iOutlier = 1:nOutliers - R_outlier(iOutlierArray(iOutlier), iOutlier) = 1; - end - - R = [R, R_outlier]; - -end diff --git a/PhysIO/code/tapas_physio_create_noise_rois_regressors.m b/PhysIO/code/tapas_physio_create_noise_rois_regressors.m deleted file mode 100644 index daf53660..00000000 --- a/PhysIO/code/tapas_physio_create_noise_rois_regressors.m +++ /dev/null @@ -1,213 +0,0 @@ -function [R_noise_rois, noise_rois, verbose] = tapas_physio_create_noise_rois_regressors(... - noise_rois, verbose) -% Compute physiological regressors by extracting principal components of -% the preprocessed fMRI time series from anatomically defined noise ROIS (e.g. CSF) -% -% [R_noise_rois, verbose] = tapas_physio_create_noise_rois_regressors(... -% noise_rois, verbose) -% -% NOTE: The mean of all time series in each ROI will also be added as a regressor -% automatically -% -% Approach similar to the one described as aCompCor: -% Behzadi, Y., Restom, K., Liau, J., Liu, T.T., 2007. -% A component based noise correction method (CompCor) for BOLD and -% perfusion based fMRI. NeuroImage 37, 90?101. -% doi:10.1016/j.neuroimage.2007.04.042 -% -% IN -% physio.model.noise_rois -% OUT -% -% EXAMPLE -% tapas_physio_create_noise_rois_regressors -% -% See also spm_ov_roi -% -% Author: Lars Kasper -% Created: 2015-07-22 -% Copyright (C) 2015 TNU, Institute for Biomedical Engineering, -% University of Zurich and ETH Zurich. -% -% This file is part of the TAPAS PhysIO Toolbox, which is released under the terms of the GNU General Public -% License (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL -% (either version 3 or, at your option, any later version). For further details, see the file -% COPYING or . -% -% $Id$ - -% TODO: visualization of Rois, components/mean and spatial loads in ROIs - -tapas_physio_strip_fields(noise_rois); - -if isempty(fmri_files) || isempty(roi_files) - R_noise_rois = []; - return -end - -if ~iscell(fmri_files) - fmri_files = cellstr(fmri_files); -end - -if ~iscell(roi_files) - roi_files = cellstr(roi_files); -end - -nRois = numel(roi_files); - -if numel(thresholds) == 1 - thresholds = repmat(thresholds, 1, nRois); -end - -if numel(n_voxel_crop) == 1 - n_voxel_crop = repmat(n_voxel_crop, 1, nRois); -end - -if numel(n_components) == 1 - n_components = repmat(n_components, 1, nRois); -end - -% TODO: what if different geometry of mask and fmri data? -% or several fmri files given? -Vimg = spm_vol(fmri_files{1}); -Yimg = spm_read_vols(Vimg); - -nVolumes = size(Yimg, 4); -Yimg = reshape(Yimg, [], nVolumes); - -R_noise_rois = []; -for r = 1:nRois - - Vroi = spm_vol(roi_files{r}); - - hasRoiDifferentGeometry = any(any(abs(Vroi.mat - Vimg(1).mat) > 1e-5)) | ... - any(Vroi.dim-Vimg(1).dim(1:3)); - hasRoiDifferentGeometry = true; - - if hasRoiDifferentGeometry - - % reslice to same geometry - matlabbatch{1}.spm.spatial.coreg.write.ref = fmri_files(1); - matlabbatch{1}.spm.spatial.coreg.write.source = roi_files(r); - matlabbatch{1}.spm.spatial.coreg.write.roptions.interp = 4; - matlabbatch{1}.spm.spatial.coreg.write.roptions.wrap = [0 0 0]; - matlabbatch{1}.spm.spatial.coreg.write.roptions.mask = 0; - matlabbatch{1}.spm.spatial.coreg.write.roptions.prefix = 'r'; - - spm_jobman('run', matlabbatch); - - % update header link to new reslice mask-file - Vroi = spm_vol(spm_file(roi_files{r}, 'prefix', 'r')); - end - - - roi = spm_read_vols(Vroi); - roi(roi < thresholds(r)) = 0; - roi(roi >= thresholds(r)) = 1; - - % crop pixel, if desired - if n_voxel_crop(r) - nSlices = size(roi,3); - for s = 1:nSlices - roi(:,:,s) = imerode(roi(:,:,s), strel('disk', n_voxel_crop(r))); - end - end - - Yroi = Yimg(roi(:)==1, :); - - %% mean and linear trend removal according to CompCor pub - % design matrix - X = ones(nVolumes,1); - X(:,2) = 1:nVolumes; - % fit 1st order polynomial to time series data in each voxel - for n_roi_voxel = 1:size(Yroi,1) - % extract data - raw_Y = Yroi(n_roi_voxel,:)'; - % estimate betas - beta = X\raw_Y; - % fitted data - fit_Y = X*beta; - % detrend data - detrend_Y = raw_Y - fit_Y; - - % overwrite Yroi - Yroi(n_roi_voxel,:) = detrend_Y; - end - - - - - - %% - - - % COEFF = [nVolumes, nPCs] principal components (PCs) ordered by variance - % explained - % SCORE = [nVoxel, nPCs] loads of each component in each voxel, i.e. - % specific contribution of each component in - % a voxel's variance - % LATENT = [nPCs, 1] eigenvalues of data covariance matrix, - % stating how much variance was explained by - % each PC overall - % TSQUARED = [nVoxels,1] Hotelling's T-Squared test whether PC - % explained significant variance in a voxel - % EXPLAINED = [nPCs, 1] relative amount of variance explained (in - % percent) by each component - % MU = [1, nVolumes] mean of all time series - [COEFF, SCORE, LATENT, TSQUARED, EXPLAINED, MU] = pca(Yroi); - - % components either via number or threshold of variance explained - if n_components(r) >= 1 - nComponents = n_components(r); - else - nComponents = find(cumsum(EXPLAINED)/100 > n_components(r), 1, 'first'); - end - - % save to return - noise_rois.n_components(r) = nComponents + 1; % + 1 for mean - - % Take mean and some components into noise regressor - R = MU'; - R = [R, COEFF(:,1:nComponents)]; - - nRegressors = size(R,2); - - % z-transform - stringLegend = cell(nRegressors,1); - for c = 1:nRegressors - R(:,c) = (R(:,c) - mean(R(:,c)))/std(R(:,c)); - - if c > 1 - stringLegend{c} = ... - sprintf('Principal component %d (%7.4f %% variance explained)', ... - c-1, EXPLAINED(c-1)); - else - stringLegend{c} = 'Mean time series of all voxels'; - end - end - - if verbose.level >=2 - stringFig = sprintf('Noise_rois: Extracted principal components for ROI %d', r); - verbose.fig_handles(end+1) = tapas_physio_get_default_fig_params(); - set(gcf, 'Name', stringFig); - plot(R); - legend(stringLegend); - title(stringFig); - xlabel('scans'); - ylabel('PCs, mean-centred and std-scaled'); - end - - % write away extracted PC-loads & roi of extraction - [tmp,fnRoi] = fileparts(Vroi(1).fname); - fpFmri = fileparts(Vimg(1).fname); - for c = 1:nComponents - Vpc = Vroi; - Vpc.fname = fullfile(fpFmri, sprintf('pc%02d_scores_%s.nii',c, fnRoi)); - % saved as float, since was masked before - Vpc.dt = [spm_type('float32') 1]; - pcScores = zeros(Vpc.dim); - pcScores(roi(:)==1) = SCORE(:, c); - spm_write_vol(Vpc, pcScores); - end - R_noise_rois = [R_noise_rois, R]; -end \ No newline at end of file diff --git a/PhysIO/code/tapas_physio_create_retroicor_regressors.m b/PhysIO/code/tapas_physio_create_retroicor_regressors.m deleted file mode 100644 index 998e50ef..00000000 --- a/PhysIO/code/tapas_physio_create_retroicor_regressors.m +++ /dev/null @@ -1,173 +0,0 @@ -function [cardiac_sess, respire_sess, mult_sess, ons_secs, order, verbose, ... - c_sample_phase, r_sample_phase] ... - = tapas_physio_create_retroicor_regressors(ons_secs, sqpar, order, verbose) -% calculation of regressors for physiological motion correction using RETROICOR (Glover, MRM, 2000) -% -% USAGE: -% [cardiac_sess, respire_sess, mult_sess, verbose, c_sample_phase, r_sample_phase] ... -% = tapas_physio_create_retroicor_regressors(ons_secs, sqpar, thresh, slicenum, order, verbose) -% -% NOTE: also updates order of models to 0, if some data does not exist -% (cardiac or respiratory) -% -% INPUT: -% ons_secs - onsets of all physlog events in seconds -% .spulse = onsets of slice scan acquisition -% .cpulse = onsets of cardiac R-wave peaks -% .r = time series of respiration -% .svolpulse = onsets of volume scan acquisition -% .t = time vector of logfile rows -% -% sqpar - sequence timing parameters -% .Nslices - number of slices per volume in fMRI scan -% .NslicesPerBeat - usually equals Nslices, unless you trigger -% with the heart beat -% .TR - repetition time in seconds -% .Ndummies - number of dummy volumes -% .Nscans - number of full volumes saved (volumes in nifti file, -% usually rows in your design matrix) -% onset_slice - slice whose scan onset determines the adjustment of the -% regressor timing to a particular slice for the whole volume -% -% -% ------------------------------------------------------------------------- -% Lars Kasper, March 2012 -% Copyright (C) 2013 Institute for Biomedical Engineering, ETH/Uni Zurich. -% -% This file is part of the PhysIO toolbox, which is released under the terms of the GNU General Public -% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL -% (either version 3 or, at your option, any later version). For further details, see the file -% COPYING or . -% -% $Id$ - -%% variable renaming -if ~exist('verbose', 'var') - verbose.level = 1; -end - -svolpulse = ons_secs.svolpulse; -cpulse = ons_secs.cpulse; -r = ons_secs.r; -spulse = ons_secs.spulse; -t = ons_secs.t; - -hasRespData = ~isempty(r); -hasPhaseData = isfield(ons_secs, 'c_sample_phase') && ~isempty(ons_secs.c_sample_phase); -hasCardiacData = hasPhaseData || ~isempty(cpulse); - -if ~hasPhaseData - - % update model order, if resp/cardiac data do not exist - if ~hasCardiacData - order.c = 0; - end - - if ~hasRespData - order.r = 0; - end - - if ~hasRespData || ~hasCardiacData - order.cr = 0; - end - - % compute phases from pulse data - - % compute differently, i.e. separate regressors for multiple slice - % generation - nSampleSlices = numel(sqpar.onset_slice); - hasMultipleSlices = nSampleSlices>1; - - %parameters for resampling - rsampint = t(2)-t(1); - - %% Get phase, downsample and Fourier-expand - sample_points = tapas_physio_get_sample_points(ons_secs, sqpar); - - % cardiac phase estimation and Fourier expansion - if (order.c || order.cr) - - [c_phase, verbose] = ... - tapas_physio_get_cardiac_phase(cpulse, spulse, verbose, svolpulse); - c_sample_phase = tapas_physio_downsample_phase(spulse, c_phase, sample_points, rsampint); - cardiac_sess = tapas_physio_get_fourier_expansion(c_sample_phase,order.c); - - cardiac_sess = tapas_physio_split_regressor_slices(cardiac_sess, ... - nSampleSlices); - - - - else - cardiac_sess = []; - c_sample_phase = []; - end - - % Respiratory phase estimation and Fourier expansion - if (order.r || order.cr) - - fr = ons_secs.fr; - - if verbose.level >=3 - [r_phase, verbose.fig_handles(end+1)] = ... - tapas_physio_get_respiratory_phase( ... - fr,rsampint, verbose.level); - else - r_phase = tapas_physio_get_respiratory_phase(fr,rsampint, 0); - end - r_sample_phase = tapas_physio_downsample_phase(t, r_phase, sample_points, rsampint); - - respire_sess = tapas_physio_get_fourier_expansion(r_sample_phase,order.r); - respire_sess = tapas_physio_split_regressor_slices(respire_sess, ... - nSampleSlices); - - else - respire_sess = []; - r_sample_phase =[]; - end - -else % compute Fourier expansion directly from cardiac/respiratory phases - % select subset of slice-wise sampled phase for current onset_slice - c_sample_phase = ons_secs.c_sample_phase; - r_sample_phase = ons_secs.r_sample_phase; - if (order.c || order.cr) - cardiac_sess = tapas_physio_get_fourier_expansion(... - c_sample_phase, order.c); - respire_sess = tapas_physio_get_fourier_expansion(... - r_sample_phase, order.r); - else - cardiac_sess = []; - respire_sess = []; - end - -end - - -% Multiplicative terms as specified in Harvey et al., 2008 -if order.cr && hasRespData && hasCardiacData - crplus_sess = tapas_physio_get_fourier_expansion(c_sample_phase+r_sample_phase,order.cr); - crdiff_sess = tapas_physio_get_fourier_expansion(c_sample_phase-r_sample_phase,order.cr); - mult_sess = [crplus_sess crdiff_sess]; - mult_sess = tapas_physio_split_regressor_slices(mult_sess, ... - nSampleSlices); -else - mult_sess = []; -end - -ons_secs.c_sample_phase = tapas_physio_split_regressor_slices(... - c_sample_phase, nSampleSlices); -ons_secs.r_sample_phase = tapas_physio_split_regressor_slices(... - r_sample_phase, nSampleSlices); - - -%% plot cardiac & resp. regressors -if verbose.level >=2 - R = [cardiac_sess, respire_sess, mult_sess]; - - hasCardiacData = ~isempty(ons_secs.c); - hasRespData = ~isempty(ons_secs.r); - - verbose.fig_handles(end+1) = ... - tapas_physio_plot_retroicor_regressors(R, order, hasCardiacData, ... - hasRespData); - -end diff --git a/PhysIO/code/tapas_physio_create_rvt_regressors.m b/PhysIO/code/tapas_physio_create_rvt_regressors.m deleted file mode 100644 index fddfceac..00000000 --- a/PhysIO/code/tapas_physio_create_rvt_regressors.m +++ /dev/null @@ -1,140 +0,0 @@ -function [convRVTOut, rvtOut, verbose] = tapas_physio_create_rvt_regressors(... - ons_secs, sqpar, model_rvt, verbose) -% computes respiratory response function regressor and respiratory volume per time -% -% [convHRV, hr] = tapas_physio_create_rvt_regressors(ons_secs, sqpar ) -% -% Reference: -% Birn, R.M., Smith, M.A., Jones, T.B., Bandettini, P.A., 2008. -% The respiration response function: The temporal dynamics of -% fMRI signal fluctuations related to changes in respiration. -% NeuroImage 40, 644-654. -% -% IN -% ons_secs. -% fr filtered respiratory signal time series -% spulse_per_vol See also tapas_physio_get_sample_points -% sqpar. -% onset_slice -% OUT -% convRVTOut [nScans, nDelays, nSampleSlices] -% respiratory response function regressor after -% convolution for specified delays and downsampled -% to given slices. -% EXAMPLE -% [convHRV, hr] = tapas_physio_create_hrv_regressor(physio_out.ons_secs, physio_out.sqpar); -% -% See also tapas_physio_rvt tapas_physio_rrf -% -% Author: Lars Kasper -% Created: 2014-01-20 -% Copyright (C) 2014 TNU, Institute for Biomedical Engineering, University of Zurich and ETH Zurich. -% -% This file is part of the physIO toolbox, which is released under the terms of the GNU General Public -% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL -% (either version 3 or, at your option, any later version). For further details, see the file -% COPYING or . -% -% $Id$ -if nargin < 3 - physio = tapas_physio_new; - model_rvt = physio.model.rvt; -end - -delays = model_rvt.delays; - - -if nargin < 4 - verbose.level = []; - verbose.fig_handles = []; -end - -slicenum = 1:sqpar.Nslices; - -sample_points = tapas_physio_get_sample_points(ons_secs, sqpar, slicenum); -rvt = tapas_physio_rvt(ons_secs.fr, ons_secs.t, sample_points, verbose); -rvt = rvt/max(rvt); % normalize for reasonable range of regressor - -if verbose.level >=2 - verbose.fig_handles(end+1) = tapas_physio_get_default_fig_params(); - set(gcf, 'Name', 'Convolution Respiration RVT X RRF'); - subplot(2,2,1) - plot(sample_points,rvt, 'g');xlabel('time (seconds)'); - title('Respiratory volume per time'); - ylabel('a.u.'); -end - -% create convolution for whole time series first... -dt = sqpar.TR/sqpar.Nslices; -t = 0:dt:50; % 50 seconds regressor -rrf = tapas_physio_rrf(t); -rrf = rrf/max(abs(rrf)); -% crf = spm_hrf(dt); - -if verbose.level >= 2 - subplot(2,2,2) - plot(t, rrf,'g');xlabel('time (seconds)'); - title('Respiratory response function'); -end - -% NOTE: the removal of the mean was implemented to avoid over/undershoots -% at the 1st and last scans of the session due to convolution -convRVT = conv(rvt-mean(rvt), rrf, 'same'); - -% rescaling to -1...1 for display purposes -convRVT = convRVT./max(abs(convRVT)); - - -if verbose.level >= 2 - subplot(2,2,3) - plot(sample_points, convRVT,'g');xlabel('time (seconds)'); - title('Resp vol time X resp response function'); -end - -% create shifted regressors convolved time series, which is equivalent to -% delayed response functions according to Wikipedia (convoution) -% -% "Translation invariance[edit] -% The convolution commutes with translations, meaning that -% -% \tau_x ({f}*g) = (\tau_x f)*g = {f}*(\tau_x g)\, -% where \tau_x fis the translation of the function f by x defined by -% (\tau_x f)(y) = f(y-x).\. - -% remove mean and linear trend to fulfill periodicity condition for -% shifting -convRVT = detrend(convRVT); - - -% TODO: what happens at the end/beginning of shifted convolutions? -nDelays = numel(delays); -nShiftSamples = ceil(delays/dt); - -% resample to slices needed -nSampleSlices = numel(sqpar.onset_slice); -nScans = numel(sample_points(sqpar.onset_slice:sqpar.Nslices:end)); - -rvtOut = zeros(nScans,nSampleSlices); -convRVTOut = zeros(nScans,nDelays,nSampleSlices); -samplePointsOut = zeros(nScans,nSampleSlices); - -for iDelay = 1:nDelays - convRVTShifted = circshift(convRVT, nShiftSamples(iDelay)); - for iSlice = 1:nSampleSlices - onset_slice = sqpar.onset_slice(iSlice); - rvtOut(:,iSlice) = rvt(onset_slice:sqpar.Nslices:end)'; - convRVTOut(:,iDelay,iSlice) = convRVTShifted(onset_slice:sqpar.Nslices:end); - samplePointsOut(:,iSlice) = sample_points(onset_slice:sqpar.Nslices:end); - end -end - -if verbose.level >= 2 - subplot(2,2,4) - [tmp, iShiftMin] = min(abs(delays)); - hp{1} = plot(samplePointsOut, rvtOut,'k--');hold all; - hp{2} = plot(samplePointsOut, squeeze(convRVTOut(:,iShiftMin,:)),'g'); - xlabel('time (seconds)'); - title('RVT regessor'); - legend([hp{1}(1), hp{2}(1)], 'respiratory volume / time (a. u.)', ... - 'respiratory response regressor'); -end diff --git a/PhysIO/code/tapas_physio_create_scan_timing.m b/PhysIO/code/tapas_physio_create_scan_timing.m deleted file mode 100644 index f459aa80..00000000 --- a/PhysIO/code/tapas_physio_create_scan_timing.m +++ /dev/null @@ -1,115 +0,0 @@ -function [ons_secs, VOLLOCS, LOCS, verbose] = tapas_physio_create_scan_timing(... - log_files, scan_timing, ons_secs, verbose) -% Extracts slice and volume scan onsets (triggers) from different vendor formats -% -% [ons_secs, VOLLOCS, LOCS, verbose] = tapas_physio_create_scan_timing(... -% log_files, scan_timing, ons_secs, verbose); -% -% -% IN -% NOTE: The detailed description of all input structures can be found as -% comments in tapas_physio_new -% -% log_files - file names (physiology and scan timing) and sampling rates -% -% ons_secs - structure for time-dependent variables, i.e. onsets, -% specified in seconds, in particular -% .t - time vector of phys time course -% -% scan_timing - Parameters for sequence timing & synchronization -% scan_tming.sqpar - slice and volume acquisition starts, TR, -% number of scans etc. -% scan_timing.sync - synchronization options -% (e.g. from gradients, trigger, tics phys -% logfile to scan acquisition) -% -% sqpar - sequence timing parameters, used for computation -% of scan events from 'nominal' timing -% .Nslices - number of slices per volume in fMRI scan -% .TR - repetition time in seconds -% .Ndummies - number of dummy volumes -% .Nscans - number of full volumes saved (volumes in nifti file, -% usually rows in your design matrix) -% .Nprep - number of non-dummy, volume like preparation pulses -% before 1st dummy scan. If set, logfile is read from beginning, -% otherwise volumes are counted from last detected volume in the logfile -% .time_slice_to_slice - time between the acquisition of 2 subsequent -% slices; typically TR/Nslices or -% minTR/Nslices, if minimal temporal slice -% spacing was chosen -% -% verbose - defines output level (which graphics to plot -% and whether to save them) -% -% OUT -% ons_secs - structure for time-dependent variables, i.e. onsets, -% specified in seconds, updated fields -% .spulse - scan slice trigger events -% .svolpulse - scan volume trigger events -% .spulse_per_vol -% - cell(nVolumes,1) of slice triggers per -% volume -% .acq_codes - acquisition codes (e.g. triggers) within -% phys log files (e.g. Philips, Biopac) -% -% VOLLOCS - index locations in time vector (of physiological recordings), -% when volume scan events started -% LOCS - locations in time vector, when slice or volume scan -% events started -% -% See also tapas_physio_new tapas_physio_main_create_regresssors -% -% Author: Lars Kasper -% Created: 2013-08-23 -% Copyright (C) 2016 Institute for Biomedical Engineering, ETH/Uni Zurich. -% -% This file is part of the PhysIO toolbox, which is released under the terms of the GNU General Public -% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL -% (either version 3 or, at your option, any later version). For further details, see the file -% COPYING or . -% -% $Id$ - -sqpar = scan_timing.sqpar; - -% TODO: introduce auto that takes time stamps from default locations -% for different vendors -switch lower(scan_timing.sync.method) - case 'nominal' - [VOLLOCS, LOCS] = ... - tapas_physio_create_scan_timing_nominal(ons_secs.t, ... - sqpar, log_files.align_scan); - case {'gradient', 'gradient_log'} - [VOLLOCS, LOCS, verbose] = ... - tapas_physio_create_scan_timing_from_gradients_philips( ... - log_files, scan_timing, verbose); - case {'gradient_auto', 'gradient_log_auto'} - [VOLLOCS, LOCS, verbose] = ... - tapas_physio_create_scan_timing_from_gradients_auto_philips( ... - log_files, scan_timing, verbose); - case 'scan_timing_log' - switch lower(log_files.vendor) - case 'siemens' - [VOLLOCS, LOCS] = ... - tapas_physio_create_scan_timing_nominal(ons_secs.t, ... - sqpar, log_files.align_scan); - case 'siemens_tics' - [VOLLOCS, LOCS, verbose] = ... - tapas_physio_create_scan_timing_from_tics_siemens( ... - ons_secs.t, log_files, verbose); - case 'biopac_mat' - [VOLLOCS, LOCS, verbose] = ... - tapas_physio_create_scan_timing_from_acq_codes( ... - ons_secs.t, ons_secs.acq_codes, sqpar, ... - log_files.align_scan, verbose); - end -end - - -% remove arbitrary offset in time vector now, since all timings have now -% been aligned to ons_secs.t -% ons_secs.t = ons_secs.t - ons_secs.t(1); - -[ons_secs.svolpulse, ons_secs.spulse, ons_secs.spulse_per_vol, verbose] = ... - tapas_physio_get_onsets_from_locs(... - ons_secs.t, VOLLOCS, LOCS, sqpar, verbose); diff --git a/PhysIO/code/tapas_physio_create_scan_timing_from_acq_codes.m b/PhysIO/code/tapas_physio_create_scan_timing_from_acq_codes.m deleted file mode 100644 index 9bb063cf..00000000 --- a/PhysIO/code/tapas_physio_create_scan_timing_from_acq_codes.m +++ /dev/null @@ -1,104 +0,0 @@ -function [VOLLOCS, LOCS, verbose] = ... - tapas_physio_create_scan_timing_from_acq_codes(t, acq_codes, sqpar, align_scan, verbose) -% Creates scan volume and slice onsets location in time series vector from acq_codes (Philips, Biopac), -% i.e. trigger signals on same time scale as physiological recordings -% -% [VOLLOCS, LOCS, verbose] = ... -% tapas_physio_create_scan_timing_from_acq_codes(t, acq_codes, verbose) -% IN -% t - [nSamples,1] timing vector of physiological logfiles -% acq_codes - [nSamples,1] event vector of acquisition codes (e.g. -% 1 for slice trigger, 10 for volume trigger), as -% provided by tapas_physio_read_physlogfiles* -% sqpar - sequence timing parameters -% .Nslices - number of slices per volume in fMRI scan -% .Ndummies - number of dummy volumes -% .Nscans - number of full volumes saved (volumes in nifti file, -% align_scan - 'first' or 'last', for abundant triggers read from -% acq_codes, defines whether the scan-related ones are counted from the -% start ('first') or end ('last') of the sequence -% verbose - structure defining which output figures to plot and -% storing their handles -% -% OUT -% VOLLOCS - locations in time vector, when volume scan -% events started -% LOCS - locations in time vector, when slice or volume scan -% events started -% -% See also tapas_physio_create_scan_timing tapas_physio_read_physlogfiles -% -% Author: Lars Kasper -% Created: 2016-08-20 -% Copyright (C) 2014 Institute for Biomedical Engineering, ETH/Uni Zurich. -% -% This file is part of the PhysIO toolbox, which is released under the terms of the GNU General Public -% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL -% (either version 3 or, at your option, any later version). For further details, see the file -% COPYING or . -% -% $Id$ -DEBUG = verbose.level >=3; - -nVolumes = sqpar.Nscans; -nDummies = sqpar.Ndummies; -nSlices = sqpar.Nslices; - -nTotalVolumes = nVolumes + nDummies; -nTotalSlices = nTotalVolumes*nSlices; - - -doCountFromStart = strcmpi(align_scan, 'first'); - -LOCS = find(acq_codes == 1); -VOLLOCS = find(acq_codes == 10); - -isValidVOLLOCS = numel(VOLLOCS) >= nTotalVolumes; -isValidLOCS = numel(LOCS) >= nTotalSlices; - -%% crop valid location vectors to needed number of volumes or slices*volumes -if isValidLOCS - if doCountFromStart - LOCS = LOCS(1:nTotalSlices); - else - LOCS = LOCS((end-(nTotalSlices+1)):end); - end -end - -if isValidVOLLOCS - if doCountFromStart - VOLLOCS = VOLLOCS(1:nTotalVolumes); - else - VOLLOCS = VOLLOCS((end-(nTotalVolumes+1)):end); - end -end - - -%% Something invalid? try calculating from other acq_codes or use nominal -% timing! -if ~isValidLOCS || ~isValidVOLLOCS - - if isValidLOCS - % no vol events, therefore take every nSlices-th slice event - VOLLOCS = LOCS(1:nSlices:nTotalSlices); - elseif isValidVOLLOCS - LOCS = tapas_physio_create_LOCS_from_VOLLOCS(VOLLOCS, t, sqpar); - else % resort to nominal timing - [VOLLOCS, LOCS] = tapas_physio_create_scan_timing_nominal(t, sqpar); - verbose = tapas_physio_log('No gradient timecourse was logged in the logfile. Using nominal timing from sqpar instead', ... - verbose, 1); - end - -end - - -if DEBUG - verbose.fig_handles(end+1) = tapas_physio_get_default_fig_params(); - stringTitle = 'Extracted Scan Slice/Volume Onsets from Acq Codes'; - set(gcf, 'Name', stringTitle); - stem(t(LOCS), ones(size(LOCS))); hold all; - stem(t(VOLLOCS), 1.5*ones(size(VOLLOCS))); - legend('slice start events', 'volume start events'); - xlabel('t (seconds)'); - title(stringTitle); -end \ No newline at end of file diff --git a/PhysIO/code/tapas_physio_create_scan_timing_from_end_marker_philips.m b/PhysIO/code/tapas_physio_create_scan_timing_from_end_marker_philips.m deleted file mode 100644 index faa898c8..00000000 --- a/PhysIO/code/tapas_physio_create_scan_timing_from_end_marker_philips.m +++ /dev/null @@ -1,172 +0,0 @@ -function [VOLLOCS, LOCS, verbose] = tapas_physio_create_scan_timing_from_end_marker_philips(log_files, thresh, sqpar, verbose) -%extracts slice and volume scan events from gradients timecourse of Philips -% SCANPHYSLOG file -% -% [VOLLOCS, LOCS] = tapas_physio_create_scan_timing_from_gradients_philips(logfile, -% thresh); -% -% IN -% log_files is a structure containing the following filenames (with full -% path) -% .log_cardiac contains ECG or pulse oximeter time course -% for Philips: 'SCANPHYSLOG.log'; -% can be found on scanner in G:/log/scanphyslog- -% directory, one file is created per scan, make sure to take -% the one with the time stamp corresponding to your PAR/REC -% files -% .log_respiration contains breathing belt amplitude time course -% for Philips: same as .log_cardiac -% -% thresh gradient amplitude thresholds to detect slice and volume events -% -% thresh is a structure with the following elements -% .zero - gradient values below this value are set to zero; -% should be those which are unrelated to slice acquisition start -% .slice - minimum gradient amplitude to be exceeded when a slice -% scan starts -% .vol - minimum gradient amplitude to be exceeded when a new -% volume scan starts; -% leave [], if volume events shall be determined as -% every Nslices-th scan event -% .grad_direction -% - leave empty to use nominal timing; -% if set, sequence timing is calculated from logged gradient timecourse; -% - value determines which gradient direction timecourse is used to -% identify scan volume/slice start events ('x', 'y', 'z') -% .vol_spacing -% - duration (in seconds) from last slice acq to -% first slice of next volume; -% leave [], if .vol-threshold shall be used -% -% sqpar - sequence timing parameters -% .Nslices - number of slices per volume in fMRI scan -% .NslicesPerBeat - usually equals Nslices, unless you trigger with the heart beat -% .TR - repetition time in seconds -% .Ndummies - number of dummy volumes -% .Nscans - number of full volumes saved (volumes in nifti file, -% usually rows in your design matrix) -% .Nprep - number of non-dummy, volume like preparation pulses -% before 1st dummy scan. If set, logfile is read from beginning, -% otherwise volumes are counted from last detected volume in the logfile -% .TimeSliceToSlice - time between the acquisition of 2 subsequent -% slices; typically TR/Nslices or -% minTR/Nslices, if minimal temporal slice -% spacing was chosen -% onset_slice - slice whose scan onset determines the adjustment of the -% regressor timing to a particular slice for the whole volume -% -% NOTE: only necessary, if thresh.grad_direction is empty -% verbose -% -% OUT -% -% EXAMPLE -% [VOLLOCS, LOCS] = tapas_physio_create_scan_timing_from_gradients_philips(logfile, -% thresh.scan_timing); -% -% See also -% -% Author: Lars Kasper -% Created: 2013-02-16 -% Copyright (C) 2013 Institute for Biomedical Engineering, ETH/Uni Zurich. -% -% This file is part of the PhysIO toolbox, which is released under the terms of the GNU General Public -% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL -% (either version 3 or, at your option, any later version). For further details, see the file -% COPYING or . -% -% $Id$ - -% everything stored in 1 logfile -if ~isfield(log_files, 'cardiac') || isempty(log_files.cardiac) - logfile = log_files.respiration; -else - logfile = log_files.cardiac; -end - -do_detect_vol_events_by_count = (~isfield(thresh, 'vol') || isempty(thresh.vol)) && (~isfield(thresh, 'vol_spacing') || isempty(thresh.vol_spacing)); -do_detect_vol_events_by_grad_height = ~do_detect_vol_events_by_count && (~isfield(thresh, 'vol_spacing') || isempty(thresh.vol_spacing)); - -% check consistency of thresh-values - -if thresh.slice <= thresh.zero - error('Please set thresh.scan_timing.slice > thresh.scan_timing.zero'); -end - -if do_detect_vol_events_by_grad_height && (thresh.slice > thresh.vol) - error('Please set thresh.scan_timing.vol > thresh.scan_timing.slice'); -end - - - -Nscans = sqpar.Nscans; -Ndummies = sqpar.Ndummies; -NslicesPerBeat = sqpar.NslicesPerBeat; -Nslices = sqpar.Nslices; - -[z{1:10}]=textread(logfile,'%d %d %d %d %d %d %d %d %d %d','commentstyle', 'shell'); -y = cell2mat(z); - -Nsamples=size(y,1); - -dt = log_files.sampling_interval; - -%default: 500 Hz sampling frequency -if isempty(dt) - dt = 2e-3; -end - -t = -log_files.startScanSeconds + ((0:(Nsamples-1))*dt)'; - -LOC_END_MARKER = find( y(:,end) == 20 ); - -TA = 1/dt; - -NallVols = (Ndummies+Nscans); -VOLLOCS = zeros(NallVols,1); -LOCS = zeros(NallVols*Nslices,1); -TR = sqpar.TR; - -VOLLOCS = [LOC_END_MARKER-TR*TA : -TR*TA : LOC_END_MARKER-TR*TA*NallVols]; -VOLLOCS = VOLLOCS(end:-1:1); -LOCS = []; -for V=1:NallVols - - LOCS_PER_SLICE = sqpar.PtsSliceToSlice; - - for SL=1:Nslices - LOCS(end+1) = [VOLLOCS(V) + (SL-1)*LOCS_PER_SLICE ]; - end - -end - - LOCS = reshape(LOCS,length(LOCS),1); - VOLLOCS = reshape(VOLLOCS,length(VOLLOCS),1); - - - if verbose.level>=1 - verbose.fig_handles(end+1) = tapas_physio_get_default_fig_params(); - set(gcf,'Name', 'Thresholding Gradient for slice acq start detection'); - fs(1) = subplot(1,1,1); - plot(t, y(:,7:9)); - legend('gradient x', 'gradient y', 'gradient z'); - title('Raw Gradient Time-courses'); - hold on, - ylims = ylim; - - plot( [(VOLLOCS(1)-1)/TA (VOLLOCS(1)-1)/TA] , ylims, 'k' ) - plot( [(VOLLOCS(1+Ndummies)-1)/TA (VOLLOCS(1+Ndummies)-1)/TA] , ylims, 'g' ) - plot( [(VOLLOCS(end)-1)/TA (VOLLOCS(end)-1)/TA], ylims, 'k' ) - plot( [(LOCS(end)-1)/TA (LOCS(end)-1)/TA] , ylims, 'k' ) - plot( [(VOLLOCS(end-1)-1)/TA (VOLLOCS(end-1)-1)/TA] , ylims, 'k' ) - - plot( [(LOC_END_MARKER-1)/TA (LOC_END_MARKER-1)/TA], ylims, 'g' ) - end - - % VOLLOCS = find(abs(diff(z2))>thresh.vol); - if isempty(VOLLOCS) || isempty(LOCS) - error('No volume start events found, Decrease thresh.vol or thresh.slice after considering the Thresholding figure'); - elseif length(LOCS) < NslicesPerBeat - error('Too few slice start events found. Decrease thresh.slice after considering the Thresholding figure'); - end - diff --git a/PhysIO/code/tapas_physio_create_scan_timing_from_gradients_auto_philips.m b/PhysIO/code/tapas_physio_create_scan_timing_from_gradients_auto_philips.m deleted file mode 100644 index f6aeeaa6..00000000 --- a/PhysIO/code/tapas_physio_create_scan_timing_from_gradients_auto_philips.m +++ /dev/null @@ -1,479 +0,0 @@ -function [VOLLOCS, LOCS, verbose] = ... - tapas_physio_create_scan_timing_from_gradients_auto_philips(log_files, ... - scan_timing, verbose) -% Automatically extracts slice/volume scan events from gradient timecourse -% SCANPHYSLOG file -% -% [VOLLOCS, LOCS, verbose] = ... -% tapas_physio_create_scan_timing_from_gradients_auto_philips(log_files, ... -% scan_timing, verbose) -% -% This function determines slice/volume starts from the gradient time course -% automatically by assuming a regularity of them from the -% sequence timing parameters, in particular TR and number of slices -% -% Therefore, unlike tapas_physio_create_scan_timing_from_gradients_philips, -% no thresholds of slice/volume starts have to be given to determine the -% timing and are inferred on iteratively instead. -% -% Workflow -% 1. Determine template of each volume gradient time course by using TR. -% 2. Determine volume events using volume template (counting either from -% start or end of the time series) -% 3. Determine slice events between all detected volumes -% (again, creating a slice template and matching it to time series -% between consecutive volumes...) -% -% IN -% log_files is a structure containing the following filenames (with full -% path) -% .log_cardiac contains ECG or pulse oximeter time course -% for Philips: 'SCANPHYSLOG.log'; -% can be found on scanner in G:/log/scanphyslog- -% directory, one file is created per scan, make sure to take -% the one with the time stamp corresponding to your PAR/REC -% files -% .log_respiration contains breathing belt amplitude time course -% for Philips: same as .log_cardiac -% -% scan_timing - -% Parameters for sequence timing & synchronization -% scan_tming.sqpar = slice and volume acquisition starts, TR, -% number of scans etc. -% scan_timing.sync = synchronize phys logfile to scan acquisition -% -% sync gradient amplitude thresholds to detect slice and volume events -% -% sync is a structure with the following elements -% .zero - gradient values below this value are set to zero; -% should be those which are unrelated to slice acquisition start -% .slice - minimum gradient amplitude to be exceeded when a slice -% scan starts -% .vol - minimum gradient amplitude to be exceeded when a new -% volume scan starts; -% leave [], if volume events shall be determined as -% every Nslices-th scan event -% .grad_direction -% - leave empty to use nominal timing; -% if set, sequence timing is calculated from logged gradient timecourse; -% - value determines which gradient direction timecourse is used to -% identify scan volume/slice start events ('x', 'y', 'z') -% .vol_spacing -% - duration (in seconds) from last slice acq to -% first slice of next volume; -% leave [], if .vol-threshold shall be used -% sqpar - sequence timing parameters -% .nSlices - number of slices per volume in fMRI scan -% .nSlicesPerBeat - usually equals nSlices, unless you trigger with the heart beat -% .TR - repetition time in seconds -% .nDummies - number of dummy volumes -% .nScans - number of full volumes saved (volumes in nifti file, -% usually rows in your design matrix) -% .Nprep - number of non-dummy, volume like preparation pulses -% before 1st dummy scan. If set, logfile is read from beginning, -% otherwise volumes are counted from last detected volume in the logfile -% .time_slice_to_slice - time between the acquisition of 2 subsequent -% slices; typically TR/nSlices or -% minTR/nSlices, if minimal temporal slice -% spacing was chosen -% onset_slice - slice whose scan onset determines the adjustment of the -% regressor timing to a particular slice for the whole volume -% -% NOTE: only necessary, if sync.grad_direction is empty -% verbose -% -% OUT -% VOLLOCS - locations in time vector, when volume scan -% events started -% LOCS - locations in time vector, when slice or volume scan -% events started -% -% EXAMPLE -% [VOLLOCS, LOCS, verbose] = tapas_physio_create_scan_timing_from_gradients_philips_auto(logfile, -% scan_timing, verbose); -% -% See also tapas_physio_create_scan_timing_from_gradients_philips -% -% Author: Lars Kasper -% Created: 2015-01-09 -% Copyright (C) 2013 Institute for Biomedical Engineering, ETH/Uni Zurich. -% -% This file is part of the PhysIO toolbox, which is released under the terms of the GNU General Public -% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL -% (either version 3 or, at your option, any later version). For further details, see the file -% COPYING or . -% -% $Id: tapas_physio_create_scan_timing_from_gradients_philips.m 632 2015-01-09 12:36:12Z kasperla $ - - -% smaller than typical single shot EPI slice duration (including waiting -% for TE) - -sqpar = scan_timing.sqpar; -sync = scan_timing.sync; - -debug = verbose.level >=2 ; - -minSliceDuration = 0.040; - -doCountSliceEventsFromLogfileStart = ... - strcmpi(log_files.align_scan, 'first'); - - -% everything stored in 1 logfile -if ~isfield(log_files, 'cardiac') || isempty(log_files.cardiac) - logfile = log_files.respiration; -else - logfile = log_files.cardiac; -end - -nScans = sqpar.Nscans; -nDummies = sqpar.Ndummies; -nSlicesPerBeat = sqpar.NslicesPerBeat; -nSlices = sqpar.Nslices; - - -if doCountSliceEventsFromLogfileStart - Nprep = sqpar.Nprep; -end - -y = tapas_physio_read_physlogfiles_philips_matrix(logfile); - -acq_codes = y(:,10); -nSamples = size(y,1); - -dt = log_files.sampling_interval(1); - -%default: 500 Hz sampling frequency -if isempty(dt) - dt = 2e-3; -end - -t = -log_files.relative_start_acquisition + ((0:(nSamples-1))*dt)'; - - - -% finding scan volume starts (svolpulse) - -G = y(:,7:9); -switch lower(sync.grad_direction) - case 'x' - gradient_choice = y(:,7); - case 'y' - gradient_choice = y(:,8); - case 'z' - gradient_choice = y(:,9); - case {'xyz', 'abs'} - gradient_choice = sqrt(sum(y(:,7:9).^2,2)); - otherwise - gradient_choice = sqrt(sum(y(:,7:9).^2,2)); -end -gradient_choice = reshape(gradient_choice, [] ,1); - -if debug - tapas_physio_plot_gradient(G); -end - - -if verbose.level>=1 - % Depict all gradients, raw - [fh, fs] = plot_raw_gradients(t, y, acq_codes); - verbose.fig_handles(end+1) = fh; -end - - -% For new Ingenia log-files, recorded gradient strength may change after a -% certain time and introduce steps that are bad for recording - -minStepDistanceSamples = 1.5*ceil(sqpar.TR/dt); -gradient_choice = tapas_physio_rescale_gradient_gain_fluctuations(... - gradient_choice, minStepDistanceSamples, verbose); - - -%% 1. Create slice template -% take template from end of readout to avoid problems with initial -% values... -minVolumeDistanceSamples = ceil(sqpar.TR*0.95/dt); -minSliceDistanceSamples = ceil(minSliceDuration/dt); - -if doCountSliceEventsFromLogfileStart - rangeTemplateDetermination = 1:(nDummies + nScans)*... - minVolumeDistanceSamples; -else - rangeTemplateDetermination = (nSamples-(nDummies+nScans) * ... - minVolumeDistanceSamples + 1):nSamples; -end - -templateTime = t(rangeTemplateDetermination); -templateG = gradient_choice(rangeTemplateDetermination); -thresh_min = tapas_physio_prctile(abs(templateG), 95); - -[templateGradientSlice, secondGuessLOCS, averageTRSliceSamples] = ... - tapas_physio_get_cardiac_pulse_template(templateTime, templateG, ... - verbose, ... - 'thresh_min', thresh_min, ... - 'minCycleSamples', minSliceDistanceSamples, ... - 'shortenTemplateFactor', 0.7); - -if debug - verbose.fig_handles(end+1) = plot_template(t, templateGradientSlice); -end - - - -%% 2. Determine slice events from template using cross-correlation - -[LOCS, verbose] = tapas_physio_findpeaks_template_correlation(... - gradient_choice, templateGradientSlice, secondGuessLOCS,... - averageTRSliceSamples, verbose); - -% Template-correlation algorithm expects regular cycles, but in certain -% parts of the gradient time course (preparation, end of log-file), there -% are no cycles expected and the gradient is silent. Therefore, remove -% those ill-detected LOCS - -minTemplateAmplitude = max(abs(templateGradientSlice))*0.1; -LOCS(abs(gradient_choice(LOCS))=1 - % Plot gradient thresholding for slice timing determination - axes(fs(2)); - - % Plot gradient thresholding for slice timing determination - plot_gradients_thresholds_events(t, gradient_choice, VOLLOCS, LOCS); - linkaxes(fs(1:2), 'x'); -end - -if verbose.level>=1 - axes(fs(3)); - plot_diff_LOCS(t, LOCS, dt); - linkaxes(fs,'x'); -end - -if debug - % verbose.fig_handles(end+1) = plot_slice_events( LOCS, t, ... - % gradient_choice, templateGradientSlice, secondGuessLOCS); - % - % plot_diff_LOCS(t, LOCS, dt) -end - - - -%% Select relevant events from detected ones using sequence parameter info - - -% VOLLOCS-detection via spacing or counting -nLocs = numel(LOCS); - -if isfield(sync, 'vol_spacing') && ~isempty(sync.vol_spacing) - iVolLocs = find((diff(LOCS) > sync.vol_spacing/dt)) + 1; -else - - - if doCountSliceEventsFromLogfileStart - iVolLocs = Nprep*nSlices + (1:nSlices:(nDummies+nScans)*nSlices); - else % count from end - iVolLocs = (nLocs-(nDummies+nScans)*nSlices+1):nSlices:nLocs; - end - -end - -VOLLOCS = LOCS(iVolLocs(iVolLocs<=nLocs)); -LOCS = reshape(LOCS, [], 1); -VOLLOCS = reshape(VOLLOCS, [], 1); - - -if verbose.level>=1 - % Plot gradient thresholding for slice timing determination - axes(fs(2)); - - % Plot gradient thresholding for slice timing determination - plot_gradients_thresholds_events(t, gradient_choice, VOLLOCS, LOCS); - linkaxes(fs, 'x'); -end - -%% Repair scan event problems... - -% remove erroneous volume events, i.e. those due to slice gaps that are -% assumed to be volume end gaps, and create new volumes with not enough -% slices -minVolumeDistanceSamplesError = ceil((1-1/sqpar.Nslices) * sqpar.TR/dt); -idxVolError = find(diff(VOLLOCS) < minVolumeDistanceSamplesError); -if ~isempty(idxVolError) - VOLLOCS(idxVolError(2:2:end)) = []; -end - -ons.acq_slice_all = LOCS; -ons.vol_all = VOLLOCS; -[any_scanevent_repaired, ons] = ... - tapas_physio_repair_scan_events_PHILIPS(ons, sqpar, verbose.level > 2); - -LOCS = ons.acq_slice_all; -VOLLOCS = ons.vol_all; - -%% Return error if not enough events flund -% VOLLOCS = find(abs(diff(z2))>sync.vol); -if isempty(VOLLOCS) || isempty(LOCS) - error('No volume start events found, Decrease sync.vol or sync.slice after considering the Thresholding figure'); -elseif length(LOCS) < nSlicesPerBeat - error('Too few slice start events found. Decrease sync.slice after considering the Thresholding figure'); -end - -if doCountSliceEventsFromLogfileStart - if length(VOLLOCS)< (Nprep+nScans+nDummies) - error(['Not enough volume events found. \n\tFound: %d\n ' ... - '\tNeeded: %d+%d+%d (Nprep+nDummies+nScans)\n' ... - 'Please lower sync.vol or sync.vol_spacing\n'], ... - length(VOLLOCS), Nprep, nDummies, nScans); - end -else - if length(VOLLOCS)< (nScans+nDummies) - error(['Not enough volume events found. \n\tFound: %d\n ' ... - '\tNeeded: %d+%d (nDummies+nScans)\n' ... - 'Please lower sync.vol or sync.vol_spacing\n'], ... - length(VOLLOCS), nDummies, nScans); - end -end - -end - -%% local functions for debugging plots - -%% Plot template for volume repetition -function fh = plot_template(t, templateGradientVolume) - -stringTitle = 'Template Gradient Timecourse during 1 Slice'; -fh = tapas_physio_get_default_fig_params(); -set(gcf, 'Name', stringTitle); - -nSamplesTemplate = numel(templateGradientVolume); -plot(t(1:nSamplesTemplate), templateGradientVolume); -xlabel('t (s)') -title(stringTitle); -end - - -%% Plot Detected volume events -function fh = plot_volume_events(VOLLOCS, t, G, ... - templateGradientVolume, secondGuessVOLLOCS) - -stringTitle = 'Template Gradient Timecourse during 1 Volume'; -fh = tapas_physio_get_default_fig_params(); -set(gcf, 'Name', stringTitle); - -ampl = max(abs(G)); -sG = conv(G, templateGradientVolume/sum(abs(templateGradientVolume)), ... - 'same'); - -plot(t, G); hold all; -plot(t, sG); -stem(t(VOLLOCS), ampl*ones(size(VOLLOCS))); -stem(t(secondGuessVOLLOCS), ampl*ones(size(secondGuessVOLLOCS))); - -stringLegend = { - 'Gradient Timecourse' - 'Gradient Timecourse convolved with Template' - 'Final detected volume events' - 'Prior detected volume events' - }; - -xlabel('t (s)') -title(stringTitle); -legend(stringLegend) -end - - -%% Plot Detected slice events -function fh = plot_slice_events(LOCS, t, G, ... - templateGradientSlice, secondGuessLOCS) - -stringTitle = 'Template Gradient Timecourse during 1 Slice'; -fh = tapas_physio_get_default_fig_params(); -set(gcf, 'Name', stringTitle); - -ampl = max(abs(G)); -sG = conv(G, templateGradientSlice/sum(abs(templateGradientSlice)), ... - 'same'); - -plot(t, G); hold all; -plot(t, sG); -stem(t(LOCS), ampl*ones(size(LOCS))); -stem(t(secondGuessLOCS), ampl*ones(size(secondGuessLOCS))); - -stringLegend = { - 'Gradient Timecourse' - 'Gradient Timecourse convolved with Template' - 'Final detected slice events' - 'Prior detected slice events' - }; - -xlabel('t (s)') -title(stringTitle); -legend(stringLegend) -end - - -%% Plot difference between detected slice events -function plot_diff_LOCS(t, LOCS, dt) - - -dLocsSecs = diff(LOCS)*dt*1000; -ymin = tapas_physio_prctile(dLocsSecs, 25); -ymax = tapas_physio_prctile(dLocsSecs, 99); - - -plot(t(LOCS(1:end-1)), dLocsSecs); -title('duration between scan events - search for bad peaks here!'); -xlabel('t (s)'); -ylabel('t (ms)'); -ylim([0.9*ymin, 1.1*ymax]); - -end - -%% Plot all raw gradient time-courses -function [fh, fs] = plot_raw_gradients(t, y, acq_codes) -fh = tapas_physio_get_default_fig_params(); -set(gcf,'Name', 'Thresholding Gradient for slice acq start detection'); -for i=1:3, fs(i) = subplot(3,1,i); end - -axes(fs(1)); - -plot(t, sqrt(sum(y(:,7:9).^2,2)), '--k'); -hold all; -plot(t, y(:,7:9)); - - -if ismember(8,acq_codes) - hold all; - stem(t, acq_codes*max(max(abs(y(:,7:9))))/20); -end - - -legend('abs(G_x^2+G_y^2+G_z^2)', 'gradient x', 'gradient y', 'gradient z'); -title('Raw Gradient Time-courses'); -end - - -% Plot gradient thresholding for slice timing determination -function plot_gradients_thresholds_events(t, gradient_choice, VOLLOCS, LOCS) - -hp = plot(t,gradient_choice); hold all; -lg = {'Chosen gradient for thresholding'}; - -title({'Thresholding Gradient for slice acq start detection', '- found scan events -'}); -legend(hp, lg); -xlabel('t(s)'); - -if ~isempty(VOLLOCS) - hp(end+1) = stem(t(VOLLOCS), 1.25*ones(size(VOLLOCS))); hold all - lg{end+1} = sprintf('Found volume events (N = %d)', numel(VOLLOCS)); -end - -if ~isempty(LOCS) - hp(end+1) = stem(t(LOCS), ones(size(LOCS))); hold all - lg{end+1} = sprintf('Found slice events (N = %d)', numel(LOCS)); -end -legend(hp, lg); -end \ No newline at end of file diff --git a/PhysIO/code/tapas_physio_create_scan_timing_from_gradients_philips.m b/PhysIO/code/tapas_physio_create_scan_timing_from_gradients_philips.m deleted file mode 100644 index e47145ec..00000000 --- a/PhysIO/code/tapas_physio_create_scan_timing_from_gradients_philips.m +++ /dev/null @@ -1,316 +0,0 @@ -function [VOLLOCS, LOCS, verbose] = ... - tapas_physio_create_scan_timing_from_gradients_philips(log_files, ... - scan_timing, verbose) -% Extracts slice and volume scan events from gradients timecourse of Philips -% SCANPHYSLOG file -% -% [VOLLOCS, LOCS] = ... -% tapas_physio_create_scan_timing_from_gradients_philips(logfile, -% scan_timing, verbose); -% -% IN -% log_files is a structure containing the following filenames (with full -% path) -% .log_cardiac contains ECG or pulse oximeter time course -% for Philips: 'SCANPHYSLOG.log'; -% can be found on scanner in G:/log/scanphyslog- -% directory, one file is created per scan, make sure to take -% the one with the time stamp corresponding to your PAR/REC -% files -% .log_respiration contains breathing belt amplitude time course -% for Philips: same as .log_cardiac -% -% scan_timing - -% Parameters for sequence timing & synchronization -% scan_tming.sqpar = slice and volume acquisition starts, TR, -% number of scans etc. -% scan_timing.sync = synchronize phys logfile to scan acquisition -% -% sync gradient amplitude thresholds to detect slice and volume events -% -% sync is a structure with the following elements -% .zero - gradient values below this value are set to zero; -% should be those which are unrelated to slice acquisition start -% .slice - minimum gradient amplitude to be exceeded when a slice -% scan starts -% .vol - minimum gradient amplitude to be exceeded when a new -% volume scan starts; -% leave [], if volume events shall be determined as -% every Nslices-th scan event -% .grad_direction -% - leave empty to use nominal timing; -% if set, sequence timing is calculated from logged gradient timecourse; -% - value determines which gradient direction timecourse is used to -% identify scan volume/slice start events ('x', 'y', 'z') -% .vol_spacing -% - duration (in seconds) from last slice acq to -% first slice of next volume; -% leave [], if .vol-threshold shall be used -% -% sqpar - sequence timing parameters -% .Nslices - number of slices per volume in fMRI scan -% .TR - repetition time in seconds -% .Ndummies - number of dummy volumes -% .Nscans - number of full volumes saved (volumes in nifti file, -% usually rows in your design matrix) -% .Nprep - number of non-dummy, volume like preparation pulses -% before 1st dummy scan. If set, logfile is read from beginning, -% otherwise volumes are counted from last detected volume in the logfile -% .time_slice_to_slice - time between the acquisition of 2 subsequent -% slices; typically TR/Nslices or -% minTR/Nslices, if minimal temporal slice -% spacing was chosen -% onset_slice - slice whose scan onset determines the adjustment of the -% regressor timing to a particular slice for the whole volume -% -% NOTE: only necessary, if sync.grad_direction is empty -% verbose -% -% OUT -% VOLLOCS - locations in time vector, when volume scan -% events started -% LOCS - locations in time vector, when slice or volume scan -% events started -% -% EXAMPLE -% [VOLLOCS, LOCS, verbose] = tapas_physio_create_scan_timing_from_gradients_philips(logfile, -% scan_timing, verbose); -% -% See also -% -% Author: Lars Kasper -% Created: 2013-02-16 -% Copyright (C) 2013 Institute for Biomedical Engineering, ETH/Uni Zurich. -% -% This file is part of the PhysIO toolbox, which is released under the terms of the GNU General Public -% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL -% (either version 3 or, at your option, any later version). For further details, see the file -% COPYING or . -% -% $Id$ - -sqpar = scan_timing.sqpar; -sync = scan_timing.sync; - -% smaller than typical single shot EPI slice duration (including waiting -% for TE) -minSliceDuration = 0.040; - -doDetectVolumesByCounting = (~isfield(sync, 'vol') || ... - isempty(sync.vol)) && (~isfield(sync, 'vol_spacing') || ... - isempty(sync.vol_spacing)); -doDetectVolumesByGradientAmplitude = ~doDetectVolumesByCounting && ... - (~isfield(sync, 'vol_spacing') || isempty(sync.vol_spacing)); -doCountSliceEventsFromLogfileStart = ... - strcmpi(log_files.align_scan, 'first'); - - -%% check consistency of sync-values - -if sync.slice <= sync.zero - verbose = tapas_physio_log('Please set sync.scan_timing.slice > sync.scan_timing.zero', ... - verbose, 2); -end - -if doDetectVolumesByGradientAmplitude && (sync.slice > sync.vol) - verbose = tapas_physio_log('Please set sync.scan_timing.vol > sync.scan_timing.slice', ... - verbose, 2); -end - - - -% everything stored in 1 logfile -if ~isfield(log_files, 'cardiac') || isempty(log_files.cardiac) - logfile = log_files.respiration; -else - logfile = log_files.cardiac; -end - -Nscans = sqpar.Nscans; -Ndummies = sqpar.Ndummies; -Nslices = sqpar.Nslices; - - -if doCountSliceEventsFromLogfileStart - Nprep = sqpar.Nprep; -end - -y = tapas_physio_read_physlogfiles_philips_matrix(logfile); - -acq_codes = y(:,10); -nSamples = size(y,1); - -dt = log_files.sampling_interval(1); - -%default: 500 Hz sampling frequency -if isempty(dt) - dt = 2e-3; -end - -t = -log_files.relative_start_acquisition + ((0:(nSamples-1))*dt)'; - - - -% finding scan volume starts (svolpulse) -switch lower(sync.grad_direction) - case 'x' - gradient_choice = y(:,7); - case 'y' - gradient_choice = y(:,8); - case 'z' - gradient_choice = y(:,9); - case {'xyz', 'abs'} - gradient_choice = sqrt(sum(y(:,7:9).^2,2)); -end -gradient_choice = reshape(gradient_choice, [] ,1); - -% For new Ingenia log-files, recorded gradient strength may change after a -% certain time and introduce steps that are bad for recording - -minStepDistanceSamples = ceil(1.5*sqpar.TR/dt); - -% Normalize gradients, if thresholds are smaller than 1, i.e. relative -doNormalize = max([sync.slice, sync.vol, sync.zero]) < 1; - -[gradient_choice, gainArray, normFactor] = ... - tapas_physio_rescale_gradient_gain_fluctuations(... - gradient_choice, minStepDistanceSamples, verbose, 'doNormalize', ... - doNormalize); - - - -% if no gradient timecourse was recorded in the logfile (due to insufficient -% Philips software keys), return nominal timing instead -if ~any(gradient_choice) % all values zero - [VOLLOCS, LOCS] = tapas_physio_create_scan_timing_nominal(t, sqpar); - verbose = tapas_physio_log('No gradient timecourse was logged in the logfile. Using nominal timing from sqpar instead', ... - verbose, 1); - return -end - -z2 = gradient_choice; z2(z2 sync.vol_spacing/dt)) + 1); - end - end - LOCS = reshape(LOCS, [], 1); - VOLLOCS = reshape(VOLLOCS, [], 1); -catch - VOLLOCS = []; -end - -if verbose.level>=1 - - % Depict all gradients, raw - verbose.fig_handles(end+1) = tapas_physio_get_default_fig_params(); - set(gcf,'Name', 'Thresholding Gradient for slice acq start detection'); - fs(1) = subplot(3,1,1); - - plot(t, sqrt(sum(y(:,7:9).^2,2)), '--k'); - hold all; - plot(t, y(:,7:9)); - - - if ismember(8,acq_codes) - hold all; - ampl = max(abs(y(~isinf(y)))); - stem(t, ampl/20*acq_codes); - end - - - legend('abs(G_x^2+G_y^2+G_z^2)', 'gradient x', 'gradient y', 'gradient z'); - title('Raw Gradient Time-courses'); - - % Plot gradient thresholding for slice timing determination - fs(2) = subplot(3,1,2); - hp = plot(t,[gradient_choice z2]); hold all; - hp(end+1) = plot(t, repmat(sync.zero, nSamples, 1)); - hp(end+1) = plot(t, repmat(sync.slice, nSamples, 1)); - lg = {'Chosen gradient for thresholding', ... - 'Gradient with values < sync.zero set to 0', ... - 'sync.zero', 'sync.slice'}; - - if doDetectVolumesByGradientAmplitude - hp(end+1) = plot(t, repmat(sync.vol, nSamples, 1)); - lg{end+1} = 'sync.vol'; - end - title({'Thresholding Gradient for slice acq start detection', '- found scan events -'}); - legend(hp, lg); - xlabel('t(s)'); - - % Plot gradient thresholding for slice timing determination - - - if ~isempty(VOLLOCS) - hp(end+1) = stem(t(VOLLOCS), 1.25*ones(size(VOLLOCS))); hold all - lg{end+1} = sprintf('Found volume events (N = %d)', numel(VOLLOCS)); - end - - if ~isempty(LOCS) - hp(end+1) = stem(t(LOCS), ones(size(LOCS))); hold all - lg{end+1} = sprintf('Found slice events (N = %d)', numel(LOCS)); - - dLocsSecs = diff(LOCS)*dt*1000; - ymin = tapas_physio_prctile(dLocsSecs, 25); - ymax = tapas_physio_prctile(dLocsSecs, 99); - - fs(3) = subplot(3,1,3); - plot(t(LOCS(1:end-1)), dLocsSecs); title('duration between scan events - search for bad peaks here!'); - xlabel('t (s)'); - ylabel('t (ms)'); - ylim([0.9*ymin, 1.1*ymax]); - linkaxes(fs,'x'); - - end - subplot(3,1,2); - legend(hp, lg); - -end - - -%% Return error if not enough events flund -% VOLLOCS = find(abs(diff(z2))>sync.vol); -if isempty(VOLLOCS) || isempty(LOCS) - verbose = tapas_physio_log('No volume start events found, Decrease sync.vol or sync.slice after considering the Thresholding figure', ... - verbose, 2); -elseif length(LOCS) < Nslices - verbose = tapas_physio_log('Too few slice start events found. Decrease sync.slice after considering the Thresholding figure', ... - verbose, 2); -end - -if doCountSliceEventsFromLogfileStart - if length(VOLLOCS)< (Nprep+Nscans+Ndummies) - verbose = tapas_physio_log(sprintf(['Not enough volume events found. \n\tFound: %d\n ' ... - '\tNeeded: %d+%d+%d (Nprep+Ndummies+Nscans)\n' ... - 'Please lower sync.vol or sync.vol_spacing\n'], ... - length(VOLLOCS), Nprep, Ndummies, Nscans), verbose, 2); - end -else - if length(VOLLOCS)< (Nscans+Ndummies) - verbose = tapas_physio_log(sprintf(['Not enough volume events found. \n\tFound: %d\n ' ... - '\tNeeded: %d+%d (Ndummies+Nscans)\n' ... - 'Please lower sync.vol or sync.vol_spacing\n'], ... - length(VOLLOCS), Ndummies, Nscans), verbose, 2); - end -end - diff --git a/PhysIO/code/tapas_physio_create_scan_timing_from_tics_siemens.m b/PhysIO/code/tapas_physio_create_scan_timing_from_tics_siemens.m deleted file mode 100644 index 573afb67..00000000 --- a/PhysIO/code/tapas_physio_create_scan_timing_from_tics_siemens.m +++ /dev/null @@ -1,83 +0,0 @@ -function [VOLLOCS, LOCS, verbose] = ... - tapas_physio_create_scan_timing_from_tics_siemens(t, log_files, verbose) -% Creates locations of scan volume and slice events in physiological time series vector from Siemens Tics-file -% (_