9 function to write out preprocessed traces as bids

.gitignore

@@ -17,7 +17,8 @@ lib
 *.DS_Store
 .idea
 .idea/*
-examples
+physio_out
+.gitignore
 
 # Python
 # See e.g. https://github.com/github/gitignore/blob/master/Python.gitignore

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.

Contributor-License-Agreement.md

@@ -12,7 +12,7 @@ and
 each GitHub User listed in the following table:  
 
 Name                     | Company/Institution/Address        | City      | Country | E-Mail/GitHub Username | CLA version accepted
------------------------- | ---------------------------------- | --------- | ------- | ---------------------- | ------------------
+------------------------ | ---------------------------------- | --------- | ------- | ---------------------- | --------------------
 Lars Kasper              | TNU, University of Zurich          | Zurich    | CH      | mrikasper              | 1.0
 Eduardo Aponte           | TNU, University of Zurich          | Zurich    | CH      | tnutapas               | 1.0
 Daniel Hoffmann Ayala    | Technical University               | Munich    | D       | DanielHoffmannAyala    | 1.0
@@ -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
 

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
 
+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
+- 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)
 ----------------------------

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

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

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,30 +22,36 @@
 % (either version 3 or, at your option, any later version). For further details, see the file
 % COPYING or <http://www.gnu.org/licenses/>.
 
-tapas_physio_strip_fields(dataCardiac);
+tapas_physio_strip_fields(dataPhysio);
 
 stringTitle = 'Read-In: Raw Siemens physlog data';
 fh = tapas_physio_get_default_fig_params();
 set(gcf, 'Name', stringTitle);
 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);

PhysIO/code/preproc/tapas_physio_filter_respiratory.m

@@ -61,9 +61,13 @@
 
 %% 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
 
@@ -178,4 +182,4 @@
     legend(handles, labels);
 end
 
-end
+end

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!

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);