Merge pull request #16 from likeajumprope/development

Update after merge of Physio read-in for Siemens logfile version 3 #14 from development into working branches

PhysIO/CHANGELOG.md

@@ -8,6 +8,15 @@ Current Release
 
 April 5th, 2022
 
+Upcoming Release Notes (v8.2.0-beta)
+------------------------------------
+
+### Added
+- Added suport for logfile version 3 of Siemens physio recordings 
+    - multi ECG/Resp channels and interleaved status messages
+### Fixed
+- Removed dependence on `nanmean` (Statistics Toolbox)
+    - See GitHub issue #205 and 
 
 Minor Release Notes (v8.1.0)
 ----------------------------

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
 

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)

PhysIO/code/readin/tapas_physio_read_physlogfiles_siemens.m

@@ -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
 
 

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;