diff --git a/.gitignore b/.gitignore index 365f7d40..8dfde4cc 100644 --- a/.gitignore +++ b/.gitignore @@ -18,6 +18,7 @@ lib .idea .idea/* examples +.MATLABDriveTag # Python # See e.g. https://github.com/github/gitignore/blob/master/Python.gitignore diff --git a/CHANGELOG.md b/CHANGELOG.md index 4823565d..4904231a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,20 @@ # Changelog TAPAS toolbox +## [6.1.0] 2022-11-18 +### Added +- Added `tapas_test` testing function. Started progressive introduction of testing framework. +- PhysIO [Version 8.2.0](PhysIO/CHANGELOG.md) + - Interface `tapas_physio_test` to TAPAS-generic `tapas_test` function + - Added suport for logfile version 3 of Siemens physio recordings + - multi ECG/Resp channels and interleaved status messages + - new integration test for Siemens VB Logversion 3 +### Fixed +- PhysIO [Version 8.2.0](PhysIO/CHANGELOG.md) + - Removed dependence on `nanmean` (Statistics Toolbox) + - See [GitHub issue #205](https://github.com/translationalneuromodeling/tapas/issues/205) + - Compatibility with multiple SPM toolbox locations for `lmod` ([GitHub issue #211](https://github.com/translationalneuromodeling/tapas/issues/211)) + - as listed in `spm_get_defaults('tbx')` + - Refactoring of Philips read-in to support novel 12-column logfile version, see [GitHub issue #207](https://github.com/translationalneuromodeling/tapas/issues/207#issuecomment-1246078600) ## [6.0.1] 2022-05-17 @@ -11,7 +26,7 @@ For this major release of TAPAS, we have decided to move currently unmaintained to a dedicated [TAPAS Legacy repository](https://github.com/translationalneuromodeling/tapas_legacy). ### Added -- PhysIO [Version 8.1.0](Physio/CHANGELOG.md): BIDS Read-in for separate cardiac/respiratory trace files (e.g., due to different sampling frequencies) +- PhysIO [Version 8.1.0](PhysIO/CHANGELOG.md): BIDS Read-in for separate cardiac/respiratory trace files (e.g., due to different sampling frequencies) - Compatibility of whole code base with Matlab compiler in order to run `spm_make_standalone`. ### Changed - Moved `MICP`,`VBLM`,`h2gf` and `mpdcm` toolboxes to TAPAS Legacy repository. diff --git a/Contributor-License-Agreement.md b/Contributor-License-Agreement.md index 3d6e2584..ab770aeb 100644 --- a/Contributor-License-Agreement.md +++ b/Contributor-License-Agreement.md @@ -21,6 +21,7 @@ Stephan Heunis | Eindhoven University of Technology | Eindhoven | NL Niklas Bürgi | SNS, University of Zurich | Zurich | CH | nbuergi | 1.0 Alexandre Sayal | CIBIT, University of Coimbra | Coimbra | PT | alexsayal | 1.0 Matthias Müller-Schrader | TNU, University of Zurich | Zurich | CH | mms-neuro | 1.1 +Johanna M. M. Bayer | The University of Melbourne | Melbourne | AU | likeajumprope | 1.1 Saskia Bollmann | The University of Queensland | Brisbane | AUS | SaskiaBollmann | 1.1 **- Add Entry here -** | **- Add Entry here -** | **Add** | **Add** | **Add** | 1.1 diff --git a/PhysIO/CHANGELOG.md b/PhysIO/CHANGELOG.md index f44a9145..1d07a81d 100644 --- a/PhysIO/CHANGELOG.md +++ b/PhysIO/CHANGELOG.md @@ -4,10 +4,29 @@ RELEASE INFORMATION Current Release --------------- -*Current version: PhysIO Toolbox Release R2022a, v8.1.0* +*Current version: PhysIO Toolbox Release R2022b, v8.2.0* -April 5th, 2022 +November 22nd, 2022 +Upcoming Release Notes (v8.2.0-beta) +------------------------------------ + +### Added +- Interface `tapas_physio_test` to TAPAS-generic `tapas_test` function +- Added suport for logfile version 3 of Siemens physio recordings + - multi ECG/Resp channels and interleaved status messages + - new integration test for Siemens VB Logversion 3 +- Added support for ADInstruments/LabChart Txt-export format (see + [CUBRIC Seminar Example](https://github.com/BRAIN-TO/cubric-physio) and + gitlab branch #107) + +### Fixed +- Removed dependence on `nanmean` (Statistics Toolbox) + - See [GitHub issue #205](https://github.com/translationalneuromodeling/tapas/issues/205) +- Compatibility with multiple SPM toolbox locations for `lmod` ([GitHub issue #211](https://github.com/translationalneuromodeling/tapas/issues/211)) + - as listed in `spm_get_defaults('tbx')` +- Refactoring of Philips read-in to support novel 12-column logfile version, see [GitHub issue #207](https://github.com/translationalneuromodeling/tapas/issues/207#issuecomment-1246078600) +- Unit/Integration tests for filtered traces (cardiac and respiratory) switched to absolute tolerances (relative problematic for traces close to zero) Minor Release Notes (v8.1.0) ---------------------------- diff --git a/PhysIO/README.md b/PhysIO/README.md index a4ef3bf0..53b5d736 100644 --- a/PhysIO/README.md +++ b/PhysIO/README.md @@ -1,7 +1,7 @@ TAPAS PhysIO Toolbox ==================== -*Current version: Release 2022a, v8.1.0* +*Current version: Release 2022b, v8.2.0* > Copyright (C) 2012-2022 > Lars Kasper @@ -246,7 +246,7 @@ Contributors ------------ - Lead Programmer: - - [Lars Kasper](https://www.tnu.ethz.ch/en/team/faculty-and-scientific-staff/kasper.html), + - [Lars Kasper](https://brain-to.ca/content/lars-kasper), TNU & MR-Technology Group, IBT, University of Zurich & ETH Zurich - Project Team: - Steffen Bollmann, Centre for Advanced Imaging, University of Queensland, Australia diff --git a/PhysIO/code/assess/tapas_physio_review.m b/PhysIO/code/assess/tapas_physio_review.m index 5e41f6cd..859a0f2e 100644 --- a/PhysIO/code/assess/tapas_physio_review.m +++ b/PhysIO/code/assess/tapas_physio_review.m @@ -41,8 +41,8 @@ % % See also -% Author: Lars Kasper -% Created: 2016-10-27 +% Author: Johanna Bayer, Lars Kasper +% Created: 2016-10-27, completed 2023 % Copyright (C) 2016 TNU, Institute for Biomedical Engineering, % University of Zurich and ETH Zurich. % @@ -72,6 +72,7 @@ sync = scan_timing.sync; model = physio.model; verbose = physio.verbose; +review = physio.verbose.review; % Compatibility with old versions if ~isfield(model, 'R_column_names') @@ -118,28 +119,126 @@ verbose.close_figs = false; end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Figure: Raw data +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + verbose = tapas_physio_plot_raw_physdata(ons_secs.raw, verbose); -% tapas_physio_get_onsets_from_locs -> create plot function out of -% sub-function -% tapas_physio_get_cardiac_pulses_auto_matched -> subfunction for plot, only called -% if in this sub-branch + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Figure: Peak detection +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +if ismember(preproc.cardiac.initial_cpulse_select.method, {'auto','auto_template', 'auto_matched'}) + + if verbose.level >=2 + if isfield(review, 'peak') + + [verbose] = tapas_physio_plot_peak_detection_from_automatically_generated(review.peak.t, review.peak.c, ... + review.peak.cpulse, verbose); + end + end +end + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Figure: Iterative template +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% First and second figure + +if verbose.level >= 3 + + if isfield(review, 'iter_temp') + [verbose] = tapas_physio_plot_iterative_template_creation(review.iter_temp.hasFirstGuessPeaks,... + review.iter_temp.t, review.iter_temp.c, review.iter_temp.cpulse1stGuess, review.iter_temp.nPulses1, ... + review.iter_temp.nPulses2, review.iter_temp.cpulse2ndGuess, review.iter_temp.meanLag1, ... + review.iter_temp.meanLag2, verbose); + end + +% thrid figure + if isfield(review, 'temp_cyc') + [verbose] = tapas_physio_plot_templates_of_cycle_time(review.temp_cyc.tTemplate, ... + review.temp_cyc.template, review.temp_cyc.pulseTemplate, ... + review.temp_cyc.pulseCleanedTemplate, verbose) + end +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Figure: Sync Bundles +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +if verbose.level >= 3 + if isfield(review, 'sync_bundles') + + [verbose] = tapas_physio_plot_sync_bundles(review.sync_bundles.Nallvols, review.sync_bundles.t, ... + review.sync_bundles.SLICELOCS, verbose); + end +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Figure: Get cardiac +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +if verbose.level >= 3 + + if isfield(review, 'get_cardiac') + + [verbose] = tapas_physio_plot_get_cardiac_phase(review.get_cardiac.scannert, ... + review.get_cardiac.cardiac_phase, review.get_cardiac.pulset, ... + review.get_cardiac.svolpulse, verbose); + end + +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Figure: Preproc Coutcout actual scans - all events and gradients +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if verbose.level >= 2 verbose.fig_handles(end+1) = ... tapas_physio_plot_cropped_phys_to_acqwindow(ons_secs, sqpar, verbose); end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Figure: Preproc Respiratory filtering +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +if verbose.level>=3 + + if isfield(review, 'resp_filter') + + [verbose] = tapas_physio_plot_filter_respiratory(review.resp_filter.rpulset, ... + review.resp_filter.m, review.resp_filter.s, review.resp_filter.t, ... + review.resp_filter.rpulset_out, review.resp_filter.rpulset_out_trend,... + review.resp_filter.trend,review.resp_filter.rpulset_out_trend_filt, verbose); + + end +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Figure: Preproc Diagnostics for raw physiological time series +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + [verbose, ons_secs.c_outliers_low, ons_secs.c_outliers_high, ... ons_secs.r_hist] = ... tapas_physio_plot_raw_physdata_diagnostics(ons_secs.cpulse, ... ons_secs.r, preproc.cardiac.posthoc_cpulse_select, verbose, ... ons_secs.t, ons_secs.c); -% in tapas_physio_create_retroicor_regressors: -% tapas_physio_get_respiratory_phase -% function fh = plot_traces(pulset, rsampint, rout, resp_max, ... -% cumsumh, sumh, h, npulse, dpulse, rphase) + +if verbose.level >=3 + + if isfield(review, 'traces') + fh = tapas_physio_plot_traces(review.traces.pulset, review.traces.rsampint, ... + review.traces.rout, review.traces.resp_max, review.traces.cumsumh,... + review.traces.sumh, review.traces.h, review.traces.npulse, review.traces.dpulse, ... + review.traces.r_phase); + end + +end %% RETROICOR @@ -177,9 +276,25 @@ verbose.fig_handles(end+1) = tapas_physio_plot_rvt(ons_secs, sqpar); end -%% tapas_physio_create_hrv_regressors, tapas_physio_create_rvt_regressors -% tapas_physio_create_noise_rois_regressors -% => create functions out of inline-plotting +if verbose.level >= 2 + if model.rvt.include + [verbose] = tapas_physio_plot_rvt_hilbert(review.rvt_hilbert.t,review.rvt_hilbert.fr, ... + review.rvt_hilbert.fr_lp, review.rvt_hilbert.fr_mag, review.rvt_hilbert.fr_rv, ... + review.rvt_hilbert.fr_phase, review.rvt_hilbert.fr_if, verbose); + end +end + +%% tapas_physio_create_hrv_regressors, +if verbose.level>=2 + if model.rvt.include + [verbose] = tapas_physio_plot_create_hrv_regressors(review.create_hrv_regressors.sample_points, ... + review.create_hrv_regressors.hrOut, review.create_hrv_regressors.hr, review.create_hrv_regressors.t, ... + review.create_hrv_regressors.crf, review.create_hrv_regressors.convHRV, ... + review.create_hrv_regressors.delays,review.create_hrv_regressors.samplePointsOut,... + review.create_hrv_regressors.convHRVOut, verbose); + end +end + %% Overall regressors diff --git a/PhysIO/code/code b/PhysIO/code/code new file mode 100644 index 00000000..cdaa36a3 Binary files /dev/null and b/PhysIO/code/code differ diff --git a/PhysIO/code/model/tapas_physio_create_hrv_regressors.m b/PhysIO/code/model/tapas_physio_create_hrv_regressors.m index f89f7ddb..b5c7cbf0 100644 --- a/PhysIO/code/model/tapas_physio_create_hrv_regressors.m +++ b/PhysIO/code/model/tapas_physio_create_hrv_regressors.m @@ -55,40 +55,16 @@ 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', 'Model: 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 - - % Generate CRF dt = sqpar.TR/sqpar.Nslices; t = 0:dt:30; % seconds crf = tapas_physio_crf(t); crf = crf / max(abs(crf)); -if verbose.level>=2 - subplot(2,2,2) - plot(t, crf,'r');xlabel('time (seconds)'); - title('Cardiac response function'); -end - - % Convolve and rescale for display purposes convHRV = tapas_physio_conv(hr, crf, 'causal'); 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 (convolution) % @@ -101,7 +77,7 @@ % remove mean and linear trend to fulfill periodicity condition for % shifting -convHRV = detrend(convHRV); +convHRV_detrend = detrend(convHRV); % TODO: what happens at the end/beginning of shifted convolutions? nDelays = numel(delays); @@ -116,7 +92,7 @@ samplePointsOut = zeros(nScans,nSampleSlices); for iDelay = 1:nDelays - convHRVShifted = circshift(convHRV, nShiftSamples(iDelay)); + convHRVShifted = circshift(convHRV_detrend, nShiftSamples(iDelay)); for iSlice = 1:nSampleSlices onset_slice = sqpar.onset_slice(iSlice); hrOut(:,iSlice) = hr(onset_slice:sqpar.Nslices:end)'; @@ -125,14 +101,21 @@ end end +% save relevant structures +verbose.review.create_hrv_regressors.sample_points = sample_points; +verbose.review.create_hrv_regressors.hrOut = hrOut; +verbose.review.create_hrv_regressors.hr = hr; +verbose.review.create_hrv_regressors.t = t; +verbose.review.create_hrv_regressors.crf = crf; +verbose.review.create_hrv_regressors.convHRV = convHRV; +verbose.review.create_hrv_regressors.delays = delays; +verbose.review.create_hrv_regressors.samplePointsOut = samplePointsOut; +verbose.review.create_hrv_regressors.convHRVOut = convHRVOut; + 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'); + [verbose] = tapas_physio_plot_create_hrv_regressors(sample_points, hrOut, ... + hr, t, crf, convHRV, delays,samplePointsOut, convHRVOut, verbose) + end end \ No newline at end of file diff --git a/PhysIO/code/model/tapas_physio_create_noise_rois_regressors.m b/PhysIO/code/model/tapas_physio_create_noise_rois_regressors.m index 4b72dc22..1fbee3f2 100644 --- a/PhysIO/code/model/tapas_physio_create_noise_rois_regressors.m +++ b/PhysIO/code/model/tapas_physio_create_noise_rois_regressors.m @@ -178,18 +178,7 @@ [fpRoi,fnRoi] = fileparts(Vroi.fname); Vroi.fname = fullfile(fpRoi, sprintf('noiseROI_%s.nii', fnRoi)); spm_write_vol(Vroi,roi); - - % Overlay the final noise ROI (code from spm_orthviews:add_c_image) - if verbose.level >= 2 - spm_orthviews('addcolouredimage',r,Vroi.fname ,[1 0 0]); - hlabel = sprintf('%s (%s)',Vroi.fname ,'Red'); - c_handle = findobj(findobj(st.vols{r}.ax{1}.cm,'label','Overlay'),'Label','Remove coloured blobs'); - ch_c_handle = get(c_handle,'Children'); - set(c_handle,'Visible','on'); - uimenu(ch_c_handle(2),'Label',hlabel,'ForegroundColor',[1 0 0],... - 'Callback','c = get(gcbo,''UserData'');spm_orthviews(''context_menu'',''remove_c_blobs'',2,c);'); - spm_orthviews('redraw') - end + final_roi_files{r} = Vroi.fname; Yroi = Yimg(roi(:)==1, :); % Time series of the fMRI volume in the noise ROIs @@ -274,7 +263,7 @@ % plot if verbose.level >=2 - stringFig = sprintf('Model: Noise\\_rois: Extracted principal components for ROI %d', r); + stringFig = sprintf('Model: 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); @@ -299,6 +288,29 @@ R_noise_rois = [R_noise_rois, R]; -end % nROI +end % nROIs + + +%% Summary plot noise roise location +% Overlay intial and final noise ROI (code from spm_orthviews:add_c_image) +% (before/after reslice, threshold and erosion) +if verbose.level >= 2 + spm_check_registration( roi_files{:} ); + spm_orthviews('context_menu','interpolation',3); % disable interpolation // 3->NN , 2->Trilin , 1->Sinc + for r = 1:nRois + Vroi.fname = final_roi_files{r}; + spm_orthviews('addcolouredimage',r,Vroi.fname ,[1 0 0]); + hlabel = sprintf('%s (%s)',Vroi.fname ,'Red'); + + if isprop(st.vols{r}.ax{1}, 'cm') + c_handle = findobj(findobj(st.vols{r}.ax{1}.cm,'label','Overlay'),'Label','Remove coloured blobs'); + ch_c_handle = get(c_handle,'Children'); + set(c_handle,'Visible','on'); + uimenu(ch_c_handle(2),'Label',hlabel,'ForegroundColor',[1 0 0],... + 'Callback','c = get(gcbo,''UserData'');spm_orthviews(''context_menu'',''remove_c_blobs'',2,c);'); + end + end + spm_orthviews('redraw') +end end % function diff --git a/PhysIO/code/model/tapas_physio_create_retroicor_regressors.m b/PhysIO/code/model/tapas_physio_create_retroicor_regressors.m index 23b5b711..4c84cba3 100644 --- a/PhysIO/code/model/tapas_physio_create_retroicor_regressors.m +++ b/PhysIO/code/model/tapas_physio_create_retroicor_regressors.m @@ -105,8 +105,21 @@ fr = ons_secs.fr; + % save variables in verbose + verbose.review.traces.pulset = pulset; + verbose.review.traces.rsampint = rsampint; + verbose.review.traces.rout = rout; + verbose.review.traces.resp_max = resp_max; + verbose.review.traces.cumsumh = cumsumh; + verbose.review.traces.sumh = sumh; + verbose.review.traces.h = h; + verbose.review.traces.npulse = npulse; + verbose.review.traces.dpulse = dpulse; + verbose.review.traces.r_phase = r_phase; + if verbose.level >=3 - [r_phase, verbose.fig_handles(end+1)] = ... + [pulset, rsampint, rout, resp_max, cumsumh, sumh, h, ... + npulse, dpulse, r_phase, verbose.fig_handles(end+1)] = ... tapas_physio_get_respiratory_phase( ... fr,rsampint, verbose.level); else @@ -122,7 +135,8 @@ 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; diff --git a/PhysIO/code/model/tapas_physio_rvt_hilbert.m b/PhysIO/code/model/tapas_physio_rvt_hilbert.m index bb2a4697..c8993844 100644 --- a/PhysIO/code/model/tapas_physio_rvt_hilbert.m +++ b/PhysIO/code/model/tapas_physio_rvt_hilbert.m @@ -142,42 +142,21 @@ %% Plot figures %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% save relevant structures + +verbose.review.rvt_hilbert.t = t; +verbose.review.rvt_hilbert.fr = fr; +verbose.review.rvt_hilbert.fr_lp = fr_lp; +verbose.review.rvt_hilbert.fr_mag = fr_mag; +verbose.review.rvt_hilbert.fr_rv = fr_rv +verbose.review.rvt_hilbert.fr_phase = fr_phase; +verbose.review.rvt_hilbert.fr_if = fr_if; + if verbose.level>=2 - verbose.fig_handles(end+1) = tapas_physio_get_default_fig_params(); - set(gcf, 'Name', 'Model: Hilbert RVT (Respiratory Volume per Time)'); + [verbose] = tapas_physio_plot_rvt_hilbert(t,fr, fr_lp, fr_mag, fr_rv, ... + fr_phase, fr_if, verbose) + - hs(1) = subplot(2,1,1); - hold on; - plot([t(1), t(end)], [0.0, 0.0], 'Color', [0.7, 0.7, 0.7]); - hp(1) = plot(t, fr); - hp(2) = plot(t, fr_lp); - hp(3) = plot(t, fr_mag); - hp(4) = plot(t, fr_rv / 2.0); - strLegend = { - 'Filtered breathing signal', ... - '... after low pass-filter', ... - 'Breathing signal envelope', ... - 'Respiratory volume'}; - legend(hp, strLegend); - xlim([t(1), t(end)]); - title('Respiratory Volume (from Hilbert Transform)'); - - hs(2) = subplot(2,1,2); - hold on - plot([t(1), t(end)], [0.0, 0.0], 'Color', [0.7, 0.7, 0.7]); - hp(1) = plot(t, fr); - hp(2) = plot(t, fr_lp); - hp(3) = plot(t, std(fr) * cos(fr_phase)); - hp(4) = plot(t, fr_if); - strLegend = { - 'Filtered breathing signal', ... - '... after low pass-filter', ... - '... after removing amplitude', ... - 'Instantaneous breathing rate'}; - legend(hp, strLegend); - xlim([t(1), t(end)]); - title('Instantaneous Breathing Rate (from Hilbert Transform)'); - linkaxes(hs, 'x') end %% Downsample %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/PhysIO/code/plot/tapas_physio_plot_create_hrv_regressors.m b/PhysIO/code/plot/tapas_physio_plot_create_hrv_regressors.m new file mode 100644 index 00000000..b551a665 --- /dev/null +++ b/PhysIO/code/plot/tapas_physio_plot_create_hrv_regressors.m @@ -0,0 +1,30 @@ +function [verbose] = tapas_physio_plot_create_hrv_regressors(sample_points, hrOut,... + hr, t, crf, convHRV, delays,samplePointsOut, convHRVOut, verbose) +%UNTITLED2 Summary of this function goes here +% Detailed explanation goes here + + + verbose.fig_handles(end+1) = tapas_physio_get_default_fig_params(); + set(gcf, 'Name', 'Model: 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)'); + + subplot(2,2,2) + plot(t, crf,'r');xlabel('time (seconds)'); + title('Cardiac response function'); + + subplot(2,2,3) + plot(sample_points, convHRV,'r');xlabel('time (seconds)'); + title('Heart rate X cardiac response function'); + + 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 \ No newline at end of file diff --git a/PhysIO/code/plot/tapas_physio_plot_create_scan_timing_philips.m b/PhysIO/code/plot/tapas_physio_plot_create_scan_timing_philips.m new file mode 100644 index 00000000..14a483ee --- /dev/null +++ b/PhysIO/code/plot/tapas_physio_plot_create_scan_timing_philips.m @@ -0,0 +1,22 @@ +function [verbose] = tapas_physio_plot_create_scan_timing_philips(t, y, ... + VOLLOCS, Ndummies, LOCS, LOC_END_MARKER) +%UNTITLED Summary of this function goes here +% Detailed explanation goes here + +verbose.fig_handles(end+1) = tapas_physio_get_default_fig_params(); +set(gcf,'Name', 'Sync: 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 \ No newline at end of file diff --git a/PhysIO/code/plot/tapas_physio_plot_filter_respiratory.m b/PhysIO/code/plot/tapas_physio_plot_filter_respiratory.m new file mode 100644 index 00000000..a5ec9f8a --- /dev/null +++ b/PhysIO/code/plot/tapas_physio_plot_filter_respiratory.m @@ -0,0 +1,33 @@ +function [verbose, handles, labels] = tapas_physio_plot_filter_respiratory(rpulset,m, s, t, ... + rpulset_out, rpulset_out_trend, trend, rpulset_out_trend_filt, verbose) +% plot repiratory filtering +% (C) 2023 Johanna Bayer + verbose.fig_handles(end+1) = tapas_physio_get_default_fig_params(); + set(gcf, 'Name', 'Preproc: Respiratory filtering'); + hold on; + handles = []; labels = {}; + plot([t(1), t(end)], [0.0, 0.0], 'Color', [0.7, 0.7, 0.7]); + handles(end+1) = plot(t, (rpulset - m) / s); + labels{end+1} = 'Raw respiratory signal'; + + % without outliers + handles(end+1) = plot(t, (rpulset_out - m) / s); + labels{end+1} = '... without outliers'; + + % detrend + figure(verbose.fig_handles(end)); + handles(end+1) = plot(t, (trend - m) / s); + labels{end+1} = '... low frequency trend'; + plot([t(1), t(end)], [-5.0, -5.0], 'Color', [0.7, 0.7, 0.7]); + handles(end+1) = plot(t, (rpulset_out_trend - m) / s - 5.0); + labels{end+1} = '... detrended'; + + % Low pass filtered to remove noise + handles(end+1) = plot(t, (rpulset_out_trend_filt - m) / s - 5.0); + labels{end+1} = '... after low-pass filter'; + xlim([t(1), t(end)]); + xlabel('Time (s)'); + yticks([]); + legend(handles, labels); + +end \ No newline at end of file diff --git a/PhysIO/code/plot/tapas_physio_plot_get_cardiac_phase.m b/PhysIO/code/plot/tapas_physio_plot_get_cardiac_phase.m new file mode 100644 index 00000000..b624c36f --- /dev/null +++ b/PhysIO/code/plot/tapas_physio_plot_get_cardiac_phase.m @@ -0,0 +1,25 @@ +function [verbose] = tapas_physio_plot_get_cardiac_phase(scannert,cardiac_phase, pulset, ... + svolpulse, verbose) + + % 1. plot chosen slice start event + % 2. plot chosen c_sample phase on top of chosen slice scan start, (as a stem + % and line of phases) + % 3. plot all detected cardiac r-wave peaks + % 4. plot volume start event + stringTitle = 'Preproc: tapas_physio_get_cardiac_phase: scanner and R-wave pulses - output phase'; + fh = tapas_physio_get_default_fig_params(); + set(fh, 'Name', stringTitle); + stem(scannert, cardiac_phase, 'k'); hold on; + plot(scannert, cardiac_phase, 'k'); + stem(pulset,3*ones(size(pulset)),'r', 'LineWidth',2); + stem(svolpulse,7*ones(size(svolpulse)),'g', 'LineWidth',2); + legend('estimated phase at slice events', ... + '', ... + 'heart beat R-peak', ... + 'scan volume start'); + title(regexprep(stringTitle,'_', '\\_')); + xlabel('t (seconds)'); + %stem(scannertpriorpulse,ones(size(scannertpriorpulse))*2,'g'); + %stem(scannertafterpulse,ones(size(scannertafterpulse))*2,'b'); + +end \ No newline at end of file diff --git a/PhysIO/code/plot/tapas_physio_plot_iterative_template_creation.m b/PhysIO/code/plot/tapas_physio_plot_iterative_template_creation.m new file mode 100644 index 00000000..7350c575 --- /dev/null +++ b/PhysIO/code/plot/tapas_physio_plot_iterative_template_creation.m @@ -0,0 +1,57 @@ +function [verbose] = tapas_physio_plot_iterative_template_creation(hasFirstGuessPeaks,... + t, c, cpulse1stGuess, nPulses1, nPulses2, cpulse2ndGuess, meanLag1, meanLag2, verbose) +% Plot refined heartrate templates +% (C) 2023 Johanna Bayer + +stringTitle = 'Preproc: Iterative Template Creation Single Cycle'; + +if hasFirstGuessPeaks + + fh = tapas_physio_get_default_fig_params(); + set(fh, 'Name', stringTitle); + verbose.fig_handles(end+1) = fh; + subplot(3,1,1); + hold off + hp(3) = plot(t, c, 'k'); + hold all; + hp(1) = stem(t(cpulse1stGuess), ... + 4*ones(nPulses1,1),'b'); + + hp(2) = stem(t(cpulse2ndGuess),... + 4*ones(nPulses2,1),'r'); + + stringLegend = { + sprintf('1st guess peaks (N =%d)', nPulses1), ... + sprintf('2nd guess peaks (N =%d)', nPulses2), ... + 'raw time series' + }; + + legend(hp, stringLegend); + title('Finding first peak (cycle start), backwards') + + + + %% Plot difference between detected events + subplot(3,1,2); + + + plot(t(cpulse1stGuess(2:end)), diff(t(cpulse1stGuess))); + hold all + plot(t(cpulse2ndGuess(2:end)), diff(t(cpulse2ndGuess))); + title('Temporal lag between events') + + stringLegend = { + sprintf('1st Guess (Mean lag duration %3.1f s)', meanLag1), ... + sprintf('2nd Guess (Mean lag duration %3.1f s)', meanLag2) ... + }; + + legend(stringLegend); + +else + fh = tapas_physio_get_default_fig_params(); + verbose.fig_handles(end+1) = fh; + subplot(3,1,1); + plot(t, c, 'k'); title('Preproc: Finding first peak of cycle, backwards') + +end +end \ No newline at end of file diff --git a/PhysIO/code/plot/tapas_physio_plot_peak_detection_from_automatically_generated.m b/PhysIO/code/plot/tapas_physio_plot_peak_detection_from_automatically_generated.m new file mode 100644 index 00000000..3e5fcca8 --- /dev/null +++ b/PhysIO/code/plot/tapas_physio_plot_peak_detection_from_automatically_generated.m @@ -0,0 +1,19 @@ +function [verbose] = tapas_physio_plot_peak_detection_from_automatically_generated(t,c,... + cpulse, verbose) +% Helper function to create peak detection plot. +% IN: t, c, cpulse +% Author Johanna Bayer, (C) 2023 + + verbose.fig_handles(end+1) = tapas_physio_get_default_fig_params(); + stringTitle = 'Preproc: Peak Detection from Automatically Generated Template'; + set(gcf, 'Name', stringTitle); + plot(t, c, 'k'); + hold all; + stem(cpulse,4*ones(size(cpulse)), 'r'); + legend('Raw time course',... + 'Detected maxima (cardiac pulses / max inhalations)'); + title(stringTitle); + + +end + diff --git a/PhysIO/code/plot/tapas_physio_plot_raw_physdata_siemens.m b/PhysIO/code/plot/tapas_physio_plot_raw_physdata_siemens.m index 564d6299..3df8c2b1 100644 --- a/PhysIO/code/plot/tapas_physio_plot_raw_physdata_siemens.m +++ b/PhysIO/code/plot/tapas_physio_plot_raw_physdata_siemens.m @@ -1,17 +1,16 @@ -function fh = tapas_physio_plot_raw_physdata_siemens(dataCardiac) +function fh = tapas_physio_plot_raw_physdata_siemens(dataPhysio) % plots cardiac data as extracted from Siemens log file % -% output = tapas_physio_plot_raw_physdata_siemens(input) +% fh = tapas_physio_plot_raw_physdata_siemens(dataCardiac) % % IN -% +% dataPhysio output struct from tapas_physio_siemens_table2cardiac % OUT % % EXAMPLE % tapas_physio_plot_raw_physdata_siemens % -% See also tapas_physio_read_physlogfiles_siemens -% See also tapas_physio_siemens_table2cardiac +% See also tapas_physio_read_physlogfiles_siemens tapas_physio_siemens_table2cardiac % Author: Lars Kasper % Created: 2016-02-29 @@ -23,7 +22,7 @@ % (either version 3 or, at your option, any later version). For further details, see the file % COPYING or . -tapas_physio_strip_fields(dataCardiac); +tapas_physio_strip_fields(dataPhysio); stringTitle = 'Read-In: Raw Siemens physlog data'; fh = tapas_physio_get_default_fig_params(); @@ -31,22 +30,28 @@ stem(cpulse_on, ampl*ones(size(cpulse_on)), 'g'); hold all; stem(cpulse_off, ampl*ones(size(cpulse_off)), 'r'); stem(t(stopSample), ampl , 'm'); -plot(t, channel_1); -plot(t, channel_AVF); -plot(t, meanChannel); -stringLegend = { ... - 'cpulse on', 'cpulse off', 'assumed last sample of last scan volume', ... - 'channel_1', 'channel_{AVF}', 'mean of channels'}; +stringLegend = {'physio trigger on', 'physio trigger off', ... + 'assumed last sample of last scan volume'}; + +nChannels = size(recordingChannels, 2); + +for iChannel = 1:nChannels + plot(t, recordingChannels(:,iChannel)); + stringLegend{end+1} = sprintf('channel %d', iChannel); +end + +plot(t, meanChannel); +stringLegend{end+1} = 'mean of channels'; if ~isempty(recording_on) stem(recording_on, ampl*ones(size(recording_on)), 'k'); - stringLegend{end+1} = 'phys recording on'; + stringLegend{end+1} = 'physio recording on'; end if ~isempty(recording_off) stem(recording_off, ampl*ones(size(recording_off)), 'k'); - stringLegend{end+1} = 'phys recording off'; + stringLegend{end+1} = 'physio recording off'; end legend(stringLegend); title(stringTitle); diff --git a/PhysIO/code/plot/tapas_physio_plot_rvt_hilbert.m b/PhysIO/code/plot/tapas_physio_plot_rvt_hilbert.m new file mode 100644 index 00000000..35826295 --- /dev/null +++ b/PhysIO/code/plot/tapas_physio_plot_rvt_hilbert.m @@ -0,0 +1,40 @@ +function [verbose] = tapas_physio_plot_rvt_hilbert(t,fr, fr_lp, fr_mag, fr_rv, ... + fr_phase, fr_if, verbose); +%UNTITLED7 Summary of this function goes here +% Detailed explanation goes here + verbose.fig_handles(end+1) = tapas_physio_get_default_fig_params(); + set(gcf, 'Name', 'Model: Hilbert RVT (Respiratory Volume per Time)'); + + hs(1) = subplot(2,1,1); + hold on; + plot([t(1), t(end)], [0.0, 0.0], 'Color', [0.7, 0.7, 0.7]); + hp(1) = plot(t, fr); + hp(2) = plot(t, fr_lp); + hp(3) = plot(t, fr_mag); + hp(4) = plot(t, fr_rv / 2.0); + strLegend = { + 'Filtered breathing signal', ... + '... after low pass-filter', ... + 'Breathing signal envelope', ... + 'Respiratory volume'}; + legend(hp, strLegend); + xlim([t(1), t(end)]); + title('Respiratory Volume (from Hilbert Transform)'); + + hs(2) = subplot(2,1,2); + hold on + plot([t(1), t(end)], [0.0, 0.0], 'Color', [0.7, 0.7, 0.7]); + hp(1) = plot(t, fr); + hp(2) = plot(t, fr_lp); + hp(3) = plot(t, std(fr) * cos(fr_phase)); + hp(4) = plot(t, fr_if); + strLegend = { + 'Filtered breathing signal', ... + '... after low pass-filter', ... + '... after removing amplitude', ... + 'Instantaneous breathing rate'}; + legend(hp, strLegend); + xlim([t(1), t(end)]); + title('Instantaneous Breathing Rate (from Hilbert Transform)'); + linkaxes(hs, 'x') +end \ No newline at end of file diff --git a/PhysIO/code/plot/tapas_physio_plot_sync_bundles.m b/PhysIO/code/plot/tapas_physio_plot_sync_bundles.m new file mode 100644 index 00000000..f3c0d1ed --- /dev/null +++ b/PhysIO/code/plot/tapas_physio_plot_sync_bundles.m @@ -0,0 +1,10 @@ +function [verbose] = tapas_physio_plot_sync_bundles(Nallvols, t, SLICELOCS, verbose) +%UNTITLED Summary of this function goes here +% Detailed explanation goes here + stringTitle = 'Sync: Slice bundles belonging to 1 volume'; + verbose.fig_handles(end+1) = tapas_physio_get_default_fig_params(); + set(gcf, 'Name', stringTitle); + for v=1:Nallvols-1, stem(t(SLICELOCS{v}),ones(size(SLICELOCS{v})));hold all;end + title(stringTitle); + xlabel('t (seconds since SCANPHYSLOG-start)'); +end \ No newline at end of file diff --git a/PhysIO/code/plot/tapas_physio_plot_templates_of_cycle_time.m b/PhysIO/code/plot/tapas_physio_plot_templates_of_cycle_time.m new file mode 100644 index 00000000..91220271 --- /dev/null +++ b/PhysIO/code/plot/tapas_physio_plot_templates_of_cycle_time.m @@ -0,0 +1,25 @@ +function [verbose] = tapas_physio_plot_templates_of_cycle_time(tTemplate, ... + template, pulseTemplate, pulseCleanedTemplate, verbose) +% (c) 2023 Johanna Bayer + + % First plot + fh = verbose.fig_handles(end); + figure(fh); + subplot(3,1,3); + + plot(tTemplate, template'); + hold all; + hp(1) = plot(tTemplate, pulseTemplate', '.--g', 'LineWidth', 3, 'Marker', ... + 'o'); + xlabel('t (seconds)'); + title('Templates of cycle time course and mean template'); + + % Second plot + stringTitle = 'Preproc: Iterative Template Creation Single Cycle'; + hp(2) = plot(tTemplate, pulseCleanedTemplate, '.-g', 'LineWidth', 4, ... + 'Marker', 'x'); + legend(hp, 'mean of templates', 'mean of most similar, chosen templates'); + tapas_physio_suptitle(stringTitle); + + +end \ No newline at end of file diff --git a/PhysIO/code/plot/tapas_physio_plot_traces.m b/PhysIO/code/plot/tapas_physio_plot_traces.m new file mode 100644 index 00000000..dae53188 --- /dev/null +++ b/PhysIO/code/plot/tapas_physio_plot_traces.m @@ -0,0 +1,54 @@ +function fh = tapas_physio_plot_traces(pulset, rsampint, rout, resp_max, ... + cumsumh, sumh, h, npulse, dpulse, rphase) + +nsamples = length(pulset); +t = (0:nsamples-1)*rsampint; +feqht = cumsumh/sumh*pi; + +fh = tapas_physio_get_default_fig_params(); +set(fh, 'Name', ... + 'Preproc: get_respiratory_phase: histogram for respiratory phase estimation'); + +hs(1) = subplot(2,2,1); +plot(t,pulset); +xlabel('t (s)'); +ylabel('breathing amplitude (a. u.)'); +title('(filtered) breathing time series'); + +if resp_max < inf + hold on; + plot(t, ones(size(t)) * resp_max, 'k--'); + hold on; + hp = plot(t, -ones(size(t)) * resp_max, 'k--'); + legend(hp, ... + 'threshold for maximum amplitude to be considered in histogram'); + set(gcf, 'Name', ... + [get(gcf, 'Name') ' - with amplitude overshoot-correction']); +end + +hs(2) = subplot(2,2,2); +bar(rout, h); +xlabel('normalized breathing amplitude'); +ylabel('counts'); +title('histogram for phase mapping'); +xlim([-0.1 1.1]); + +hs(3) = subplot(2,2,3); plot(rout, [feqht, cos(feqht), sin(feqht)]); +xlabel('normalized breathing amplitude'); +title(... + 'equalized histogram bin amplitude -> phase transfer function (f_{eqht})'); +legend('f: normalized amplitude -> phase transfer function', 'cos(f)', ... + 'sin(f)', 'Location', 'NorthWest'); + +%figure('Name', 'Histogram: Respiration phase estimation'); +hs(4) = subplot(2,2,4); +plot(t, [npulse*10, dpulse, (rphase-pi)]); +legend('10*normalized breathing belt amplitude', ... + '-1 = exhale, 1 = inhale', 'estimated respiratory phase'); +ylim([-0.2 10.2]); +title('Histogram-based respiration phase estimation'); + +linkaxes(hs([1 4]), 'x'); +linkaxes(hs([2 3]), 'x'); + +end \ No newline at end of file diff --git a/PhysIO/code/preproc/tapas_physio_filter_respiratory.m b/PhysIO/code/preproc/tapas_physio_filter_respiratory.m index 6c3d0305..2a85ae3e 100644 --- a/PhysIO/code/preproc/tapas_physio_filter_respiratory.m +++ b/PhysIO/code/preproc/tapas_physio_filter_respiratory.m @@ -61,23 +61,19 @@ %% Basic preproc and outlier removal -% If rpulset has nans, replace them with zeros -rpulsetOffset = nanmean(rpulset); -rpulset(isnan(rpulset)) = nanmean(rpulset); +% If rpulset has nans, replace them with mean of valid values +try + rpulsetOffset = mean(rpulset, 'omitnan'); +catch % for backwards compatibility < Matlab 2016a + rpulsetOffset = nanmean(rpulset); +end +rpulset(isnan(rpulset)) = rpulsetOffset; rpulset = detrend(rpulset, 3); % Demean / detrend to reduce edge effects -if verbose.level>=3 - verbose.fig_handles(end+1) = tapas_physio_get_default_fig_params(); - set(gcf, 'Name', 'Preproc: Respiratory filtering'); - hold on; - handles = []; labels = {}; - t = linspace(0.0, rsampint * (length(rpulset) - 1), length(rpulset)); - plot([t(1), t(end)], [0.0, 0.0], 'Color', [0.7, 0.7, 0.7]); - m = mean(rpulset); s = std(rpulset); - handles(end+1) = plot(t, (rpulset - m) / s); - labels{end+1} = 'Raw respiratory signal'; -end + m = mean(rpulset); s = std(rpulset); + t = linspace(0.0, rsampint * (length(rpulset) - 1), length(rpulset)); + % Now do a check for any gross outliers relative to the statistics of the % whole timeseries @@ -86,9 +82,10 @@ mpulse = mean(rpulset); stdpulse = 1.4826 * mad(rpulset, 1); % Robust to outliers: https://en.wikipedia.org/wiki/Robust_measures_of_scale outliers = (rpulset > (mpulse + (z_thresh * stdpulse))); -rpulset(outliers) = mpulse + (z_thresh * stdpulse); -outliers = (rpulset < (mpulse - (z_thresh * stdpulse))); -rpulset(outliers) = mpulse - (z_thresh * stdpulse); +rpulset_out = rpulset; +rpulset_out(outliers) = mpulse + (z_thresh * stdpulse); +outliers = (rpulset_out < (mpulse - (z_thresh * stdpulse))); +rpulset_out(outliers) = mpulse - (z_thresh * stdpulse); % if verbose.level>=3 % plot([t(1), t(end)], [z_thresh, z_thresh], 'Color', [0.7, 0.7, 0.7]); % plot([t(1), t(end)], [-z_thresh, -z_thresh], 'Color', [0.7, 0.7, 0.7]); @@ -98,11 +95,11 @@ if despike mad_thresh = 5.0; % Again, relatively high so only get large spikes (low-pass filter gets the rest) n_pad = ceil(0.25 / rsampint); % 0.5 s total window length - rpulset_padded = padarray(rpulset, n_pad, 'symmetric'); + rpulset_padded = padarray(rpulset_out, n_pad, 'symmetric'); medians = movmedian(rpulset_padded, 2 * n_pad + 1, 'Endpoints', 'discard'); mads = movmad(rpulset_padded, 2 * n_pad + 1, 'Endpoints', 'discard'); - outliers = (abs(rpulset - medians) > mad_thresh * mads); - rpulset(outliers) = medians(outliers); + outliers = (abs(rpulset_out - medians) > mad_thresh * mads); + rpulset_out(outliers) = medians(outliers); % if verbose.level>=3 % plot(t, (medians - m) / s, 'Color', [0.7, 0.7, 0.7]); % plot(t, (medians + mad_thresh * mads - m) / s, 'Color', [0.7, 0.7, 0.7]); @@ -110,10 +107,6 @@ % end end -if verbose.level>=3 - handles(end+1) = plot(t, (rpulset - m) / s); - labels{end+1} = '... without outliers'; -end %% Detrend and remove noise via filtering @@ -130,7 +123,7 @@ % Use a large padding, and window so tapers back to mean naturally padding_window = window(@blackmanharris, 2 * n_pad + 1); -rpulset_padded = padarray(rpulset, n_pad, 'symmetric'); +rpulset_padded = padarray(rpulset_out, n_pad, 'symmetric'); rpulset_padded(1:n_pad) = padding_window(1:n_pad) .* rpulset_padded(1:n_pad); rpulset_padded(end-n_pad+1:end) = padding_window(end-n_pad+1:end) .* rpulset_padded(end-n_pad+1:end); @@ -140,42 +133,42 @@ % figure('Name', 'rpulset_padded'); plot(rpulset_padded); hold on; plot(trend) trend = trend(n_pad+1:end-n_pad); -rpulset = rpulset - trend; +rpulset_out_trend = rpulset_out - trend; -if verbose.level>=3 - figure(verbose.fig_handles(end)); - handles(end+1) = plot(t, (trend - m) / s); - labels{end+1} = '... low frequency trend'; - plot([t(1), t(end)], [-5.0, -5.0], 'Color', [0.7, 0.7, 0.7]); - handles(end+1) = plot(t, (rpulset - m) / s - 5.0); - labels{end+1} = '... detrended'; -end % Low-pass filter to remove noise d = designfilt( ... 'lowpassiir', 'HalfPowerFrequency', cutoff_freqs(2), ... 'FilterOrder', 20, 'SampleRate', sampfreq); -rpulset = filtfilt(d, padarray(rpulset, n_pad, 'symmetric')); -rpulset = rpulset(n_pad+1:end-n_pad); +rpulset_out_trend_filt = filtfilt(d, padarray(rpulset_out_trend, n_pad, 'symmetric')); +rpulset_out_trend_filt = rpulset_out_trend_filt(n_pad+1:end-n_pad); -if verbose.level>=3 - handles(end+1) = plot(t, (rpulset - m) / s - 5.0); - labels{end+1} = '... after low-pass filter'; -end %% Normalise, if requested if normalize - rpulset = rpulset/max(abs(rpulset)); + rpulset_out_trend_filt = rpulset_out_trend_filt/max(abs(rpulset_out_trend_filt)); end %% + % save relevant variavles for retrospective plotting + verbose.review.resp_filter.rpulset = rpulset; + verbose.review.resp_filter.rsampint = rsampint; + verbose.review.resp_filter.m = m; + verbose.review.resp_filter.s = s; + verbose.review.resp_filter.t = t; + verbose.review.resp_filter.rpulset_out = rpulset_out; + verbose.review.resp_filter.rpulset_out_trend = rpulset_out_trend; + verbose.review.resp_filter.trend = trend; + verbose.review.resp_filter.rpulset_out_trend_filt = rpulset_out_trend_filt; + +% Debug and plot if verbose.level>=3 - xlim([t(1), t(end)]); - xlabel('Time (s)'); - yticks([]); - legend(handles, labels); + [verbose] = tapas_physio_plot_filter_respiratory(rpulset,m, s, t, ... + rpulset_out, rpulset_out_trend,trend,rpulset_out_trend_filt, verbose); + end +rpulset = rpulset_out_trend_filt; end \ No newline at end of file diff --git a/PhysIO/code/preproc/tapas_physio_get_cardiac_phase.m b/PhysIO/code/preproc/tapas_physio_get_cardiac_phase.m index 6d71ee61..b70fa4d6 100644 --- a/PhysIO/code/preproc/tapas_physio_get_cardiac_phase.m +++ b/PhysIO/code/preproc/tapas_physio_get_cardiac_phase.m @@ -80,24 +80,13 @@ cardiac_phase=(2*pi*(scannert'-scannertpriorpulse)./(scannertafterpulse-scannertpriorpulse))'; if isVerbose - % 1. plot chosen slice start event - % 2. plot chosen c_sample phase on top of chosen slice scan start, (as a stem - % and line of phases) - % 3. plot all detected cardiac r-wave peaks - % 4. plot volume start event - stringTitle = 'Preproc: tapas_physio_get_cardiac_phase: scanner and R-wave pulses - output phase'; - fh = tapas_physio_get_default_fig_params(); - set(fh, 'Name', stringTitle); - stem(scannert, cardiac_phase, 'k'); hold on; - plot(scannert, cardiac_phase, 'k'); - stem(pulset,3*ones(size(pulset)),'r', 'LineWidth',2); - stem(svolpulse,7*ones(size(svolpulse)),'g', 'LineWidth',2); - legend('estimated phase at slice events', ... - '', ... - 'heart beat R-peak', ... - 'scan volume start'); - title(regexprep(stringTitle,'_', '\\_')); - xlabel('t (seconds)'); - %stem(scannertpriorpulse,ones(size(scannertpriorpulse))*2,'g'); - %stem(scannertafterpulse,ones(size(scannertafterpulse))*2,'b'); + [verbose] = tapas_physio_plot_get_cardiac_phase(scannert,cardiac_phase, pulset, ... + svolpulse, verbose) + + % save relevent stuctures + verbose.review.get_cardiac.scannert = scannert; + verbose.review.get_cardiac.cardiac_phase = cardiac_phase; + verbose.review.get_cardiac.pulset = pulset; + verbose.review.get_cardiac.svolpulse = svolpulse; + end \ No newline at end of file diff --git a/PhysIO/code/preproc/tapas_physio_get_cardiac_pulse_template.m b/PhysIO/code/preproc/tapas_physio_get_cardiac_pulse_template.m index fbdb10f8..96f924fc 100644 --- a/PhysIO/code/preproc/tapas_physio_get_cardiac_pulse_template.m +++ b/PhysIO/code/preproc/tapas_physio_get_cardiac_pulse_template.m @@ -1,4 +1,4 @@ -function [pulseCleanedTemplate, cpulse2ndGuess, averageHeartRateInSamples] = ... +function [pulseCleanedTemplate, cpulse2ndGuess, averageHeartRateInSamples, verbose] = ... tapas_physio_get_cardiac_pulse_template(t, c, verbose, ... varargin) % determines cardiac template by a 2-pass peak detection and averaging of @@ -52,76 +52,50 @@ hasFirstGuessPeaks = ~isempty(cpulse1stGuess); -%% Second step, refined heart rate estimate - -stringTitle = 'Preproc: Iterative Template Creation Single Cycle'; - if hasFirstGuessPeaks - + averageHeartRateInSamples = round(mean(diff(cpulse1stGuess))); [tmp, cpulse2ndGuess] = tapas_physio_findpeaks(c,... 'minpeakheight', thresh_min,... 'minpeakdistance', round(shortenTemplateFactor*... averageHeartRateInSamples)); - + if doDebug - + +%% Second step, refined heart rate estimate + nPulses1 = length(cpulse1stGuess); nPulses2 = length(cpulse2ndGuess); - fh = tapas_physio_get_default_fig_params(); - set(fh, 'Name', stringTitle); - verbose.fig_handles(end+1) = fh; - subplot(3,1,1); - hold off - hp(3) = plot(t, c, 'k'); - hold all; - hp(1) = stem(t(cpulse1stGuess), ... - 4*ones(nPulses1,1),'b'); - - hp(2) = stem(t(cpulse2ndGuess),... - 4*ones(nPulses2,1),'r'); - - stringLegend = { - sprintf('1st guess peaks (N =%d)', nPulses1), ... - sprintf('2nd guess peaks (N =%d)', nPulses2), ... - 'raw time series' - }; - - legend(hp, stringLegend); - title('Finding first peak (cycle start), backwards') - - - - %% Plot difference between detected events - subplot(3,1,2); - + meanLag1 = mean(diff(t(cpulse1stGuess))); meanLag2 = mean(diff(t(cpulse2ndGuess))); - - plot(t(cpulse1stGuess(2:end)), diff(t(cpulse1stGuess))); - hold all - plot(t(cpulse2ndGuess(2:end)), diff(t(cpulse2ndGuess))); - title('Temporal lag between events') - - stringLegend = { - sprintf('1st Guess (Mean lag duration %3.1f s)', meanLag1), ... - sprintf('2nd Guess (Mean lag duration %3.1f s)', meanLag2) ... - }; - - legend(stringLegend); - end -else - if doDebug - fh = tapas_physio_get_default_fig_params(); - verbose.fig_handles(end+1) = fh; - subplot(3,1,1); - plot(t, c, 'k'); title('Preproc: Finding first peak of cycle, backwards') + end + +else verbose = tapas_physio_log(['No peaks found in raw cardiac time series. Check raw ' ... 'physiological recordings figure whether there is any non-constant ' ... 'cardiac data'], verbose, 2); % error! -end +end + +%% Plot in case of Debugging (verbose =>3) + +if doDebug + + [verbose] = tapas_physio_plot_iterative_template_creation(hasFirstGuessPeaks,... + t, c, cpulse1stGuess, nPulses1, nPulses2, cpulse2ndGuess, meanLag1, meanLag2, verbose); + %save relevant functions + verbose.review.iter_temp.hasFirstGuessPeaks = hasFirstGuessPeaks; + verbose.review.iter_temp.t = t; + verbose.review.iter_temp.c = c; + verbose.review.iter_temp.cpulse1stGuess = cpulse1stGuess; + verbose.review.iter_temp.nPulses1 = nPulses1; + verbose.review.iter_temp.nPulses2 = nPulses2; + verbose.review.iter_temp.cpulse2ndGuess = cpulse2ndGuess; + verbose.review.iter_temp.meanLag1 = meanLag1; + verbose.review.iter_temp.meanLag2 = meanLag2; +end %% Build template based on the guessed peaks: @@ -133,7 +107,7 @@ (averageHeartRateInSamples / 2)); -[pulseCleanedTemplate, pulseTemplate] = tapas_physio_get_template_from_pulses(... +[pulseCleanedTemplate, pulseTemplate, verbose] = tapas_physio_get_template_from_pulses(... c, cpulse2ndGuess, halfTemplateWidthInSamples, ... verbose, 'doNormalizeTemplate', doNormalizeTemplate, ... 'thresholdHighQualityCorrelation', thresholdHighQualityCorrelation, ... diff --git a/PhysIO/code/preproc/tapas_physio_get_cardiac_pulses_auto_matched.m b/PhysIO/code/preproc/tapas_physio_get_cardiac_pulses_auto_matched.m index d9c2c1a2..b33fcd4b 100644 --- a/PhysIO/code/preproc/tapas_physio_get_cardiac_pulses_auto_matched.m +++ b/PhysIO/code/preproc/tapas_physio_get_cardiac_pulses_auto_matched.m @@ -51,7 +51,7 @@ c = c-mean(c); c = c./std(c); % normalize time series %% Determine template for QRS-wave (or pulse) -[pulseCleanedTemplate, cpulseSecondGuess, averageHeartRateInSamples] = ... +[pulseCleanedTemplate, cpulseSecondGuess, averageHeartRateInSamples, verbose] = ... tapas_physio_get_cardiac_pulse_template(t, c, verbose, ... 'thresh_min', thresh_min, ... 'minCycleSamples', minPulseDistanceSamples, ... @@ -99,15 +99,14 @@ cpulse = t(cpulse); if verbose.level >=2 - verbose.fig_handles(end+1) = tapas_physio_get_default_fig_params(); - stringTitle = 'Preproc: Peak Detection from Automatically Generated Template'; - set(gcf, 'Name', stringTitle); - plot(t, c, 'k'); - hold all; - stem(cpulse,4*ones(size(cpulse)), 'r'); - legend('Raw time course',... - 'Detected maxima (cardiac pulses / max inhalations)'); - title(stringTitle); + + [verbose] = tapas_physio_plot_peak_detection_from_automatically_generated(t, c, ... + cpulse, verbose); + + % freeze versions of t, c and cpulse for review + verbose.review.peak.t = t; + verbose.review.peak.c = c; + verbose.review.peak.cpulse = cpulse; end end diff --git a/PhysIO/code/preproc/tapas_physio_get_respiratory_phase.m b/PhysIO/code/preproc/tapas_physio_get_respiratory_phase.m index c0019dd1..1260ea4c 100644 --- a/PhysIO/code/preproc/tapas_physio_get_respiratory_phase.m +++ b/PhysIO/code/preproc/tapas_physio_get_respiratory_phase.m @@ -1,4 +1,5 @@ -function [rphase, fh] = tapas_physio_get_respiratory_phase(pulset, rsampint,... +function [pulset, rsampint, rout, resp_max, cumsumh, sumh, h, ... + npulse, dpulse, rphase, fh] = tapas_physio_get_respiratory_phase(pulset, rsampint,... verbose, thresh) % get_respiratory_phase is a function for creating respiratory phase regressor. @@ -102,7 +103,7 @@ rphase = pi*(cumsumh(binnum)/sumh).*dpulse+pi; if verbose - fh = plot_traces(pulset, rsampint, rout, resp_max, cumsumh, sumh, h, ... + fh = tapas_physio_plot_traces(pulset, rsampint, rout, resp_max, cumsumh, sumh, h, ... npulse, dpulse, rphase); else fh = []; @@ -110,57 +111,4 @@ end -function fh = plot_traces(pulset, rsampint, rout, resp_max, ... - cumsumh, sumh, h, npulse, dpulse, rphase) - -nsamples = length(pulset); -t = (0:nsamples-1)*rsampint; -feqht = cumsumh/sumh*pi; - -fh = tapas_physio_get_default_fig_params(); -set(fh, 'Name', ... - 'Preproc: get_respiratory_phase: histogram for respiratory phase estimation'); - -hs(1) = subplot(2,2,1); -plot(t,pulset); -xlabel('t (s)'); -ylabel('breathing amplitude (a. u.)'); -title('(filtered) breathing time series'); - -if resp_max < inf - hold on; - plot(t, ones(size(t)) * resp_max, 'k--'); - hold on; - hp = plot(t, -ones(size(t)) * resp_max, 'k--'); - legend(hp, ... - 'threshold for maximum amplitude to be considered in histogram'); - set(gcf, 'Name', ... - [get(gcf, 'Name') ' - with amplitude overshoot-correction']); -end - -hs(2) = subplot(2,2,2); -bar(rout, h); -xlabel('normalized breathing amplitude'); -ylabel('counts'); -title('histogram for phase mapping'); -xlim([-0.1 1.1]); - -hs(3) = subplot(2,2,3); plot(rout, [feqht, cos(feqht), sin(feqht)]); -xlabel('normalized breathing amplitude'); -title(... - 'equalized histogram bin amplitude -> phase transfer function (f_{eqht})'); -legend('f: normalized amplitude -> phase transfer function', 'cos(f)', ... - 'sin(f)', 'Location', 'NorthWest'); - -%figure('Name', 'Histogram: Respiration phase estimation'); -hs(4) = subplot(2,2,4); -plot(t, [npulse*10, dpulse, (rphase-pi)]); -legend('10*normalized breathing belt amplitude', ... - '-1 = exhale, 1 = inhale', 'estimated respiratory phase'); -ylim([-0.2 10.2]); -title('Histogram-based respiration phase estimation'); - -linkaxes(hs([1 4]), 'x'); -linkaxes(hs([2 3]), 'x'); -end diff --git a/PhysIO/code/preproc/tapas_physio_get_template_from_pulses.m b/PhysIO/code/preproc/tapas_physio_get_template_from_pulses.m index 86170dde..bc82742a 100644 --- a/PhysIO/code/preproc/tapas_physio_get_template_from_pulses.m +++ b/PhysIO/code/preproc/tapas_physio_get_template_from_pulses.m @@ -1,4 +1,4 @@ -function [pulseCleanedTemplate, pulseTemplate] = tapas_physio_get_template_from_pulses(... +function [pulseCleanedTemplate, pulseTemplate, verbose] = tapas_physio_get_template_from_pulses(... c, cpulse, halfTemplateWidthInSamples, verbose, varargin) % Computes mean pulse template given time stamp of detected pulses and raw % time series; removes outlier pulses with correlation to initial guess of @@ -94,19 +94,8 @@ % template as average of the found representative waves pulseTemplate = mean(template); +tTemplate = dt*(0:2*halfTemplateWidthInSamples); -if doDebug - fh = verbose.fig_handles(end); - figure(fh); - subplot(3,1,3); - tTemplate = dt*(0:2*halfTemplateWidthInSamples); - plot(tTemplate, template'); - hold all; - hp(1) = plot(tTemplate, pulseTemplate', '.--g', 'LineWidth', 3, 'Marker', ... - 'o'); - xlabel('t (seconds)'); - title('Templates of cycle time course and mean template'); -end %% Final template for peak search from best-matching templates % delete the peaks deviating from the mean too @@ -140,9 +129,15 @@ pulseCleanedTemplate = mean(template(indHighQualityTemplates, :)); if doDebug - stringTitle = 'Preproc: Iterative Template Creation Single Cycle'; - hp(2) = plot(tTemplate, pulseCleanedTemplate, '.-g', 'LineWidth', 4, ... - 'Marker', 'x'); - legend(hp, 'mean of templates', 'mean of most similar, chosen templates'); - tapas_physio_suptitle(stringTitle); + + [verbose] = tapas_physio_plot_templates_of_cycle_time(tTemplate, ... + template, pulseTemplate, pulseCleanedTemplate, verbose); + + %save relevant structures for review + + verbose.review.temp_cyc.tTemplate = tTemplate; + verbose.review.temp_cyc.template = template; + verbose.review.temp_cyc.pulseTemplate = pulseTemplate; + verbose.review.temp_cyc.pulseCleanedTemplate = pulseCleanedTemplate; + end \ No newline at end of file diff --git a/PhysIO/code/readin/tapas_physio_read_columnar_textfiles.m b/PhysIO/code/readin/tapas_physio_read_columnar_textfiles.m index c46bd02c..e3dfc9c8 100644 --- a/PhysIO/code/readin/tapas_physio_read_columnar_textfiles.m +++ b/PhysIO/code/readin/tapas_physio_read_columnar_textfiles.m @@ -2,11 +2,11 @@ % Reads _PULS, _RESP, _ECG, _Info-files from Siemens tics format with % multiple numbers of columns and different column headers % -% output = tapas_physio_read_columnar_textfiles(input) +% [C, columnNames] = tapas_physio_read_columnar_textfiles(fileName, fileType) % % IN % fileName *.log from Siemens VD/VE tics file format -% fileType 'ECG', 'PULS', 'RESP', 'Info', 'BIOPAC_TXT' +% fileType 'ECG', 'PULS', 'RESP', 'Info', 'BIOPAC_TXT', 'PHILIPS' % If not specified, this is read from the last part of the % filename after the last underscore, e.g. % Physio_*_ECG.log -> log @@ -56,6 +56,16 @@ switch upper(fileType) + case {'ADINSTRUMENTS_TXT', 'LABCHART_TXT'} + strColumnHeader = 'ChannelTitle='; + + %sub-02 from CUBRIC + parsePatternPerNColumns{4} = '%f %f %f %f'; + nEmptyLinesAfterHeader(4) = 1; + + %sub-01 from CUBRIC + parsePatternPerNColumns{7} = '%f %f %f %f %f %f %f'; + nEmptyLinesAfterHeader(7) = 4; case 'BIDS' strColumnHeader = ''; parsePatternPerNColumns{3} = '%f %f %f'; @@ -80,6 +90,14 @@ nEmptyLinesAfterHeader(3) = 0; nEmptyLinesAfterHeader(4) = 1; nEmptyLinesAfterHeader(5) = 1; + case 'PHILIPS' + strColumnHeader = 'v1raw'; % Philips header similar to: # v1raw v2raw v1 v2 ppu resp vsc gx gy gz mark mark2 + parsePatternPerNColumns{10} = '%d %d %d %d %d %d %d %d %d %s'; + parsePatternPerNColumns{11} = '%d %d %d %d %d %d %d %d %d %s %s'; % log file version 2 with two marker columns + parsePatternPerNColumns{12} = '%d %d %d %d %d %d %d %d %d %d %s %s'; % logfile version 3(?, from Nottingham) with 'vsc' column + nEmptyLinesAfterHeader(10) = 0; + nEmptyLinesAfterHeader(11) = 0; + nEmptyLinesAfterHeader(12) = 0; case {'PULS', 'RESP', 'ECG'} % have similar format % Cologne (RESP/PULS/ECG for 2nd column): % Time_tics RESP Signal @@ -104,20 +122,27 @@ while ~haveFoundColumnHeader nHeaderLines = nHeaderLines + 1; strLine = fgets(fid); - haveFoundColumnHeader = any(regexp(upper(strLine), strColumnHeader)); + haveFoundColumnHeader = any(regexpi(strLine, strColumnHeader)); end switch upper(fileType) + case {'ADINSTRUMENTS_TXT', 'LABCHART_TXT'} + columnNames = regexp(strLine, '([\t])', 'split'); + nColumns = numel(columnNames); case 'BIDS' % will be in separate json-file columnNames = {}; nColumns = 3; case 'BIOPAC_TXT' % bad column names with spaces...e.g. 'RESP - RSP100C' columnNames = regexp(strLine, '([\t])', 'split'); nColumns = numel(columnNames); + case 'PHILIPS' + columnNames = regexp(strLine, '([\w]*)', 'tokens'); + columnNames = [columnNames{:}]; % cell of cell into cell of strings + nColumns = numel(columnNames); otherwise columnNames = regexp(strLine, '([\w]*)', 'tokens'); columnNames = [columnNames{:}]; % cell of cell into cell of strings - nColumns = numel(regexp(strLine, ' *')) + 1; % columns are separated by arbitrary number of spaces + nColumns = numel(regexp(strLine, ' *')) + 1; % columns are separated by arbitrary number of spaces...TODO: Why + 1? end fclose(fid); @@ -126,8 +151,10 @@ % now read the rest of the file fid = fopen(fileName); switch upper(fileType) - case {'BIDS', 'BIOPAC_TXT', 'INFO', 'PULS', 'RESP'} + case {'ADINSTRUMENTS_TXT', 'LABCHART_TXT', 'BIDS', 'BIOPAC_TXT', 'INFO', 'PULS', 'RESP'} C = textscan(fid, parsePatternPerNColumns{nColumns}, 'HeaderLines', nHeaderLines); + case 'PHILIPS' % sometimes odd lines with single # occur within text file + C = textscan(fid, parsePatternPerNColumns{nColumns}, 'HeaderLines', nHeaderLines, 'CommentStyle', '#'); case {'ECG'} C = textscan(fid, parsePatternPerNColumns{nColumns}, 'HeaderLines', nHeaderLines); if nColumns == 4 % CMRR, different ECG channels possible! diff --git a/PhysIO/code/readin/tapas_physio_read_physlogfiles.m b/PhysIO/code/readin/tapas_physio_read_physlogfiles.m index 8b3519a6..840a5925 100755 --- a/PhysIO/code/readin/tapas_physio_read_physlogfiles.m +++ b/PhysIO/code/readin/tapas_physio_read_physlogfiles.m @@ -59,6 +59,9 @@ end switch lower(log_files.vendor) + case {'adinstruments_txt', 'labchart_txt'} + [c, r, t, cpulse, acq_codes] = ... + tapas_physio_read_physlogfiles_adinstruments_txt(log_files, cardiac_modality, verbose); case 'bids' [c, r, t, cpulse, acq_codes] = ... tapas_physio_read_physlogfiles_bids(log_files, cardiac_modality, verbose); diff --git a/PhysIO/code/readin/tapas_physio_read_physlogfiles_adinstruments_txt.m b/PhysIO/code/readin/tapas_physio_read_physlogfiles_adinstruments_txt.m new file mode 100644 index 00000000..56ff2031 --- /dev/null +++ b/PhysIO/code/readin/tapas_physio_read_physlogfiles_adinstruments_txt.m @@ -0,0 +1,250 @@ +function [c, r, t, cpulse, acq_codes, verbose, gsr] = tapas_physio_read_physlogfiles_adinstruments_txt(... + log_files, cardiac_modality, verbose, varargin) +% Reads in 4-column txt-export from ADInstruments/LabChart data +% (e.g., channels titled O2 CO2 Pulse Respiration Volume trigger +% trans_force) +% +% [c, r, t, cpulse, acq_codes, verbose] = tapas_physio_read_physlogfiles_adinstruments_txt(... +% log_files, cardiac_modality, verbose, varargin) +% +% IN log_files +% .log_cardiac *.txt file, contains header and several columns of the form +% ChannelTitle= Pulse Respiration Volume trigger +% Range= 20.000 mV 10.000 V 10.000 V +% 0 8.953 1.245 0.006 +% 0.001 9.609 1.255 0.006 +% 0.002 8.684 1.244 0.007 +% 0.003 9.508 1.250 0.006 +% 0.004 7.490 1.244 0.007 +% 0.005 7.830 1.250 0.006 +% 0.006 3.668 1.239 0.006 +% 0.007 3.283 1.245 0.007 +% 0.008 1.652 1.258 0.007 +% .log_respiration same as .log_cardiac +% .sampling_interval sampling interval (in seconds) +% default: 1 ms (1000 Hz) +% cardiac_modality 'ECG' or 'PULS'/'PPU'/'OXY' to determine +% which channel data to be returned +% UNUSED, is always pulse plethysmographic unit +% for BioPac +% verbose +% .level debugging plots are created if level >=3 +% .fig_handles appended by handle to output figure +% +% OUT +% cpulse time events of R-wave peak in cardiac time series (seconds) +% , since not written to logfile +% r respiratory time series +% t vector of time points (in seconds) +% c cardiac time series (PPU) +% acq_codes slice/volume start events marked by number <> 0 +% for time points in t +% 10/20 = scan start/end; +% 1 = ECG pulse; 2 = OXY max; 4 = Resp trigger; +% 8 = scan volume trigger (on) +% 16 = scan volume trigger (off) +% gsr galvanic skin response (not used) +% +% EXAMPLE +% tapas_physio_read_physlogfiles_biopac_txt +% +% See also tapas_physio_read_physlogfiles_siemens tapas_physio_plot_raw_physdata_siemens_hcp + +% Author: Lars Kasper +% Created: 2018-09-27 +% Copyright (C) 2018 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 . + +%% read out values +DEBUG = verbose.level >= 2; +doDebugTriggers = verbose.level >= 3; +doReplaceNans = true; + +% if true, assume that trigger switches between +5V and 0 for start of one +% volume, and back to 5V at start of next volume +hasAlternatingTriggerFlank = true; + +hasRespirationFile = ~isempty(log_files.respiration); +hasCardiacFile = ~isempty(log_files.cardiac); + +hasRespirationFile = ~isempty(log_files.respiration); +hasCardiacFile = ~isempty(log_files.cardiac); + +if hasCardiacFile + fileName = log_files.cardiac; +elseif hasRespirationFile + fileName = log_files.respiration; +end + + +[C, columnNames] = tapas_physio_read_columnar_textfiles(fileName, 'ADINSTRUMENTS_TXT'); + +% determine right column +switch lower(cardiac_modality) + case {'ecg_raw', 'ecg1_raw', 'v1raw'} + labelColCardiac = 'UNKNOWN'; %TODO: find out! + case {'oxy','oxyge', 'ppu'} + labelColCardiac = 'Pulse'; +end + +labelColResp = 'Respiration'; +labelColTrigger = 'Volume trigger'; + +% it just so works that the first column is just the time (in seconds) per +% row, but labeled "ChannelTitle=", so other column indices line up +% startsWith is needed to ignore trailing end-of-line for column names +if exist('startsWith', 'builtin') + idxColCardiac = find(startsWith(columnNames, labelColCardiac)); + idxColResp = find(startsWith(columnNames, labelColResp)); + idxColTrigger = find(startsWith(columnNames, labelColTrigger)); +else % older Matlab versions, checking chars only up to length of label + idxColCardiac = find(strncmpi(columnNames, labelColCardiac, numel(labelColCardiac))); + idxColResp = find(strncmpi(columnNames, labelColResp, numel(labelColResp))); + idxColTrigger = find(strncmpi(columnNames, labelColTrigger, numel(labelColTrigger))); +end + +thresholdTrigger = 4; % Volt, TTL trigger + +c = C{idxColCardiac}; +r = C{idxColResp}; +gsr = C{2}; % TODO: do correctly! +iAcqOn = (C{idxColTrigger}>thresholdTrigger); % trigger is 5V, but flips on/off between volumes + +% replace NaNs by max or min depending on nearest non-NaN-neighbour +% (whether it was closer to max or min) +if doReplaceNans + maxVal = max(c); + minVal = min(c); + idxNan = find(isnan(c)); % for loop + idxValid = find(~isnan(c)); + + % find nearest neighbors of valid indices, and replace with min/max, + % whatever value closest neighbor was closer to + if exist('knnsearch') + idxValidClosest = knnsearch(idxValid, idxNan); + validValClosest = c(idxValid(idxValidClosest)); + isValidValClosestCloserToMin = abs(validValClosest-maxVal) > abs(validValClosest - minVal); + c(idxNan(isValidValClosestCloserToMin)) = minVal; + c(idxNan(~isValidValClosestCloserToMin)) = maxVal; + else % slow...todo: optimize! + + nNans = numel(idxNan); + iNan = 1; + while iNan <= nNans + + if ~mod(iNan, 1000) + fprintf('%d/%d NaNs replaced\n', iNan, nNans); + end + + idx = idxNan(iNan); + + [~,iValidClosest] = min(abs(idxValid-idx)); + validValClosest = c(idxValid(iValidClosest)); + + % choose min or max valid value, whatever is closest + if abs(maxVal-validValClosest) > abs(validValClosest - minVal) + c(idx) = minVal; + else + c(idx) = maxVal; + end + + iNan = iNan + 1 + end + end +end + +%% Create timing vector from samples + +dt = log_files.sampling_interval; + +if isempty(dt) + dt = mean(diff(C{1})); % first column has timing vector for LabChart +end + +nSamples = max(numel(c), numel(r)); +t = -log_files.relative_start_acquisition + ((0:(nSamples-1))*dt)'; + +%% Recompute acq_codes as for Siemens (volume on/volume off) +acq_codes = []; + +if ~isempty(iAcqOn) % otherwise, nothing to read ... + % iAcqOn is a column of 1s and 0s, 1 whenever scan acquisition is on + % Determine 1st start and last stop directly via first/last 1 + % Determine everything else in between via difference (go 1->0 or 0->1) + iAcqStart = find(iAcqOn, 1, 'first'); + iAcqEnd = find(iAcqOn, 1, 'last'); + d_iAcqOn = diff(iAcqOn); + + % index shift + 1, since diff vector has index of differences i_(n+1) - i_n, + % and the latter of the two operands (i_(n+1)) has sought value +1 + iAcqStart = [iAcqStart; find(d_iAcqOn == 1) + 1]; + % no index shift, for the same reason + iAcqEnd = [find(d_iAcqOn == -1); iAcqEnd]; + + % remove duplicate entries + iAcqStart = unique(iAcqStart); + iAcqEnd = unique(iAcqEnd); + + if hasAlternatingTriggerFlank + % choose all rising and falling triggers as volume starts (instead + % of interpreting falling as an end of a trigger) + iAcqStartNew = sort([iAcqStart;iAcqEnd]); + iAcqEndNew = []; + else + iAcqStartNew = iAcqStart; + iAcqEndNew = iAcqEnd; + end + + + if doDebugTriggers + figure;plot(C{1},C{idxColTrigger}) + hold all; + stem(C{1}(iAcqStart), C{idxColTrigger}(iAcqStart)) + stem(C{1}(iAcqEnd), C{idxColTrigger}(iAcqEnd)) + stem(C{1}(iAcqStartNew), 0.8*C{idxColTrigger}(iAcqStartNew)) + + title('Trigger Debugging') + legend('Volume Trigger Signal', 'Trigger Rising Flank', ... + 'Trigger Falling Flank', 'Chosen Trigger Volume Start') + end + + iAcqStart = iAcqStartNew; + iAcqEnd = iAcqEndNew; + + acq_codes = zeros(nSamples,1); + acq_codes(iAcqStart) = 8; % to match Philips etc. format + acq_codes(iAcqEnd) = 16; % don't know... + + % report estimated onset gap between last slice of volume_n and 1st slice of + % volume_(n+1) + nAcqStarts = numel(iAcqStart); + nAcqEnds = numel(iAcqEnd); + nAcqs = min(nAcqStarts, nAcqEnds); + + if nAcqs >= 1 + % report time of acquisition, as defined in SPM + TA = mean(t(iAcqEnd(1:nAcqs)) - t(iAcqStart(1:nAcqs))); + verbose = tapas_physio_log(... + sprintf('TA = %.4f s (Estimated time of acquisition during one volume TR)', ... + TA), verbose, 0); + end +end + + +%% Plot, if wanted + +if DEBUG + stringTitle = 'Read-In: Raw ADInstruments/LabChart physlog data (TXT Export)'; + verbose.fig_handles(end+1) = ... + tapas_physio_plot_raw_physdata_siemens_hcp(t, c, r, acq_codes, ... + stringTitle); +end + +%% Undefined output parameters + +cpulse = []; diff --git a/PhysIO/code/readin/tapas_physio_read_physlogfiles_bids.m b/PhysIO/code/readin/tapas_physio_read_physlogfiles_bids.m index 0e8a58e6..386898a1 100644 --- a/PhysIO/code/readin/tapas_physio_read_physlogfiles_bids.m +++ b/PhysIO/code/readin/tapas_physio_read_physlogfiles_bids.m @@ -12,7 +12,7 @@ % physiological recordings can be found here: % % https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/ -% 06-physiological-and-other-continous-recordings.html +% 06-physiological-and-other-continuous-recordings.html % % [c, r, t, cpulse, acq_codes, verbose] = tapas_physio_read_physlogfiles_biopac_txt(... % log_files, cardiac_modality, verbose, varargin) diff --git a/PhysIO/code/readin/tapas_physio_read_physlogfiles_philips.m b/PhysIO/code/readin/tapas_physio_read_physlogfiles_philips.m index 861e1545..a5c66d89 100644 --- a/PhysIO/code/readin/tapas_physio_read_physlogfiles_philips.m +++ b/PhysIO/code/readin/tapas_physio_read_physlogfiles_philips.m @@ -11,7 +11,7 @@ % .cardiac % .sampling_interval % .relative_start_acquisition -% cardiac_modality +% cardiac_modality % 'ecg1_filtered' filtered 1st ECG channel signal % (Default) % 'ecg2_filtered' filteered 2nd ECG channel @@ -24,7 +24,7 @@ % Note: for ECG, append '_wifi' % for adjusted sampling rate in % wireless Ingenia Scanners -% +% % % OUT % c cardiac time series (ECG or pulse oximetry) @@ -33,8 +33,8 @@ % cpulse time events of R-wave peak in cardiac time series (seconds) % acq_codes slice/volume start events marked by number <> 0 % for time points in t -% 10/20 = scan start/end; -% 1 = ECG pulse; 2 = OXY max; 4 = Resp trigger; +% 10/20 = scan start/end; +% 1 = ECG pulse; 2 = OXY max; 4 = Resp trigger; % 8 = scan volume trigger % % EXAMPLE @@ -64,8 +64,31 @@ end if hasCardiac || hasResp - y = tapas_physio_read_physlogfiles_philips_matrix(logfile); - acq_codes = y(:,10); + % y = tapas_physio_read_physlogfiles_philips_matrix(logfile); + [C, columnNames] = tapas_physio_read_columnar_textfiles(logfile, 'PHILIPS'); + + idxColHex = find(contains(columnNames, 'mark')); % mark and mark2 are hexadecimal columns + idxColDec = setdiff(1:numel(columnNames), idxColHex); + C(idxColDec) = cellfun(@double, C(idxColDec), 'UniformOutput', false); + + + %% Convert hexadecimal acquisition codes + for iCol = idxColHex + C{iCol} = hex2dec(C{iCol}); + end + + + %% Account for incomplete rows + nMinRows = min(cellfun(@numel,C)); + C = cellfun(@(x) x(1:nMinRows), C, 'UniformOutput', false); + y = cell2mat(C); + + idxColAcqCodes = idxColHex(1); % first 'mark' column + idxColResp = find(contains(columnNames, 'resp')); % first 'mark' column + acq_codes = y(:,idxColAcqCodes); + + + else y = []; acq_codes = []; @@ -101,10 +124,10 @@ cpulse = find(acq_codes==1); if ~isempty(cpulse) cpulse = t(cpulse); -end; +end if hasResp - r = y(:,6); + r = y(:,idxColResp); else r = []; end @@ -129,18 +152,20 @@ iModality = iModality + 1; cardiac_modality = cardiacModalityArray{iModality}; switch lower(cardiac_modality) - case {'ecg_raw', 'ecg1_raw'} - c = y(:,1); - case {'ecg2_raw'} - c = y(:,2); - case {'ecg1', 'ecg_filtered', 'ecg1_filtered'} - c = y(:,3); - case { 'ecg2', 'ecg2_filtered'} - c = y(:,4); - case {'oxy','oxyge', 'ppu'} - c = y(:,5); + case {'ecg_raw', 'ecg1_raw', 'v1raw'} + labelColCardiac = 'v1raw'; + case {'ecg2_raw', 'v2raw'} + labelColCardiac = 'v2raw'; + case {'ecg1', 'ecg_filtered', 'ecg1_filtered', 'v1'} + labelColCardiac = 'v1'; + case { 'ecg2', 'ecg2_filtered', 'v2'} + labelColCardiac = 'v2'; + case {'oxy','oxyge', 'ppu'} + labelColCardiac = 'ppu'; end + idxColCardiac = find(strcmpi(columnNames, labelColCardiac)); + c = y(:, idxColCardiac); hasValidCardiacReadout = any(c); end diff --git a/PhysIO/code/readin/tapas_physio_read_physlogfiles_siemens.m b/PhysIO/code/readin/tapas_physio_read_physlogfiles_siemens.m index 0050768a..2e370de8 100755 --- a/PhysIO/code/readin/tapas_physio_read_physlogfiles_siemens.m +++ b/PhysIO/code/readin/tapas_physio_read_physlogfiles_siemens.m @@ -68,7 +68,7 @@ % used channel depends on cardiac modality switch cardiac_modality case 'ECG' - defaults.ecgChannel = 'mean'; %'mean'; 'v1'; 'v2' + defaults.ecgChannel = 'mean'; %'mean'; 'v1'; 'v2'; 'v3'; 'v4' otherwise defaults.ecgChannel = 'v1'; end @@ -139,31 +139,27 @@ [lineData, logFooter] = tapas_physio_read_physlogfiles_siemens_raw(... log_files.cardiac); - tLogTotal = logFooter.LogStopTimeSeconds - logFooter.LogStartTimeSeconds; + tLogTotal = logFooter.StopTimeSeconds - logFooter.StartTimeSeconds; if hasScanTimingDicomImage - tStartScan = tStartScanDicom; % this is just the start of the selected DICOM volume - tStopScan = tStopScanDicom; + tStartScan = tStartScanDicom; % this is the start of the DICOM volume selected for sync + tStopScan = tStopScanDicom; % this is the end time (start + TR) of the DICOM volume selected for sync else - % Just different time scale, gives bad scaling in plots, and not - % needed... - % tStartScan = logFooter.ScanStartTimeSeconds; - % tStopScan = logFooter.ScanStopTimeSeconds; - tStartScan = logFooter.LogStartTimeSeconds; - tStopScan = logFooter.LogStopTimeSeconds; + tStartScan = logFooter.StartTimeSeconds; + tStopScan = logFooter.StopTimeSeconds; end switch log_files.align_scan case 'first' relative_start_acquisition = tStartScan ... - - logFooter.LogStartTimeSeconds; + - logFooter.StartTimeSeconds; case 'last' % shift onset of first scan by knowledge of run duration and % onset of last scan in run relative_start_acquisition = ... (tStopScan - sqpar.Nscans*sqpar.TR) ... - - logFooter.LogStartTimeSeconds; + - logFooter.StartTimeSeconds; end @@ -214,30 +210,26 @@ if hasRespData [lineData, logFooter] = tapas_physio_read_physlogfiles_siemens_raw(... log_files.respiration); - tLogTotal = logFooter.LogStopTimeSeconds - logFooter.LogStartTimeSeconds; + tLogTotal = logFooter.StopTimeSeconds - logFooter.StartTimeSeconds; if hasScanTimingDicomImage - tStartScan = tStartScanDicom; - tStopScan = tStopScanDicom; % is incorrect, use tStartScan + TR! - else - % Just different time scale, gives bad scaling in plots, and not - % needed... - % tStartScan = logFooter.ScanStartTimeSeconds; - % tStopScan = logFooter.ScanStopTimeSeconds; - tStartScan = logFooter.LogStartTimeSeconds; - tStopScan = logFooter.LogStopTimeSeconds; + tStartScan = tStartScanDicom; % this is the start of the DICOM volume selected for sync + tStopScan = tStopScanDicom; % this is the end time (start + TR) of the DICOM volume selected for sync + else + tStartScan = logFooter.StartTimeSeconds; + tStopScan = logFooter.StopTimeSeconds; end switch log_files.align_scan case 'first' relative_start_acquisition = tStartScan - ... - logFooter.LogStartTimeSeconds; + logFooter.StartTimeSeconds; case 'last' % shift onset of first scan by knowledge of run duration and % onset of last scan in run relative_start_acquisition = ... - (tStartScan - (sqpar.Nscans-1)*sqpar.TR) ... - - logFooter.LogStartTimeSeconds; + (tStopScan - sqpar.Nscans*sqpar.TR) ... + - logFooter.StartTimeSeconds; end diff --git a/PhysIO/code/readin/tapas_physio_read_physlogfiles_siemens_raw.m b/PhysIO/code/readin/tapas_physio_read_physlogfiles_siemens_raw.m index 4e69458e..0f478cc3 100644 --- a/PhysIO/code/readin/tapas_physio_read_physlogfiles_siemens_raw.m +++ b/PhysIO/code/readin/tapas_physio_read_physlogfiles_siemens_raw.m @@ -47,25 +47,32 @@ linesFooter = C{1}(2:end); lineData = C{1}{1}; -%Get time stamps from footer: +% Get time stamps from footer: -logFooter.LogStartTimeSeconds = str2num(char(regexprep(linesFooter(~cellfun(@isempty,strfind(linesFooter,... +% MPCU = Computer who controls the physiological logging in real-time => physio logging happens here +% MDH = Computer who is the host (Measurement Data Header); console => DICOM time stamp here +logFooter.StartTimeSecondsScannerClock = str2num(char(regexprep(linesFooter(~cellfun(@isempty,strfind(linesFooter,... 'LogStartMDHTime'))),'\D',''))) / 1000; -logFooter.LogStopTimeSeconds = str2num(char(regexprep(linesFooter(~cellfun(@isempty,strfind(linesFooter,... +logFooter.StopTimeSecondsScannerClock = str2num(char(regexprep(linesFooter(~cellfun(@isempty,strfind(linesFooter,... 'LogStopMDHTime'))),'\D',''))) / 1000; -% MPCU = Computer who controls the scanner => physio logging happens here -% MDH = Compute who is the host; console => DICOM time stamp here! -% -% - according to Chris Rorden, PART -% (http://www.mccauslandcenter.sc.edu/crnl/tools/part) - MDH is the time we -% should use for phys logging synchronization, since DICOM conversion uses -% this clock +logFooter.StartTimeSecondsRecordingClock = str2num(char(regexprep(linesFooter(~cellfun(@isempty,strfind(linesFooter,... + 'LogStartMPCUTime'))),'\D',''))) / 1000; +logFooter.StopTimeSecondsRecordingClock = str2num(char(regexprep(linesFooter(~cellfun(@isempty,strfind(linesFooter,... + 'LogStopMPCUTime'))),'\D',''))) / 1000; -% This is just a different time-scale (of the phys log computer), it does -% definitely NOT match with the Acquisition time in the DICOM-headers -logFooter.ScanStartTimeSeconds = str2num(char(regexprep(linesFooter(~cellfun(@isempty,strfind(linesFooter,... - 'LogStartMPCUTime'))),'\D',''))); -logFooter.ScanStopTimeSeconds = str2num(char(regexprep(linesFooter(~cellfun(@isempty,strfind(linesFooter,... - 'LogStopMPCUTime'))),'\D',''))); +% +% We use the time stamp of the clock of the Measurement Data Header (MDH), +% i.e., computer that controls the scanner, to synchronize with the DICOMs, +% because this computer also controls the creation of the scan data, i.e., +% reconstructed DICOM images. This is in accordance to other packages +% reading Siemens physiological logfile data, e.g., Chris Rorden's PART +% (https://github.com/neurolabusc/Part#usage), +% with a detailed explanation on the DICOM timestamp in AcquisitionTime +% found here (https://github.com/nipy/heudiconv/issues/450#issuecomment-645003447) +% +% MPCU is the clock of the computer that controls the physiological +% recording (same as MARS?), but does not know about the scan volume and DICOM timing +logFooter.StartTimeSeconds = logFooter.StartTimeSecondsScannerClock; +logFooter.StopTimeSeconds = logFooter.StopTimeSecondsScannerClock; \ No newline at end of file diff --git a/PhysIO/code/readin/tapas_physio_siemens_line2table.m b/PhysIO/code/readin/tapas_physio_siemens_line2table.m index 103702d6..a65349e4 100755 --- a/PhysIO/code/readin/tapas_physio_siemens_line2table.m +++ b/PhysIO/code/readin/tapas_physio_siemens_line2table.m @@ -1,6 +1,6 @@ -function data_table = tapas_physio_siemens_line2table(lineData, cardiacModality) -% transforms data line of Siemens log file into table (sorting amplitude -% and trigger signals) +function [data_table, log_parts] = tapas_physio_siemens_line2table(lineData, cardiacModality) +% transforms data line of Siemens log file (before linebreak after 5003 +% signal for recording end) into table (sorting amplitude and trigger signals) % % data_table = tapas_physio_siemens_line2table(input) % @@ -9,12 +9,13 @@ % of file. See also tapas_physio_read_physlogfiles_siemens_raw % % OUT -% data_table [nSamples,3] table of channel_1, channels_AVF and trigger +% data_table [nSamples,nChannels+1] table of recording channel_1, ..., channel_N and trigger % signal with trigger codes: % 5000 = cardiac pulse on % 6000 = cardiac pulse off % 6002 = phys recording on % 6003 = phys recording off +% log_parts part of logfile according to markers by Siemens % % EXAMPLE % tapas_physio_siemens_line2table @@ -32,120 +33,160 @@ % COPYING or . -% signals start of data logging -iTrigger = regexpi(lineData, ' 6002 '); +% NOTE: The following extraction of data from the logfiles might seem +% parsimonious, but it's written in a way to support all different physio +% trace modalities and logfile versions (more detail below) + +% The typical logfile structure is as follows (all data in first line of +% logfile, the footer is in the next line (after 5003), not used in this +% file, but tapas_physio_read_physlogfiles_siemens_raw) +% +%
5002 6002 +% <[optional] training trace data> 5002 uiHwRevisionPeru ... [optional] 6002 +% 5002 6002 5002 6002 ... 5002 6002 +% ... +% 5002 6002 +% ... +% 5002 6002 ... +% ... 5003 +% +% Since the number and content of actual info regions and number of data +% channels differs by trace modality and logfile version, we have to cut +% data out carefully: +% +% 1. Cut away header +% 2. Cut away logfile version +% 3. Cut away optional training trace +% 4. Remove all infoRegions can be interleaved with trace data, e.g., +% cushionGain for RESP trace) +% 5. Remove all trigger markers (5000, 6000), but remmember position +% 6. Sort remaining trace into coresponding channels (number of channels is +% logfile-version dependent) +% 7. Re-insert trigger markers as extra column -if ~isempty(iTrigger) - % crop string after trigger - lineData = lineData((iTrigger(end)+6):end); - doCropLater = false; -else - % crop first 4 values as in UPenn manual after conversion - doCropLater = true; -end -data = textscan(lineData, '%d', 'Delimiter', ' ', 'MultipleDelimsAsOne',1); -if doCropLater - % crop first 4 values; - data{1} = data{1}(5:end); +%% 1. Header goes from start of line to first occurence of 5002 +[iStartHeader, iEndHeader, logHeader] = regexp(lineData, '^(.*?) (?=\<5002\>)', 'start', 'end', 'tokens' ); +logHeader = logHeader{1}{1}; +lineData(iStartHeader:iEndHeader) = []; % remove header, but keep 5002 for next step + + +%% 2. Logfile version (which alters no of channels etc.) +% stored in first 5002/6002 info region +% 5002 LOGVERSION 1 6002% +% 5002 LOGVERSION_RESP 3 6002% +% 5002 LOGVERSION_PULS 3 6002% +[iStartVersionInfo, iEndVersionInfo, logVersion] = regexp(lineData, '^\<5002\> LOGVERSION[^0-9]*(\d+)\s\<6002\>', 'start', 'end', 'tokens' ); +logVersion = str2double(logVersion{1}{1}); +lineData(iStartVersionInfo:iEndVersionInfo) = []; % remove version info region incl. 5002/6002 markers + + +%% 3. Optional training data (for Siemens own peak detection) is after +% "6002" of Logversion info and "5002 uiHwRevisionPeru" +[iStartTraining, iEndTraining, dataTraining] = regexp(lineData, '^\s*(.*?) (?=\<5002\> uiHwRevisionPeru)', 'start', 'end', 'tokens'); +if ~isempty(iStartTraining) % training trace does not always exist + dataTraining = dataTraining{1}{1}; + lineData(iStartTraining:iEndTraining) = []; % remove training trace, but keep following 5002 for next step end + +%% 4. Identify and remove info regions between 5002 and 6002 (may be +% interleaved with trace data (e.g., messages or cushion Gain for RESP) +[iStartInfoRegion, iEndInfoRegion, logInfoRegion] = regexp(lineData, '\<5002\>(.*?)\<6002\>', 'start', 'end', 'tokens' ); +logInfoRegion = cellfun(@(x) x{1,1}, logInfoRegion, 'UniformOutput', false)'; +traceData = regexprep(lineData, '\<5002\>(.*?)\<6002\>\s', ''); +traceData = regexprep(traceData, '\<5003\>$', ''); % remove 5003 mark of trace end + +log_parts.logHeader = logHeader; +log_parts.logInfoRegion = logInfoRegion; +log_parts.logVersion = logVersion; + +% convert remaining data (all numbers string) to number (int32) +data = textscan(traceData, '%d', 'Delimiter', ' ', 'MultipleDelimsAsOne', true); + + +%% 5. Remove all trigger markers (5000, 6000), but remmember position % Remove the systems own evaluation of triggers. cpulse = find(data{1} == 5000); % System uses identifier 5000 as trigger ON cpulse_off = find(data{1} == 6000); % System uses identifier 5000 as trigger OFF -recording_on = find(data{1} == 6002);% Scanner trigger to Stim PC? -recording_off = find(data{1} == 5003); - - % Filter the trigger markers from the ECG data % Note: depending on when the scan ends, the last size(t_off)~=size(t_on). -iNonEcgSignals = [cpulse; cpulse_off; recording_on; recording_off]; -codeNonEcgSignals = [5000*ones(size(cpulse)); ... - 6000*ones(size(cpulse_off)); ... - 6002*ones(size(recording_on)) - 5003*ones(size(recording_off))]; - -% data_stream contains only the 2 ECG-channel time courses (with -% interleaved samples +iTriggerMarker = [cpulse; cpulse_off]; +codeTriggerMarker = [5000*ones(size(cpulse)); ... + 6000*ones(size(cpulse_off))]; + +% data_stream contains only the time courses (with +% interleaved samples for each channel) data_stream = data{1}; -data_stream(iNonEcgSignals) = []; +data_stream(iTriggerMarker) = []; -%iDataStream contains the indices of all true ECG signals in the full -%data{1}-Array that contains also non-ECG-signals +%iDataStream contains the indices of all true trace signals in the full +%data{1}-Array that contains also the trigger markers iDataStream = 1:numel(data{1}); -iDataStream(iNonEcgSignals) = []; +iDataStream(iTriggerMarker) = []; -nSamples = numel(data_stream); +%% 6. Sort remaining trace into coresponding channels (number of channels is +% logfile-version dependent) +nSamples = numel(data_stream); switch upper(cardiacModality) % ecg has two channels, resp and puls only one case 'ECG' - nRows = ceil(nSamples/2); - - % create a table with channel_1, channels_AVF and trigger signal in - % different Columns - % - iData_table is a helper table that maps the original indices of the - % ECG signals in data{1} onto their new positions - data_table = zeros(nRows,3); - iData_table = zeros(nRows,3); - - data_table(1:nRows,1) = data_stream(1:2:end); - iData_table(1:nRows,1) = iDataStream(1:2:end); - - if mod(nSamples,2) == 1 - data_table(1:nRows-1,2) = data_stream(2:2:end); - iData_table(1:nRows-1,2) = iDataStream(2:2:end); - else - data_table(1:nRows,2) = data_stream(2:2:end); - iData_table(1:nRows,2) = iDataStream(2:2:end); + switch logVersion + case 1 + nChannels = 2; + case 3 + nChannels = 4; end - - % now fill up 3rd column with trigger data - % - for each trigger index in data{1}, check where ECG data with closest - % smaller index ended up in the data_table ... and put trigger code in - % same row of that table - nTriggers = numel(iNonEcgSignals); - - for iTrigger = 1:nTriggers - % find index before trigger event in data stream and - % detect it in table - iRow = find(iData_table(:,2) == iNonEcgSignals(iTrigger)-1); - - % look in 1st column as well whether maybe signal detected there - if isempty(iRow) - iRow = find(iData_table(:,1) == iNonEcgSignals(iTrigger)-1); - end - - data_table(iRow,3) = codeNonEcgSignals(iTrigger); + case 'PPU' + nChannels = 1; + case 'RESP' + switch logVersion + case 1 + nChannels = 1; + case 3 + nChannels = 5; % breathing belt plus 4 channel biomatrix end - case {'RESP', 'PPU'} % only one channel available, fill second row with zeros - nRows = nSamples; - - % create a table with channel_1 and trigger signal in - % different Columns - % - iData_table is a helper table that maps the original indices of the - % ECG signals in data{1} onto their new positions - data_table = zeros(nRows,3); - iData_table = zeros(nRows,3); - - data_table(1:nRows,1) = data_stream; - iData_table(1:nRows,1) = iDataStream; - - % now fill up 3rd column with trigger data - % - for each trigger index in data{1}, check where ECG data with closest - % smaller index ended up in the data_table ... and put trigger code in - % same row of that table - nTriggers = numel(iNonEcgSignals); - - for iTrigger = 1:nTriggers - % find index before trigger event in data stream and - % detect it in table - iRow = find(iData_table(:,1) == iNonEcgSignals(iTrigger)-1); - if ~isempty(iRow) - data_table(iRow,3) = codeNonEcgSignals(iTrigger); - end - end otherwise error('unknown cardiac/respiratory logging modality: %s', cardiacModality); +end + +nRows = ceil(nSamples/nChannels); + +% create a table with channel_1, channels_AVF and trigger signal in +% different Columns +% - iData_table is a helper table that maps the original indices of the +% ECG signals in data{1} onto their new positions +data_table = zeros(nRows,nChannels+1); +iData_table = zeros(nRows,nChannels+1); + +for iChannel = 1:nChannels + data_table(1:nRows,iChannel) = data_stream(iChannel:nChannels:end); + iData_table(1:nRows,iChannel) = iDataStream(iChannel:nChannels:end); +end + +% TODO: deal with mod(nSamples, nChannels) > 0 (incomplete data?) + + +%% 7. Re-insert trigger markers as extra column +% now fill up nChannel+1. column with trigger data +% - for each trigger index in data{1}, check where ECG data with closest +% smaller index ended up in the data_table ... and put trigger code in +% same row of that table +nTriggers = numel(iTriggerMarker); + +for iTrigger = 1:nTriggers + % find index before trigger event in data stream and + % detect it in table, look in last columns first, then go + % backwards + iRow = []; + iChannel = nChannels; + while isempty(iRow) + + iRow = find(iData_table(:,iChannel) == iTriggerMarker(iTrigger)-1); + iChannel = iChannel - 1; + end + + data_table(iRow,nChannels+1) = codeTriggerMarker(iTrigger); end \ No newline at end of file diff --git a/PhysIO/code/readin/tapas_physio_siemens_table2cardiac.m b/PhysIO/code/readin/tapas_physio_siemens_table2cardiac.m index 7f07762a..9ae6cfea 100644 --- a/PhysIO/code/readin/tapas_physio_siemens_table2cardiac.m +++ b/PhysIO/code/readin/tapas_physio_siemens_table2cardiac.m @@ -18,19 +18,24 @@ % onset of first scan (t=0) % % OUT +% % dataCardiac = struct(... % 'cpulse_on', cpulse_on, ... % 'cpulse_off', cpulse_off, ... % 'recording_on', recording_on, ... % 'recording_off', recording_off, ... -% 'channel_1', channel_1, ... -% 'channel_AVF', channel_AVF, ... +% 'recordingChannels', recordingChannels, ... % 'meanChannel', meanChannel, ... % 'c', c, ... % 't', t, ... % 'ampl', ampl, ... % 'stopSample', stopSample ... % ); +% +% AVF means "Augmented Vector Foot" and is a label for a unipolar lead +% of the ECG electrode setup. For .ecg- LOGVERSIOM one, this was the 2nd +% column ("channel") in the data (based on Siemens' Physioload function), +% but for later versions, it's not clear % % EXAMPLE % tapas_physio_siemens_table2cardiac @@ -49,27 +54,32 @@ % COPYING or . -% set new indices to actual -cpulse_on = find(data_table(:,3) == 5000); -cpulse_off = find(data_table(:,3) == 6000); -recording_on = find(data_table(:,3) == 6002); -recording_off = find(data_table(:,3) == 5003); +% set new indices to actual, last column contains +cpulse_on = find(data_table(:,end) == 5000); +cpulse_off = find(data_table(:,end) == 6000); +recording_on = find(data_table(:,end) == 6002); +recording_off = find(data_table(:,end) == 5003); % Split a single stream of ECG data into channel 1 and channel 2. -channel_1 = data_table(:,1); -channel_AVF = data_table(:,2); -meanChannel = mean([channel_1(:) channel_AVF(:)],2); +nChannels = size(data_table,2) - 1; + +for iChannel = 1:nChannels + recordingChannels(:,iChannel) = data_table(:,iChannel); +end + +meanChannel = mean(recordingChannels, 2); % Make them the same length and get time estimates. -switch ecgChannel +switch lower(ecgChannel) case 'mean' c = meanChannel - mean(meanChannel); - case 'v1' - c = channel_1 - mean(channel_1); - - case 'v2' - c = channel_AVF - mean(channel_AVF); + case {'v1', 'v2', 'v3', 'v4'} + iChannel = str2double(ecgChannel(2)); + c = recordingChannels(:,iChannel) - mean(recordingChannels(:,iChannel)); + case {'avf'} + iChannel =2; + c = recordingChannels(:,iChannel) - mean(recordingChannels(:,iChannel)); end % compute timing vector and time of detected trigger/cpulse events @@ -94,8 +104,7 @@ 'cpulse_off', cpulse_off, ... 'recording_on', recording_on, ... 'recording_off', recording_off, ... - 'channel_1', channel_1, ... - 'channel_AVF', channel_AVF, ... + 'recordingChannels', recordingChannels, ... 'meanChannel', meanChannel, ... 'c', c, ... 't', t, ... diff --git a/PhysIO/code/sync/tapas_physio_create_scan_timing.m b/PhysIO/code/sync/tapas_physio_create_scan_timing.m index a374b4ac..078a51ca 100644 --- a/PhysIO/code/sync/tapas_physio_create_scan_timing.m +++ b/PhysIO/code/sync/tapas_physio_create_scan_timing.m @@ -101,7 +101,7 @@ [VOLLOCS, LOCS, verbose] = ... tapas_physio_create_scan_timing_from_tics_siemens( ... ons_secs.t, ons_secs.t_start, log_files, verbose); - case {'biopac_mat', 'biopac_txt', 'bids'} + case {'adinstruments_txt', 'labchart_txt', 'biopac_mat', 'biopac_txt', 'bids'} [VOLLOCS, LOCS, verbose] = ... tapas_physio_create_scan_timing_from_acq_codes( ... ons_secs.t + ons_secs.t_start, ons_secs.acq_codes, ... diff --git a/PhysIO/code/sync/tapas_physio_create_scan_timing_from_acq_codes.m b/PhysIO/code/sync/tapas_physio_create_scan_timing_from_acq_codes.m index e1c57d1b..b2ca03b3 100644 --- a/PhysIO/code/sync/tapas_physio_create_scan_timing_from_acq_codes.m +++ b/PhysIO/code/sync/tapas_physio_create_scan_timing_from_acq_codes.m @@ -56,6 +56,12 @@ VOLLOCS = find(acq_codes == 8); end +% bugfix if slice-wise triggers were mistaken for volume triggers +if isempty(LOCS) && numel(VOLLOCS) >= nTotalSlices + LOCS = VOLLOCS; + VOLLOCS = []; +end + isValidVOLLOCS = numel(VOLLOCS) >= nTotalVolumes; isValidLOCS = numel(LOCS) >= nTotalSlices; diff --git a/PhysIO/code/sync/tapas_physio_create_scan_timing_from_end_marker_philips.m b/PhysIO/code/sync/tapas_physio_create_scan_timing_from_end_marker_philips.m index 17680a45..8924984b 100644 --- a/PhysIO/code/sync/tapas_physio_create_scan_timing_from_end_marker_philips.m +++ b/PhysIO/code/sync/tapas_physio_create_scan_timing_from_end_marker_philips.m @@ -144,22 +144,8 @@ if verbose.level>=1 - verbose.fig_handles(end+1) = tapas_physio_get_default_fig_params(); - set(gcf,'Name', 'Sync: 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' ) + [verbose] = tapas_physio_plot_create_scan_timing_philips(t, y, VOLLOCS, Ndummies, LOCS, ... + LOC_END_MARKER) end % VOLLOCS = find(abs(diff(z2))>thresh.vol); diff --git a/PhysIO/code/sync/tapas_physio_get_onsets_from_locs.m b/PhysIO/code/sync/tapas_physio_get_onsets_from_locs.m index 629d07f9..49bda6b3 100644 --- a/PhysIO/code/sync/tapas_physio_get_onsets_from_locs.m +++ b/PhysIO/code/sync/tapas_physio_get_onsets_from_locs.m @@ -88,12 +88,13 @@ end if verbose.level >= 3 - stringTitle = 'Sync: Slice bundles belonging to 1 volume'; - verbose.fig_handles(end+1) = tapas_physio_get_default_fig_params(); - set(gcf, 'Name', stringTitle); - for v=1:Nallvols-1, stem(t(SLICELOCS{v}),ones(size(SLICELOCS{v})));hold all;end - title(stringTitle); - xlabel('t (seconds since SCANPHYSLOG-start)'); + verbose = tapas_physio_plot_sync_bundles(Nallvols, t,SLICELOCS, verbose) + + % save relevant parameters + verbose.review.sync_bundles.Nallvols = Nallvols; + verbose.review.sync_bundles.t = t; + verbose.review.sync_bundles.SLICELOCS =SLICELOCS; + end try diff --git a/PhysIO/code/tapas_physio_cfg_matlabbatch.m b/PhysIO/code/tapas_physio_cfg_matlabbatch.m index dca37a3f..9d7d732a 100644 --- a/PhysIO/code/tapas_physio_cfg_matlabbatch.m +++ b/PhysIO/code/tapas_physio_cfg_matlabbatch.m @@ -39,6 +39,7 @@ vendor.tag = 'vendor'; vendor.name = 'vendor'; vendor.help = {' Vendor Name depending on your MR Scanner/Physiological recording system' + ' ''ADInstruments/LabChart_Txt'' - Text file (.txt) export of .adinst files, for recordings from AD Instruments devices, typically viewed and processed in LabChart' ' ''BIDS'' - Brain Imaging Data Structure (https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/06-physiological-and-other-continous-recordings.html)' ' ''Biopac_Txt'' - exported txt files from Biopac system (4 columns, [Resp PPU GSR Trigger]' ' ''Biopac_Mat'' - exported mat files from Biopac system' @@ -70,14 +71,16 @@ ' HCP-downloaded files of name format *_Physio_log.txt ' ' are already preprocessed into this simple 3-column text format' }; -vendor.labels = {'BIDS (Brain Imaging Data Structure)', 'Biopac_Txt', 'Biopac_Mat', ... +vendor.labels = { 'ADInstruments/LabChart_Txt', ... + 'BIDS (Brain Imaging Data Structure)', 'Biopac_Txt', 'Biopac_Mat', ... 'BrainProducts', 'Custom', ... 'GE', 'Philips', ... 'Siemens (VB, *.puls/*.ecg/*.resp)', ... 'Siemens_Tics (VD: *_PULS.log/*_ECG1.log/*_RESP.log/*_AcquisitionInfo*.log)', ... 'Siemens_HCP (Human Connectome Project, *Physio_log.txt, 3 column format', ... }; -vendor.values = {'BIDS', 'Biopac_Txt','Biopac_Mat', 'BrainProducts', 'Custom', ... +vendor.values = {'ADInstruments_Txt', 'BIDS', 'Biopac_Txt','Biopac_Mat', ... + 'BrainProducts', 'Custom', ... 'GE', 'Philips', 'Siemens', 'Siemens_Tics', 'Siemens_HCP', ... }; vendor.val = {'Philips'}; diff --git a/PhysIO/code/tapas_physio_init.m b/PhysIO/code/tapas_physio_init.m index fd8e5a41..9a34919e 100644 --- a/PhysIO/code/tapas_physio_init.m +++ b/PhysIO/code/tapas_physio_init.m @@ -46,6 +46,12 @@ fprintf('OK.\n'); end +% Adding test paths as well to run tapas_physio_test via tapas_test +% TODO: does not work, if code folder in SPM/toolbox +if ~exist('tapas_physio_test') + pathPhysIOTest = fullfile(fileparts(pathPhysIO), 'tests'); + addpath(genpath(pathPhysIOTest)); +end %% Check and add SPM path fprintf('Checking Matlab SPM path now...'); diff --git a/PhysIO/code/tapas_physio_new.m b/PhysIO/code/tapas_physio_new.m index d0d65627..1704fac6 100644 --- a/PhysIO/code/tapas_physio_new.m +++ b/PhysIO/code/tapas_physio_new.m @@ -101,7 +101,12 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% log_files = []; - % vendor Name depending on your MR Scanner system + % vendor Name depending on your MR Scanner / + % physiological recording system + % 'ADInstruments/LabChart_Txt'' - Text file (.txt) + % export of .adinst files, for recordings + % from AD Instruments devices, typically + % viewed and processed in LabChart % 'BIDS' - Brain Imaging Data Structure (http://bids.neuroimaging.io/bids_spec.pdf, section 8.6)' % 'Biopac_Txt' - exported txt files from Biopac system (4 columns, [Resp PPU GSR Trigger]' % 'Biopac_Mat' - exported mat files from Biopac system' diff --git a/PhysIO/code/tapas_physio_version.m b/PhysIO/code/tapas_physio_version.m index a6db2ee0..512e1d63 100644 --- a/PhysIO/code/tapas_physio_version.m +++ b/PhysIO/code/tapas_physio_version.m @@ -23,4 +23,4 @@ % version 3 or, at your option, any later version). For further details, % see the file COPYING or . % -versionPhysio = 'R2022a-v8.1.0'; +versionPhysio = 'R2022a-v8.2.0-beta'; diff --git a/PhysIO/code/utils/tapas_physio_check_spm_batch_editor_integration.m b/PhysIO/code/utils/tapas_physio_check_spm_batch_editor_integration.m index 642f4ba3..1425914e 100644 --- a/PhysIO/code/utils/tapas_physio_check_spm_batch_editor_integration.m +++ b/PhysIO/code/utils/tapas_physio_check_spm_batch_editor_integration.m @@ -32,14 +32,22 @@ [isSpmOnPath, pathSpm] = tapas_physio_check_spm(); -isPhysioVisibleForSpmBatchEditor = isSpmOnPath; % minimum requirement for integration: SPM works! +isPhysioVisibleForSpmBatchEditor = false; % checked below, need SPM visible for that % check for config matlabbatch file if isSpmOnPath - filePhysioCfgMatlabbatch = ... - dir(fullfile(pathSpm, 'toolbox', '**/tapas_physio_cfg_matlabbatch.m')); - isPhysioVisibleForSpmBatchEditor = ~isempty(filePhysioCfgMatlabbatch); + % check all possible SPM toolbox locations, as listed in its defaults: + tbx = spm_get_defaults('tbx'); % SPM toolbox parameter struct + + iDir = 1; + while (iDir <= numel(tbx.dir)) && ~isPhysioVisibleForSpmBatchEditor + filePhysioCfgMatlabbatch = ... + dir(fullfile(tbx.dir{iDir}, '**/tapas_physio_cfg_matlabbatch.m')); + + isPhysioVisibleForSpmBatchEditor = ~isempty(filePhysioCfgMatlabbatch); + iDir = iDir + 1; + end % also important to set default modality of spm to fMRI and % initialize batch editor, if not done before diff --git a/PhysIO/code/utils/tapas_physio_fill_empty_parameters.m b/PhysIO/code/utils/tapas_physio_fill_empty_parameters.m index 3ad5be5d..a8f159e6 100644 --- a/PhysIO/code/utils/tapas_physio_fill_empty_parameters.m +++ b/PhysIO/code/utils/tapas_physio_fill_empty_parameters.m @@ -52,7 +52,7 @@ switch lower(physio.log_files.vendor) case {'bids', 'biopac_mat', 'brainproducts', 'siemens'} % will be read from file later physio.log_files.sampling_interval = []; - case 'biopac_txt' + case {'adinstruments_txt', 'labchart_txt', 'biopac_txt'} physio.log_files.sampling_interval = 1/1000; case 'ge' physio.log_files.sampling_interval = 25e-3; diff --git a/PhysIO/code/utils/tapas_physio_pca.m b/PhysIO/code/utils/tapas_physio_pca.m index 34d996a9..a1a313e1 100644 --- a/PhysIO/code/utils/tapas_physio_pca.m +++ b/PhysIO/code/utils/tapas_physio_pca.m @@ -29,14 +29,17 @@ [nVoxels,nVolumes] = size(timeseries); -if nVoxels <= nVolumes - error([mfilename ':NotEnoughVoxels'], 'nVoxels <= nVolumes') -end - if nargin < 2 method = 'svd'; end +if nVoxels <= nVolumes + % TODO: reimplement using PCA of COV((X-mu)'*(X-mu)) + method = 'stats-pca'; + % error([mfilename ':NotEnoughVoxels'], 'nVoxels <= nVolumes') +end + + switch lower(method) case 'svd' @@ -56,8 +59,8 @@ LATENT = eigen_values; % [nPCs, 1] % Eigen_values -> Variance explained - vairance_explained = 100*eigen_values/sum(eigen_values); % in percent (%) - EXPLAINED = vairance_explained; % [nPCs, 1] + variance_explained = 100*eigen_values/sum(eigen_values); % in percent (%) + EXPLAINED = variance_explained; % [nPCs, 1] % Sign convention : the max(abs(PCs)) is positive [~,maxabs_idx] = max(abs(v)); diff --git a/PhysIO/docs/documentation.html b/PhysIO/docs/documentation.html index 8e66b2aa..da41aad7 100644 --- a/PhysIO/docs/documentation.html +++ b/PhysIO/docs/documentation.html @@ -825,9 +825,9 @@

Troubleshoot

Readme

TAPAS PhysIO Toolbox

-

Current version: Release 2022a, v8.1.0

+

Current version: Release 2022b, v8.2.0

-

Copyright (C) 2012-2022
Lars Kasper
kasper@biomed.ee.ethz.ch

+

Copyright (C) 2012-2022
Lars Kasper
kasper@biomed.ee.ethz.ch

Translational Neuromodeling Unit (TNU)
Institute for Biomedical Engineering
University of Zurich and ETH Zurich

Download

@@ -1523,29 +1523,21 @@

Example Datasets for PhysIO

The following datasets are available to explore the read-in and modeling capabilities of PhysIO. They can be downloaded by running the function tapas_download_example_data() in Matlab, which is located in the misc subfolder of the TAPAS software release you downloaded (probably here).

Afterwards, the examples can be found in tapas/examples/<tapasVersion>/PhysIO as different subfolders (vendor/device) and shall be run directly from within these individual folders.

-

Besides the raw physiological logfiles, each example contains example scripts to run PhysIO as

+

Besides the raw physiological logfiles, each example contains example scripts to run PhysIO as

    -
  • SPM job (*spm_job.mat)
  • -
  • editable SPM job (*spm_job.m)
  • -
  • plain matlab script (*matlab_script.m)
  • +
  • SPM job (\\\*spm_job.mat)
  • +
  • editable SPM job (\\\*spm_job.m)
  • +
  • plain matlab script (\\\*matlab_script.m)

Brain Imaging Data Structure (BIDS)

CPULSE 3T

-

Courtesy of Hrvoje Stojic, Max Planck UCL Centre for Computational Psychiatry -and Ageing Research, University College London

-

Vendor-computed (software: Spike2) cardiac pulse events from PPU (finger -plethysmograph) data, Siemens 3T scanner, Multiband CMRR sequence

-

Description: This datasets contains the (compressed) tab-separated value -(.tsv.gz) files as well as the meta-file (.json) holding sampling rate of -the physiological recording, and its relative onset to scanning, in adherence -with the BIDS standard for peripheral recordings -files.

+

Courtesy of Hrvoje Stojic, Max Planck UCL Centre for Computational Psychiatry and Ageing Research, University College London

+

Vendor-computed (software: Spike2) cardiac pulse events from PPU (finger plethysmograph) data, Siemens 3T scanner, Multiband CMRR sequence

+

Description: This datasets contains the (compressed) tab-separated value (.tsv.gz) files as well as the meta-file (.json) holding sampling rate of the physiological recording, and its relative onset to scanning, in adherence with the BIDS standard for peripheral recordings files.

PPU 3T

-

Courtesy of Hrvoje Stojic, Max Planck UCL Centre for Computational Psychiatry -and Ageing Research, University College London

+

Courtesy of Hrvoje Stojic, Max Planck UCL Centre for Computational Psychiatry and Ageing Research, University College London

PPU (finger plethysmograph) and breathing belt, Siemens 3T scanner, Multiband CMRR sequence

-

Description: Similar to CPULSE3T (same acquisition system), but now with analog -data instead of vendor-detected pulses, data from different subject

+

Description: Similar to CPULSE3T (same acquisition system), but now with analog data instead of vendor-detected pulses, data from different subject

PPU 3T Separate Files

Courtesy of Alexandre Sayal CIBIT, University of Coimbra

PPU (finger plethysmograph) and breathing belt, Siemens 3T scanner, Multiband CMRR sequence

@@ -1553,70 +1545,67 @@

PPU 3T Separate Files

General Electric

PPU 3T

Courtesy of Steffen Bollmann, Kinderspital Zurich and ETH Zurich

-

PPU (finger plethysmograph) and breathing belt, General Electric 3T -scanner

-

Description: Similar to PPU, but acquired on a GE system with two -separate output logfiles for pulse oximetry and breathing amplitude, -sampled with 40 Hz. The quality of the signal is particularly -challenging, stemming from a patient population.

+

PPU (finger plethysmograph) and breathing belt, General Electric 3T scanner

+

Description: Similar to PPU, but acquired on a GE system with two separate output logfiles for pulse oximetry and breathing amplitude, sampled with 40 Hz. The quality of the signal is particularly challenging, stemming from a patient population.

Philips

ECG 3T

-

Courtesy of Sandra Iglesias, Translational Neuromodeling Unit, ETH & -University of Zurich

+

Courtesy of Sandra Iglesias, Translational Neuromodeling Unit, ETH & University of Zurich

4-electrode ECG and breathing belt, Philips 3T Achieva scanner

-

Description: Standard example; shows how to use volume counting either -from beginning or end of run to synchronize physiological logfile with -acquisition onsets of fMRI scans.

+

Description: Standard example; shows how to use volume counting either from beginning or end of run to synchronize physiological logfile with acquisition onsets of fMRI scans.

ECG 7T

Courtesy of Zina-Mary Manjaly, University Hospital Zurich

4-electrode ECG and breathing belt, Philips 7T Achieva scanner

-

Description: The ECG data for ultra-high field data is typically much -noisier than at 3 Tesla. Therefore, R-wave peaks are frequently missed -by prospective trigger detection and not marked correctly in the -logfile. This example shows how to select typical R-wave-peaks manually -and let the algorithm find the heartbeat events.

+

Description: The ECG data for ultra-high field data is typically much noisier than at 3 Tesla. Therefore, R-wave peaks are frequently missed by prospective trigger detection and not marked correctly in the logfile. This example shows how to select typical R-wave-peaks manually and let the algorithm find the heartbeat events.

PPU 3T

Courtesy of Diana Wotruba, University and University Hospital of Zurich

-

PPU (finger plethysmograph) and breathing belt, Philips 3T Achieva -scanner

-

Description: Similar to ECG3T, but a plethysmograph instead of an ECG -was used to monitor the cardiac pulsation. Example shows how to extract -heart and breathing rate.

+

PPU (finger plethysmograph) and breathing belt, Philips 3T Achieva scanner

+

Description: Similar to ECG3T, but a plethysmograph instead of an ECG was used to monitor the cardiac pulsation. Example shows how to extract heart and breathing rate.

PPU 7T

Courtesy of Jakob Heinzle and Lars Kasper, TNU, University Zurich and ETH Zurich

-

PPU (finger plethysmograph) and breathing belt, Philips 7T Achieva -scanner

+

PPU (finger plethysmograph) and breathing belt, Philips 7T Achieva scanner

Description: Challenging cardiac data that requires bandpass-filtering during preprocessing, since it is compromised by both high frequency noise (from the scanner, modulated at every slice TR) and low frequency noise (breathing modulation).

Siemens - VB

Siemens has different physiological logfile formats, for which examples are provided here. A detailed description of these formats is on a different wiki page.

-

This is the older Siemens log file format (also available via manual recording), which is part of software release VB, and can be determined by the file extensions .resp, .ecg, .puls, in combination with an optional .dcm DICOM header file for the first (or last) acquired volume.

-

A lot of 7T scanners still use this format.

+

This is the older Siemens log file format (also available via manual recording), which is part of software release _VB_, and can be determined by the file extensions .resp, .ecg, .puls, in combination with an optional .dcm DICOM header file for the first (or last) acquired volume.

+

A lot of 7T scanners still use this format, but it is also the default on modern 3T systems, if you don't have C2P sequences for fMRI (e.g., from CMRR) or WIPs from Siemens (see below).

ECG 3T

Courtesy of Miriam Sebold, Charite Berlin, and Quentin Huys, TNU Zurich

-

4-electrode ECG data, Siemens 3T scanner

+

4-electrode ECG data, Siemens 3T scanner, logfile version 1

Description: Similar to ECG 3T, but acquired on a Siemens system with only one logfile for ECG data. The quality of the signal is challenging, stemming from a patient population.

PPU3T (Sync First and Sync Last)

Courtesy of Alexander Ritter, University of Jena, Germany

-

Siemens 3T pulse oximetry and respiratory bellows data from a complete scan session of a healthy volunteer, plus the DICOM header file of the first and last (382nd) volume of an fMRI run, respectively.

-

This showcases scan timing synchronization using the DICOM timestamps in an intricate case, where the physiological logfile spans the whole scan session (and not only the fMRI run). See TAPAS github issue #55 for further details.

+

Siemens 3T pulse oximetry and respiratory bellows data, logfile version 1 DICOM header file of first and last (382nd) volume of an fMRI run, respectively.

+

Description: This data covering a complete scan session of a healthy volunteer showcases scan timing synchronization using the DICOM timestamps in an intricate case, where the physiological logfile spans the whole scan session (and not only the fMRI run). See TAPAS github issue #55 for further details.

+

ECG 3T - Logversion 3

+

Courtesy of Shahin Safa, see TAPAS GitHub issue 204

+

4-electrode ECG data, Siemens scanner, logfile version 3 corresponding respiratory data: Resp 3T - Logversion 1

+

Description: This is an fMRI study on the auditory system of the brain, which explains the long TR (10 s), to put scanning gaps when presenting the sound to the subject.

+

Resp 3T - Logversion 1

+

Courtesy of Shahin Safa, see TAPAS GitHub issue 204

+

Respiratory bellows data, Siemens scanner, logfile version 1 corresponding cardiac data: ECG 3T - Logversion 3

+

Description: This is an fMRI study on the auditory system of the brain, which explains the long TR (10 s), to put scanning gaps when presenting the sound to the subject.

+

Resp 3T - Logversion 3

+

Courtesy of Lars Kasper, University Health Network Toronto, Canada

+

Respiratory bellows data, Siemens Prisma 3T, logfile version 3

+

Description: Short fingertapping run with logging automatically switched off after about 2 minutes (nominally 5) due to ECG channels not connected, but requested for recording. Biomatrix sensors were not available, but are logged as 4 extra channels with constant values here.

Siemens - HCP

-

The Human Connectome Project uses Siemens scanners, and the logfile format that comes with their published data seems to be pre-converted and custom (even though the documentation desribes the VB format). We have implemented an own reader for that and written a little tutorial for a single subject dataset of the HCP.

+

The Human Connectome Project uses Siemens scanners, and the logfile format that comes with their published data seems to be pre-converted and custom (even though the documentation desribes the VB format). We have implemented an own reader for that and written a little tutorial for a single subject dataset of the HCP.

https://github.com/translationalneuromodeling/tapas/issues/6#issuecomment-361001716

If you download the whole dataset (including functional image files), this example with the additional batches mentioned below also demonstrates how to use the toolbox for model assessment using statistical maps (F-contrasts).

HCP (Subject 178748)

You will have to download the dataset from the HCP yourself, we just provide the matlab batches and the physiological logfile tfMRI_MOTOR_LR_Physio_log.txt here.

For consistency with the other example files, the batch files have been renamed compared to the blog entry:

    -
  • batch_preproc.m -> batch_preproc.m
  • -
  • batch_physio.m -> siemens_hcp_ppu3t_spm_job.m
  • -
  • batch_glm.m -> batch_glm.m
  • +
  • batch_preproc.m -> batch_preproc.m
  • +
  • batch_physio.m -> siemens_hcp_ppu3t_spm_job.m
  • +
  • batch_glm.m -> batch_glm.m

If you want to run the preproc and glm batch, place them on the same level as the subject folder 178748 for the downloaded data. The physio-batch shall reside in the same folder as the physiological logfile tfMRI_MOTOR_LR_Physio_log.txt.

Siemens - VD/VE Tics

-

This is the most recent logfile format of Siemens, included in Software releases VD, VE and sometimes referred to as the Tics format, because all time stamps in all files refer to the same reference point (start of the day) and count in the same intervals or "tics" of 2.5 ms from there.

-

You will recognize this file format via the extensions _Info.log (or _AcquisitionInfo.log), _RESP.log, _ECG.log and _PULS.log. Sometimes, it is also written into the DICOM header (.dcm) file of your functional data directly. In this case, use extractCMRRPhysio.m to convert it to the above separate files before using PhysIO.

+

This is the most recent logfile format of Siemens, included in Software releases _VD_, _VE_ and sometimes referred to as the Tics format, because all time stamps in all files refer to the same reference point (start of the day) and count in the same intervals or "tics" of 2.5 ms from there.

+

You will recognize this file format via the extensions \\\_Info.log (or \\\_AcquisitionInfo.log), \\\_RESP.log, \\\_ECG.log and \\\_PULS.log. Sometimes, it is also written into the DICOM header (.dcm) file of your functional data directly. In this case, use extractCMRRPhysio.m to convert it to the above separate files before using PhysIO.

Most modern Siemens scanners, such as the Prisma or 7T Terra, use this format.

-

There are a couple of variants for this format around (e.g., with the WIP Multiband Protocol that is distributed to multiple sites), and PhysIO tries to support all of them.

+

There are a couple of variants for this format around (e.g., with the WIP Multiband Protocol that is distributed to multiple sites), and PhysIO tries to support all of them.

PPU 3T

Courtesy of Saskia Bollmann, Centre for Advanced Imaging, University of Queensland, Brisbane, Australia

Pulse oximetry and breathing belt data, Siemens Prisma 3T, logfile version EJA_1, multi-echo fMRI (3 echoes)

@@ -1624,45 +1613,41 @@

PPU 3T

PPU 3T Separate Files

Courtesy of Alexandre Sayal CIBIT, University of Coimbra

PPU (finger plethysmograph) and breathing belt, Siemens 3T scanner, Multiband CMRR sequence

-

Description: Raw data that was used to convert to two separate BIDS files above (BIDS/PPU3T_Separate_Files) for cardiac and respiratory recordings, because of differing sampling rate (5 vs 20 ms).

+

Description: Raw data that was used to convert to two separate BIDS files above (BIDS/PPU3T_Separate_Files) for cardiac and respiratory recordings, because of differing sampling rate (5 vs 20 ms).

The UUID and date/time stamps were altered for anonymization.

Technical Documentation: Read-in

Brain Imaging Data Structure (BIDS)

PhysIO supports physiological logfiles prepared according to the BIDS standard

    -
  • In brief, BIDS files are (optionally compressed) tab-separated values -(*.tsv[.gz]) files that contain raw traces of peripheral recordings from -cardiac and respiratory sources, as well as scan trigger events
  • -
  • The header of the columns of this *.tsv file, as well as meta-information, -such as sampling rate and relative onset of physiological logging to MRI scan -onset is described in an accompanying *.json file
      -
    • It is assumed to have the that this *.json file has the same name -(apart from the extension) as the *.tsv file
    • -
    • If PhysIO does not find this file, you can manually enter the timing -information in the log_files structure, and a default column order of -(cardiac, respiratory, trigger) is assumed
    • -
    -
  • -
  • Example *.tsv file (with cardiac, respiratory, trigger column -:
                               -0.949402 -0.00610382 0
    -                           -0.949402 -0.00610382 0
    -                           -0.951233 -0.00915558 0
    -                           -0.951233 -0.00915558 0
    -                           -0.953064 -0.0122073  0
    -                           -0.953064 -0.0122073  0
    -                           -0.95459  -0.0076297  1
    -                           -0.95459  -0.0076297  0
  • -
  • Example *.json file:
    {
    -  "SamplingFrequency": 50.0, 
    -  "Columns": [
    -      "cardiac", 
    -      "respiratory", 
    -      "trigger"
    -  ], 
    -  "StartTime": -255.45
    -}
  • -
  • See tapas_physio_read_physlogfiles_bids.m for more details and technical -documentation.
  • +
  • In brief, BIDS files are (optionally compressed) tab-separated values (\*.tsv\[.gz\]) files that contain raw traces of peripheral recordings from cardiac and respiratory sources, as well as scan trigger events
  • +
  • The header of the columns of this \*.tsv file, as well as meta-information, such as sampling rate and relative onset of physiological logging to MRI scan onset is described in an accompanying \*.json file
      +
    • It is assumed to have the that this \*.json file has the same name (apart from the extension) as the \*.tsv file
    • +
    • If PhysIO does not find this file, you can manually enter the timing information in the log_files structure, and a default column order of (cardiac, respiratory, trigger) is assumed
    • +
    +
  • +
  • Example \*.tsv file (with cardiac, respiratory, trigger column :
  • +
+
                             -0.949402 -0.00610382 0
+                             -0.949402 -0.00610382 0
+                             -0.951233 -0.00915558 0
+                             -0.951233 -0.00915558 0
+                             -0.953064 -0.0122073  0
+                             -0.953064 -0.0122073  0
+                             -0.95459  -0.0076297  1
+                             -0.95459  -0.0076297  0
    +
  • Example \*.json file:
  • +
+
{
+    "SamplingFrequency": 50.0, 
+    "Columns": [
+        "cardiac", 
+        "respiratory", 
+        "trigger"
+    ], 
+    "StartTime": -255.45
+}
    +
  • Note that StartTime refers to when the physiological recording started relative to the first scan volume of the fMRI run, which means that typically this value is negative, because one starts the recording before the onset of scan volumes.
  • +
  • See tapas_physio_read_physlogfiles_bids.m for more details and technical documentation.

BioPac

Mat-file Export (.mat)

@@ -1684,13 +1669,12 @@

Single Text File Export (.txt
  • Export your traces from cardiac and breathing recording devices into 2 text files and select log_files.vendor = 'Custom'. The format is explained in tapas_physio_new or the help window of the Batch Editor:

    • '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.

      +
    • 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
  • - -
  • You have to specify the sampling intervals for these log files (in seconds), via log_files.sampling_interval, e.g. [0.01 0.02] if you have 10 ms (100 Hz) and 20 ms (50 Hz) sampling intervals (frequencies) for cardiac and respiratory data, respectively
  • You will probably have to change log_files.relative_start_acquisition, if logging of your physiological recording device does not start synchronized to the first fMRI volume.
  • @@ -1698,17 +1682,19 @@

    General Electric (GE)

    • Very similar to custom format
    • One text file each for ECG, pulse oximetry and respiratory data, e.g., ECGData_epiRT_phys_0921201215_38_08 or RespData_epiRT_phys_0921201215_38_08
    • -
    • One amplitude entry per line, e.g.,
      2626
      -2649
      -2673
      -2699
      -2727
      -2755
    • +
    • One amplitude entry per line, e.g.,
    • +
    +
     2626
    + 2649
    + 2673
    + 2699
    + 2727
    + 2755
    • sampling rate is determined as a setting beforehand, has to be noted manually (not in log file)

    Philips

      -
    • Physiology automatically recorded into SCANPHYSLOG_<Date>_<Time>.log (one file per scan) as soon as ECG is connected to scanner, and scan is started
    • +
    • Physiology automatically recorded into SCANPHYSLOG\_<Date>\_<Time>.log (one file per scan) as soon as ECG is connected to scanner, and scan is started
    • tabular text (ascii) format, different columns for ECG, pulse oximetry and breathing data
      • additionally, trigger events and gradient timecourses are logged, and can be used for synchronization by the toolbox
      @@ -1729,8 +1715,10 @@

      Philips

      Siemens

      Manual Recording

      Physiological data collection on the Siemens scanners uses the physiological monitoring unit (PMU). The initial sampling is performed at 400 Hz, but through the PMU buffer the effective sampling intervals are ECG: 2.5 ms, RESP: 20 ms, PULS: 20 ms and EXT: 5 ms.

      -

      There are several ways to control the physiological data collection. The 'manual' version is available on all platforms. It uses the telnet mpcu/ideacmdtool to manually start and stop the log file acquisition. The log files (logFileName.ecg, logFileName.resp, logFileName.puls, logFileName.ext) are stored in \MedCom\log. More details on how to record these data can be found here or in the "Other Miscellaneous Topics" slides from the IDEA course.

      -

      An example of a .puls logfile is given below. The data are stored in one long line. The text between 5002 and 6002 forms the header, and the text between 5003 and 6003 the footer. Important information in the footer is the LogStartMDHTime and the LogStopMDHTime (in ms since midnight), which can be used to synchronize the logfiles with the dicom images using the AcquisitionTime in the dicom header (in hhmmss.ms). The values 5000 and 6000 are inserted into the signal trace and indicate trigger events. Note that only the modality which is selected to be displayed during the acquisition will have triggers.

      +

      There are several ways to control the physiological data collection. The 'manual' version is available on all platforms. It uses the telnet mpcu/ideacmdtool to manually start and stop the log file acquisition. The log files (logFileName.ecg, logFileName.resp, logFileName.puls, logFileName.ext) are stored in \\MedCom\\log. More details on how to record these data can be found here or in the "Other Miscellaneous Topics" slides from the IDEA course.

      +

      General Properties

      +

      An example of a .puls logfile is given below. The data are stored in one long line. The text between 5002 and 6002 forms the header, and the text between 5003 and 6003 the footer. Important information in the footer is the LogStartMDHTime and the LogStopMDHTime (in ms since midnight), which can be used to synchronize the logfiles with the DICOM images using the AcquisitionTime in the DICOM header (in hhmmss.ms). The values 5000 and 6000 are inserted into the signal trace and indicate trigger events. Note that only the modality which is selected to be displayed during the acquisition will have triggers.

      +

      We use the time stamp of the clock of the Measurement Data Header (MDH), i.e., computer that controls the scanner, to synchronize with the DICOMs, because this computer also controls the creation of the scan data, i.e., reconstructed DICOM images. This is in accordance to other packages reading Siemens physiological logfile data, e.g., Chris Rorden's PART, with a detailed explanation on the DICOM timestamp in AcquisitionTime found here.

      1 2 40 280 5002 Logging PULSE signal: reduction factor = 1, PULS_SAMPLES_PER_SECOND = 50; PULS_SAMPLE_INTERVAL = 20000 6002 1653 1593 1545 1510 1484 ...
       ACQ FINISHED
        6002 3093 3096 3064 5000 3016 2926 5003
      @@ -1747,9 +1735,25 @@ 

      Manual Recording

      LogStopMDHTime: 47654452 LogStartMPCUTime: 47030087 LogStopMPCUTime: 47652240 -6003

      CMRR Sequence

      +6003

      Updates Logversion 3

      +

      The above specification for physiological Siemens logfiles has been changed by the new logfile version 3 for both cardiac and respiratory data:

      +
        +
      1. There may be multiple info regions between 5002 and 6002 tags, interleaved with the physiological trace, e.g., to log an update in the respiratory cushion gain.
      2. +
      3. The .ecg file has now 4 instead of 2 channels.
      4. +
      5. The .resp file now contains not only the respiratory bellows signal, but also 4 datapoints per sampling interval of the biomatrix signals (integrated sensors in patient bed for breathing detection, e.g., in Siemens Vida series).
      6. +
      +

      A logfile in version 3 therefore adheres to the following template (note the logfile version indicated by keyword LOGVERSION at the start):

      +
      <Header> 5002 <LOGVERSION XX> 6002
      +<[optional] training trace data> 5002 uiHwRevisionPeru ... [optional] 6002
      +5002 <infoRegion3> 6002 5002 <infoRegion4> 6002 ... 5002 <infoRegionN> 6002
      +<trace data 1 (all channels, arbitrary number of samples, trigger markers 5000, 6000)> ...
      +5002 <infoRegionN+1> 6002
      +<trace data 2 (all channels, arbitrary number of samples, trigger markers 5000, 6000)> ...
      +5002 <infoRegionN+2> 6002 ...
      +<trace data M> ... 5003

      PhysIO automatically reads the Siemens logfiles correctly for all modalities and logfile versions 1,2 and 3.

      +

      CMRR Sequence / WIP Advanced Physio Logging

      The CMRR sequence on VD/VE also allows the automatic recording of physiological log files (to be selected in the sequence special card). For more information have a look at the manual. The physiological traces are stored in logFileName_PULS.log, logFileName_RESP, logFileName_ECG.log. Timing information is stored in logFileName_Info.log and external trigger events in logFileName_EXT.log.

      -

      An example of the current format (December 2017, Release 016a) for the logFileName_Info.log is given below:

      +

      An example of the current format (December 2017, Release 016a) for the logFileName_Info.log is given below:

      UUID        = 7a16ea95-ac36-4ee3-9b76-bbb686ac07ca
       ScanDate    = 20171206_150609
       LogVersion  = EJA_1
      @@ -1798,7 +1802,7 @@ 

      Manual Recording

      21747877 PULS 1962

      PhysIO uses the logFileName_Info.log to synchronize the physiological traces with the data acquisition. Note that the reference slice does not yet take into account the multiband slice ordering, but just assumes an even distribution. Older version of the CMRR sequence produced slightly different output files, which might work. Please log an issue if you have a very different format that is not supported.

      Human Connectome Project

      Disclaimer: Most of the information below is a best guess from the developers, but without any guarantee of accuracy.

      -

      The physiological log file (*paradigm*_Physio_log.txt) distributed with the Human Connectome Project data contains respiratory and puls-oximeter data in one file. The first column marks when data acquisition is performed, the second and third contain the respiratory and puls-oximeter traces, respectively. The files are written at a sampling rate of 400Hz and start and end with the scan. PhysIO does provide a reader, you just need to select the appropriate option in the file format tab. An example is provided below:

      +

      The physiological log file (\*paradigm\*\_Physio_log.txt) distributed with the Human Connectome Project data contains respiratory and puls-oximeter data in one file. The first column marks when data acquisition is performed, the second and third contain the respiratory and puls-oximeter traces, respectively. The files are written at a sampling rate of 400Hz and start and end with the scan. PhysIO does provide a reader, you just need to select the appropriate option in the file format tab. An example is provided below:

      1    1904    1756
       1    1904    1756
       1    1907    1754
      @@ -1826,8 +1830,27 @@ 

      Human Connectome Project

      1 1904 1780

      Version History (Changelog)

      RELEASE INFORMATION

      Current Release

      -

      Current version: PhysIO Toolbox Release R2022a, v8.1.0

      -

      April 5th, 2022

      +

      Current version: PhysIO Toolbox Release R2022b, v8.2.0

      +

      September 12th, 2022

      +

      Minor Release Notes (v8.2.0)

      +

      Added

      +
        +
      • Interface tapas_physio_test to TAPAS-generic tapas_test function
      • +
      • Added suport for logfile version 3 of Siemens physio recordings
          +
        • multi ECG/Resp channels and interleaved status messages
        • +
        • new integration test for Siemens VB Logversion 3

          Fixed

          +
        • +
        +
      • +
      • Removed dependence on nanmean (Statistics Toolbox)
          +
        • See GitHub issue #205
        • +
        +
      • +
      • Compatibility with multiple SPM toolbox locations for lmod (GitHUb issue #211)
          +
        • as listed in spm_get_defaults('tbx')
        • +
        +
      • +

      Minor Release Notes (v8.1.0)

      Added

        diff --git a/PhysIO/docs/documentation.pdf b/PhysIO/docs/documentation.pdf index d1f151e1..4e213c0c 100644 Binary files a/PhysIO/docs/documentation.pdf and b/PhysIO/docs/documentation.pdf differ diff --git a/PhysIO/tests/integration/tapas_physio_examples_test.m b/PhysIO/tests/integration/tapas_physio_examples_test.m index e2c1aec4..0a445f96 100644 --- a/PhysIO/tests/integration/tapas_physio_examples_test.m +++ b/PhysIO/tests/integration/tapas_physio_examples_test.m @@ -41,10 +41,9 @@ % path to examples, needed for all test cases function setupOnce(testCase) -% run GE example and extract physio +% Get PhysIO public repo base folder from this file's location testCase.TestData.pathPhysioPublic = fullfile(fileparts(mfilename('fullpath')), '..', '..'); -% TODO: Make generic! -testCase.TestData.pathExamples = fullfile(testCase.TestData.pathPhysioPublic, '..', 'examples'); +testCase.TestData.pathExamples = tapas_physio_get_path_examples(testCase.TestData.pathPhysioPublic); end @@ -82,6 +81,7 @@ function test_bids_ppu3t_separate_files_matlab_only(testCase) run_example_and_compare_reference(testCase, dirExample, doUseSpm) end + function test_biopac_txt_ppu3t_matlab_only(testCase) %% Compares previously saved physio-structure and multiple regressors file % to current output of re-run of BioPac_txt/PPU3T example using matlab only @@ -108,7 +108,6 @@ function test_philips_ecg3t_matlab_only(testCase) run_example_and_compare_reference(testCase, dirExample, doUseSpm) end - function test_philips_ecg3t_v2_matlab_only(testCase) %% Compares previously saved physio-structure and multiple regressors file % to current output of re-run of Philips/ECG3T_V2 example using matlab only @@ -117,7 +116,6 @@ function test_philips_ecg3t_v2_matlab_only(testCase) run_example_and_compare_reference(testCase, dirExample, doUseSpm) end - function test_philips_ecg7t_matlab_only(testCase) %% Compares previously saved physio-structure and multiple regressors file % to current output of re-run of Philips/ECG7T example using matlab only @@ -126,7 +124,6 @@ function test_philips_ecg7t_matlab_only(testCase) run_example_and_compare_reference(testCase, dirExample, doUseSpm) end - function test_philips_ppu3t_matlab_only(testCase) %% Compares previously saved physio-structure and multiple regressors file % to current output of re-run of Philips/PPU3T example using matlab only @@ -153,6 +150,13 @@ function test_siemens_vb_ecg3t_matlab_only(testCase) run_example_and_compare_reference(testCase, dirExample, doUseSpm) end +function test_siemens_vb_ecg3t_logversion_3_matlab_only(testCase) +%% Compares previously saved physio-structure and multiple regressors file +% to current output of re-run of Siemens_VB/ECG3T example using matlab only +dirExample = 'Siemens_VB/ECG3T_Logversion_3'; +doUseSpm = false; +run_example_and_compare_reference(testCase, dirExample, doUseSpm) +end function test_siemens_vb_ppu3t_sync_first_matlab_only(testCase) %% Compares previously saved physio-structure and multiple regressors file @@ -163,7 +167,6 @@ function test_siemens_vb_ppu3t_sync_first_matlab_only(testCase) run_example_and_compare_reference(testCase, dirExample, doUseSpm) end - function test_siemens_vb_ppu3t_sync_last_matlab_only(testCase) %% Compares previously saved physio-structure and multiple regressors file % to current output of re-run of Siemens_VB/PPU3T example using matlab only @@ -173,6 +176,22 @@ function test_siemens_vb_ppu3t_sync_last_matlab_only(testCase) run_example_and_compare_reference(testCase, dirExample, doUseSpm) end +function test_siemens_vb_resp3t_logversion_1_matlab_only(testCase) +%% Compares previously saved physio-structure and multiple regressors file +% to current output of re-run of Siemens_VB/ECG3T example using matlab only +dirExample = 'Siemens_VB/RESP3T_Logversion_1'; +doUseSpm = false; +run_example_and_compare_reference(testCase, dirExample, doUseSpm) +end + +function test_siemens_vb_resp3t_logversion_3_matlab_only(testCase) +%% Compares previously saved physio-structure and multiple regressors file +% to current output of re-run of Siemens_VB/ECG3T example using matlab only +dirExample = 'Siemens_VB/RESP3T_Logversion_3'; +doUseSpm = false; +run_example_and_compare_reference(testCase, dirExample, doUseSpm) +end + function test_siemens_vd_ppu3t_matlab_only(testCase) %% Compares previously saved physio-structure and multiple regressors file @@ -232,6 +251,7 @@ function test_bids_ppu3t_separate_files_with_spm(testCase) run_example_and_compare_reference(testCase, dirExample, doUseSpm) end + function test_biopac_txt_ppu3t_with_spm(testCase) %% Compares previously saved physio-structure and multiple regressors file % to current output of re-run of BioPac_txt/PPU3T example using SPM Batch Editor @@ -258,7 +278,6 @@ function test_philips_ecg3t_with_spm(testCase) run_example_and_compare_reference(testCase, dirExample, doUseSpm) end - function test_philips_ecg3t_v2_with_spm(testCase) %% Compares previously saved physio-structure and multiple regressors file % to current output of re-run of Philips/ECG3T_V2 example using SPM Batch Editor @@ -267,7 +286,6 @@ function test_philips_ecg3t_v2_with_spm(testCase) run_example_and_compare_reference(testCase, dirExample, doUseSpm) end - function test_philips_ecg7t_with_spm(testCase) %% Compares previously saved physio-structure and multiple regressors file % to current output of re-run of Philips/ECG7T example using SPM Batch Editor @@ -276,7 +294,6 @@ function test_philips_ecg7t_with_spm(testCase) run_example_and_compare_reference(testCase, dirExample, doUseSpm) end - function test_philips_ppu3t_with_spm(testCase) %% Compares previously saved physio-structure and multiple regressors file % to current output of re-run of Philips/PPU3T example using SPM Batch Editor @@ -303,6 +320,14 @@ function test_siemens_vb_ecg3t_with_spm(testCase) run_example_and_compare_reference(testCase, dirExample, doUseSpm) end +function test_siemens_vb_ecg3t_logversion_3_with_spm(testCase) +%% Compares previously saved physio-structure and multiple regressors file +% to current output of re-run of Siemens_VB/ECG3T example using SPM Batch Editor +dirExample = 'Siemens_VB/ECG3T_Logversion_3'; +doUseSpm = true; +run_example_and_compare_reference(testCase, dirExample, doUseSpm) +end + function test_siemens_vb_ppu3t_sync_first_with_spm(testCase) %% Compares previously saved physio-structure and multiple regressors file % to current output of re-run of Siemens_VB/PPU3T example @@ -321,6 +346,22 @@ function test_siemens_vb_ppu3t_sync_last_with_spm(testCase) run_example_and_compare_reference(testCase, dirExample, doUseSpm) end +function test_siemens_vb_resp3t_logversion_1_with_spm(testCase) +%% Compares previously saved physio-structure and multiple regressors file +% to current output of re-run of Siemens_VB/ECG3T example using SPM Batch Editor +dirExample = 'Siemens_VB/RESP3T_Logversion_1'; +doUseSpm = true; +run_example_and_compare_reference(testCase, dirExample, doUseSpm) +end + +function test_siemens_vb_resp3t_logversion_3_with_spm(testCase) +%% Compares previously saved physio-structure and multiple regressors file +% to current output of re-run of Siemens_VB/ECG3T example using SPM Batch Editor +dirExample = 'Siemens_VB/RESP3T_Logversion_3'; +doUseSpm = true; +run_example_and_compare_reference(testCase, dirExample, doUseSpm) +end + function test_siemens_vd_ppu3t_with_spm(testCase) %% Compares previously saved physio-structure and multiple regressors file @@ -330,7 +371,6 @@ function test_siemens_vd_ppu3t_with_spm(testCase) run_example_and_compare_reference(testCase, dirExample, doUseSpm) end - function test_siemens_vd_ppu3t_for_bids_with_spm(testCase) %% Compares previously saved physio-structure and multiple regressors file % to current output of re-run of Siemens_VD/PPU3T_For_BIDS example using SPM Batch Editor @@ -398,6 +438,7 @@ function run_example_and_compare_reference(testCase, dirExample, doUseSpm, ... % hard-coded relative tolerance relTol = 0.01; % 0.01 means 1 percent deviation from expected value allowed +absTol = 1e-6; % for time courses (e.g., breathing) that reach close to 0, relative tolerance can be misleading, use relative value to max instead %% Generic settings % methods for recursively comparing structures, see @@ -487,8 +528,21 @@ function run_example_and_compare_reference(testCase, dirExample, doUseSpm, ... IsEqualTo(expPhysio.ons_secs.raw, ... 'Using', StructComparator(NumericComparator, 'Recursively', true), ... 'Within', RelativeTolerance(relTol), ... - 'IgnoringFields', {'spulse_per_vol'}... + 'IgnoringFields', {'spulse_per_vol', 'fr'}... ), 'Comparing all numeric subfields of ons_secs.raw to check read-in and basic filtering of phys recordings'); + + % test filtered respiratory trace separetely, because of values close + % to zero in trace (relative errors can be huge even if 1e-6) + % if isempty, set arbitrary tolerance + absTolFr = absTol*max(abs(expPhysio.ons_secs.raw.fr)); + if isempty(absTolFr) + absTolFr = absTol; + end + + verifyEqual(testCase, actPhysio.ons_secs.raw.fr, expPhysio.ons_secs.raw.fr, ... + 'AbsTol', absTolFr, ... + 'Comparing ons_secs.raw.fr (raw filtered respiratory trace)'); + end % 2. Check some crucial timing parameters more vigorously diff --git a/PhysIO/tests/tapas_physio_get_path_examples.m b/PhysIO/tests/tapas_physio_get_path_examples.m new file mode 100644 index 00000000..5b6aac1e --- /dev/null +++ b/PhysIO/tests/tapas_physio_get_path_examples.m @@ -0,0 +1,45 @@ +function pathExamples = tapas_physio_get_path_examples(pathPhysioPublic) +% Returns GitLab-internal or TAPAS-public PhysIO Examples folder, based on +% location of public PhysIO directory +% +% pathExamples = tapas_physio_get_path_examples(pathPhysioPublic) +% +% IN +% pathPhysioPublic location of public PhysIO folder, e.g., +% 'tapas/PhysIO' +% OUT +% +% EXAMPLE +% tapas_physio_get_path_examples('tapas/PhysIO') +% +% See also + +% Author: Lars Kasper +% Created: 2022-09-05 +% Copyright (C) 2022 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 . + +if nargin < 1 + pathPhysioPublic = fullfile(fileparts(mfilename('fullpath')), '..'); +end + +% try PhysIO-internal GitLab examples first +pathExamples = fullfile(pathPhysioPublic, '..', 'examples'); + +% otherwise use public TAPAS examples +if ~isfolder(fullfile(pathExamples, 'BIDS')) + pathExamples = fullfile(pathPhysioPublic, ... + '..', 'examples', tapas_get_current_version(), 'PhysIO'); +end + +% If no examples folder found, suggest to download them via tapas-function +if ~isfolder(fullfile(pathExamples, 'BIDS')) + physio = tapas_physio_new(); + tapas_physio_log('No PhysIO examples data found. Please download via tapas_download_example_data()', physio.verbose, 2); +end diff --git a/PhysIO/tests/tapas_physio_test.m b/PhysIO/tests/tapas_physio_test.m new file mode 100644 index 00000000..e5a30ffa --- /dev/null +++ b/PhysIO/tests/tapas_physio_test.m @@ -0,0 +1,81 @@ +function [nTestFailed, nTestTotal, testResults] = tapas_physio_test(level) +% Interface of PhysIO-Toolbox-specific test functions to global tapas_test +% +% [nTestFailed, nTestTotal] = tapas_physio_test(level) +% +% IN +% level Depth of testing required +% Approximate run time per level +% 0: around a minute +% here: unit tests only +% 1: around 5 minutes (a coffee) +% here: matlab-only integration tests, no SPM integration +% 2: around an hour (time for lunch) +% here: matlab AND SPM integration tests, takes less +% than 10 minutes +% 3: overnight (time to freakout [deadline]) +% OUT +% nTestFailed number of failed tests +% nTestTotal total number of executed tests +% testResults detailed test results object (matlab.unittest.TestResult) +% which may be queried to retrieve error messages or +% rerun tests +% +% EXAMPLE +% tapas_physio_test(0) +% +% See also tapas_test matlab.unittest.TestResult + +% Author: Lars Kasper +% Created: 2022-09-04 +% Copyright (C) 2022 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 . + +if nargin < 1 + level = 0; +end + +tic + +testResults = []; + +% Returns an error message, if no example data found +tapas_physio_get_path_examples(); + +if level >= 0 + testResults = [testResults tapas_physio_run_unit_tests()]; +else + nTestTotal = 0; + nTestFailed = 0; + return +end + +if level == 1 + % code adapted from run_integration_tests, but chooses only tests with + % 'matlab' in the name (SPM GUI tests have SPM in the name) + import matlab.unittest.TestSuite; + + pathTests = fullfile(fileparts(mfilename('fullpath')), 'integration'); + suiteFolder = TestSuite.fromFolder(pathTests, ... + 'IncludingSubfolders', true, 'Name', '*matlab_only*'); + testResults = [testResults run(suiteFolder)]; +end + +if level >= 2 + testResults = [testResults tapas_physio_run_integration_tests()]; +end + +nTestTotal = numel(testResults); +nTestFailed = sum([testResults.Failed]); + +% pretty summary output +fprintf('\n\n\n\tTable of all executed PhysIO Tests:\n\n'); +disp(testResults.table); + +toc diff --git a/PhysIO/tests/unit/preproc/tapas_physio_filter_cardiac_test.m b/PhysIO/tests/unit/preproc/tapas_physio_filter_cardiac_test.m index 4d932985..c7f6794f 100644 --- a/PhysIO/tests/unit/preproc/tapas_physio_filter_cardiac_test.m +++ b/PhysIO/tests/unit/preproc/tapas_physio_filter_cardiac_test.m @@ -26,15 +26,26 @@ tests = functiontests(localfunctions); end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Auxiliary functions for automation of code folder structure +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% path to examples, needed for all test cases +function setupOnce(testCase) +% Get PhysIO public repo base folder from this file's location +testCase.TestData.pathPhysioPublic = fullfile(fileparts(mfilename('fullpath')), '..', '..', '..'); +testCase.TestData.pathExamples = tapas_physio_get_path_examples(testCase.TestData.pathPhysioPublic); +% for time courses (e.g., breathing) that reach close to 0, relative +% tolerance can be misleading, use relative value to max instead +testCase.TestData.absTol = 1e-6; +end + + function test_philips_ppu7t_filter_cheby2(testCase) %% Compares previously saved Chebychev Type 2 IIR-filtered cropped cardiac % time course with current re-run of same batch from Philips 7T PPU data -% run GE example and extract physio -pathPhysioPublic = fullfile(fileparts(mfilename('fullpath')), '..', '..', '..'); -pathExamples = fullfile(pathPhysioPublic, '..', 'examples'); - -pathCurrentExample = fullfile(pathExamples, 'Philips/PPU7T'); +pathCurrentExample = fullfile(testCase.TestData.pathExamples, 'Philips/PPU7T'); cd(pathCurrentExample); % for prepending absolute paths correctly fileExample = fullfile(pathCurrentExample, 'philips_ppu7t_spm_job.m'); run(fileExample); % retrieve matlabbatch @@ -53,7 +64,7 @@ function test_philips_ppu7t_filter_cheby2(testCase) actPhysio = tapas_physio_main_create_regressors(physio); % load physio from reference data -fileReferenceData = fullfile(pathExamples, 'TestReferenceResults', 'preproc', ... +fileReferenceData = fullfile(testCase.TestData.pathExamples, 'TestReferenceResults', 'preproc', ... 'physio_filter_cardiac_cheby2.mat'); load(fileReferenceData, 'physio'); expPhysio = physio; @@ -62,19 +73,17 @@ function test_philips_ppu7t_filter_cheby2(testCase) actSolution = actPhysio.ons_secs.c; expSolution = expPhysio.ons_secs.c; -verifyEqual(testCase, actSolution, expSolution); +% RelTol would be too conservative, because values close to 0 in raw +% timeseries +verifyEqual(testCase, actSolution, expSolution, ... + 'AbsTol', testCase.TestData.absTol*max(abs(expSolution))); end function test_philips_ppu7t_filter_butter(testCase) %% Compares previously saved butterworth-filtered cropped cardiac time course % with current re-run of same batch from Philips 7T PPU data -% run GE example and extract physio -pathPhysioPublic = fullfile(fileparts(mfilename('fullpath')), '..', '..', '..'); -% TODO: Make generic! -pathExamples = fullfile(pathPhysioPublic, '..', 'examples'); - -pathCurrentExample = fullfile(pathExamples, 'Philips/PPU7T'); +pathCurrentExample = fullfile(testCase.TestData.pathExamples, 'Philips/PPU7T'); cd(pathCurrentExample); % for prepending absolute paths correctly fileExample = fullfile(pathCurrentExample, 'philips_ppu7t_spm_job.m'); run(fileExample); % retrieve matlabbatch @@ -94,7 +103,7 @@ function test_philips_ppu7t_filter_butter(testCase) actPhysio = tapas_physio_main_create_regressors(physio); % load physio from reference data -fileReferenceData = fullfile(pathExamples, 'TestReferenceResults', 'preproc', ... +fileReferenceData = fullfile(testCase.TestData.pathExamples, 'TestReferenceResults', 'preproc', ... 'physio_filter_cardiac_butter.mat'); load(fileReferenceData, 'physio'); expPhysio = physio; @@ -103,5 +112,8 @@ function test_philips_ppu7t_filter_butter(testCase) actSolution = actPhysio.ons_secs.c; expSolution = expPhysio.ons_secs.c; -verifyEqual(testCase, actSolution, expSolution); +% RelTol would be too conservative, because values close to 0 in raw +% timeseries +verifyEqual(testCase, actSolution, expSolution, ... + 'AbsTol', testCase.TestData.absTol*max(abs(expSolution))); end diff --git a/PhysIO/tests/unit/readin/tapas_physio_readin_bids_test.m b/PhysIO/tests/unit/readin/tapas_physio_readin_bids_test.m index 4688a394..3ee96ca2 100644 --- a/PhysIO/tests/unit/readin/tapas_physio_readin_bids_test.m +++ b/PhysIO/tests/unit/readin/tapas_physio_readin_bids_test.m @@ -31,10 +31,10 @@ % results function test_readin_bids_ppu3t(testCase) -% run BIDS PPU example and extract physio +% Get PhysIO public repo base folder from this file's location pathPhysioPublic = fullfile(fileparts(mfilename('fullpath')), '..', '..', '..'); -% TODO: Make generic! -pathExamples = fullfile(pathPhysioPublic, '..', 'examples'); +pathExamples = tapas_physio_get_path_examples(pathPhysioPublic); + % load SPM matlabbatch, but convert to pure script before executing % remove unnecessary (beyond read-in) part from job exeuction (e.g. @@ -79,8 +79,7 @@ function test_readin_bids_cpulse3t(testCase) % run BIDS cpulse3t example and extract physio pathPhysioPublic = fullfile(fileparts(mfilename('fullpath')), '..', '..', '..'); -% TODO: Make generic! -pathExamples = fullfile(pathPhysioPublic, '..', 'examples'); +pathExamples = tapas_physio_get_path_examples(pathPhysioPublic); % load SPM matlabbatch, but convert to pure script before executing % remove unnecessary (beyond read-in) part from job exeuction (e.g. diff --git a/PhysIO/tests/unit/utils/tapas_physio_findpeaks_test.m b/PhysIO/tests/unit/utils/tapas_physio_findpeaks_test.m index 826bffa2..889f3889 100644 --- a/PhysIO/tests/unit/utils/tapas_physio_findpeaks_test.m +++ b/PhysIO/tests/unit/utils/tapas_physio_findpeaks_test.m @@ -41,8 +41,7 @@ function test_ge_ppu3t_peaks(testCase) % run GE example and extract physio pathPhysioPublic = fullfile(fileparts(mfilename('fullpath')), '..', '..', '..'); -% TODO: Make generic! -pathExamples = fullfile(pathPhysioPublic, '..', 'examples'); +pathExamples = tapas_physio_get_path_examples(pathPhysioPublic); if doUseSpm pathCurrentExample = fullfile(pathExamples, 'GE/PPU3T'); diff --git a/PhysIO/tests/unit/utils/tapas_physio_pca_test.m b/PhysIO/tests/unit/utils/tapas_physio_pca_test.m index 45414df5..2a519084 100644 --- a/PhysIO/tests/unit/utils/tapas_physio_pca_test.m +++ b/PhysIO/tests/unit/utils/tapas_physio_pca_test.m @@ -46,6 +46,9 @@ function test_pca_more_voxels_than_volumes(testCase) end % function function test_pca_more_volumes_than_voxels(testCase) +% NOTE: This test is trivially true for now, because tapas_physio_pca falls +% back on the stats-pca implementation for nVolumes >= nVoxels + % Compare the outputs of tapas_physio_pca(X,'svd') and % tapas_physio_pca(X,'stats-pca') and check is they are close enough % numericaly, within a certain tolerence. @@ -58,8 +61,15 @@ function test_pca_more_volumes_than_voxels(testCase) NVolumes = 300; timeseries = randn(nVoxels,NVolumes); -% Perform PCA with both methods -verifyError(testCase,@() tapas_physio_pca( timeseries, 'svd' ) ,'tapas_physio_pca:NotEnoughVoxels') -verifyError(testCase,@() tapas_physio_pca( timeseries, 'stats-pca' ) ,'tapas_physio_pca:NotEnoughVoxels') +[svd. COEFF, svd. SCORE, svd. LATENT, svd. EXPLAINED, svd. MU] = tapas_physio_pca( timeseries, 'svd' ); +[stats.COEFF, stats.SCORE, stats.LATENT, stats.EXPLAINED, stats.MU] = tapas_physio_pca( timeseries, 'stats-pca' ); + +% Compare both methods +verifyEqual(testCase,svd.COEFF ,stats.COEFF ) +verifyEqual(testCase,svd.SCORE ,stats.SCORE ) +verifyEqual(testCase,svd.LATENT ,stats.LATENT ) +verifyEqual(testCase,svd.EXPLAINED,stats.EXPLAINED) +verifyEqual(testCase,svd.MU ,stats.MU ) + end % function diff --git a/README.md b/README.md index 4f6fc2ce..762325c1 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ ![TAPAS Logo](misc/TapasLogo.png?raw=true "TAPAS Logo") -*Version 6.0.1* +*Version 6.1.0* T A P A S - Translational Algorithms for Psychiatry-Advancing Science. ======================================================================== diff --git a/misc/log_tapas.txt b/misc/log_tapas.txt index d80976b3..deaea13d 100644 --- a/misc/log_tapas.txt +++ b/misc/log_tapas.txt @@ -1,3 +1,4 @@ +6.1.0 https://www.tapas.tnu-zurich.com/examples_v6.1.0.zip ab01a2082cee20b5bb576b48d7f4e66c 6.0.1 https://www.tapas.tnu-zurich.com/examples_v6.0.0.zip b4bf7601b9a521159af00b00019989b6 6.0.0 https://www.tapas.tnu-zurich.com/examples_v6.0.0.zip b4bf7601b9a521159af00b00019989b6 5.1.2 https://www.tapas.tnu-zurich.com/examples_v5.0.0.zip 24c1248d3054f025b853b5ab7ce08a1a diff --git a/misc/tapas_get_toolbox_infos.m b/misc/tapas_get_toolbox_infos.m index 9ea13b3f..8b5e64ba 100644 --- a/misc/tapas_get_toolbox_infos.m +++ b/misc/tapas_get_toolbox_infos.m @@ -41,9 +41,11 @@ infos.physio.init_dir = strcat('PhysIO',filesep,'code'); infos.physio.init_function = 'tapas_physio_init'; - infos.physio.dependencies = []; + infos.physio.dependencies = {'Signal Processing Toolbox', ... + 'Image Processing Toolbox', ... + 'Statistics and Machine Learning Toolbox'}; infos.physio.diagnose_files = 'tapas_physio_main_create_regressors'; - infos.physio.test_function_name = ''; + infos.physio.test_function_name = 'tapas_physio_test'; infos = tapas_default_toolbox_info(infos,'rDCM');