Skip to content

Commit

Permalink
Merge pull request translationalneuromodeling#167 from alexsayal/patch-1
Browse files Browse the repository at this point in the history 
Allow separate cardiac and resp files in BIDS format
  • Loading branch information
@mrikasper
mrikasper authored Mar 29, 2022
2 parents 78a0193 + 8e389e9 commit 2f07036
Show file tree
Hide file tree
Showing 2 changed files with 290 additions and 14 deletions.
1 change: 1 addition & 0 deletions Contributor-License-Agreement.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@ Daniel Hoffmann Ayala | Technical University | Munich | D
Benoît Béranger | CENIR, ICM | Paris | FR | benoitberanger
Stephan Heunis | Eindhoven University of Technology | Eindhoven | NL | jsheunis
Niklas Bürgi | SNS, University of Zurich | Zurich | CH | nbuergi
Alexandre Sayal | CIBIT, University of Coimbra | Coimbra | PT | alexsayal
**- Add Entry here -** | **- Add Entry here -** | **Add** | **Add** | **Add**

(hereinafter referred to as "Contributor")
Expand Down
303 changes: 289 additions & 14 deletions PhysIO/code/readin/tapas_physio_read_physlogfiles_bids.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,9 +7,9 @@
% column order of physiological recordings will be read from there as well
% as values for sampling_interval and relative_start_acquisition, if they were
% empty before
%
%
% Details of the Brain Imaging Data Structure (BIDS) standard for peripheral
% physiological recordings can be found here:
% physiological recordings can be found here:
%
% https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/
% 06-physiological-and-other-continous-recordings.html
Expand All @@ -18,7 +18,7 @@
% log_files, cardiac_modality, verbose, varargin)
%
% IN log_files
% .log_cardiac *.tsv[.gz] file ([zipped] tab separated file,
% .log_cardiac *.tsv[.gz] file ([zipped] tab separated file,
% contains 3 columns of the form
% cardiac respiratory trigger
% -0.949402 -0.00610382 0
Expand Down Expand Up @@ -80,15 +80,277 @@
hasRespirationFile = ~isempty(log_files.respiration);
hasCardiacFile = ~isempty(log_files.cardiac);
hasExplicitJsonFile = ~isempty(log_files.scan_timing);
diffCardiacRespirationFile = ~strcmp(log_files.cardiac,log_files.respiration);

if hasCardiacFile
if hasCardiacFile && hasRespirationFile && diffCardiacRespirationFile % if the sampling rates of the two signals are different, they are registered in different files according to BIDS.
fileName{1} = log_files.cardiac;
fileName{2} = log_files.respiration;

[c,r,t,cpulse,acq_codes] = tapas_physio_read_physlogfiles_bids_separate(fileName,log_files,hasExplicitJsonFile,DEBUG,verbose);

elseif hasCardiacFile
fileName = log_files.cardiac;
[c,r,t,cpulse,acq_codes] = tapas_physio_read_physlogfiles_bids_unified(fileName,log_files,hasExplicitJsonFile,DEBUG,verbose);

elseif hasRespirationFile
fileName = log_files.respiration;
[c,r,t,cpulse,acq_codes] = tapas_physio_read_physlogfiles_bids_unified(fileName,log_files,hasExplicitJsonFile,DEBUG,verbose);
end

end

%% Function for separate physio files (cardiac and resp with different sampling rates)
function [c,r,t,cpulse,acq_codes] = tapas_physio_read_physlogfiles_bids_separate(fileName,log_files,hasExplicitJsonFile,DEBUG,verbose)

%% Read cardiac
%% check for .gz files and unzip to temporary file
[~, ~, ext] = fileparts(fileName{1});

isZipped = strcmpi(ext, '.gz');

if isZipped
fileJson = regexprep(fileName{1}, '\.tsv\.gz', '\.json');
tempFilePath = tempname; % tempname is matlab inbuilt
fileNameCardiac = gunzip(fileName{1}, tempFilePath);
fileNameCardiac = fileNameCardiac{1};
else
fileJson = regexprep(fileName{1}, '\.tsv', '\.json');
end

if hasExplicitJsonFile
fileJson = log_files.scan_timing;
end

hasJsonFile = isfile(fileJson);

if hasJsonFile
val = jsondecode(fileread(fileJson));
else
verbose = tapas_physio_log(...
['No .json file found. Please specify log_files.sampling_interval' ...
' and log_files.relative_start_acquisition explicitly.'], verbose, 1);
end

dtCardiac = log_files.sampling_interval;
if isempty(dtCardiac)
if hasJsonFile
dtCardiac = 1/val.SamplingFrequency;
else
verbose = tapas_physio_log(...
['No .json file found and empty log_files.sampling_interval. ' ...
'Please specify explicitly.'], verbose, 2);
end
end

% sum implicit (.json) and explicit relative shifts of log/scan acquisition
if isempty(log_files.relative_start_acquisition)
if hasJsonFile
% in BIDS, start of the phys logging is stated relative to the first volume scan start.
% PhysIO defines the scan acquisiton relative to the phys log start
tRelStartScan = -val.StartTime;
else
verbose = tapas_physio_log(...
['No .json file found and empty log_files.relative_start_acquisition. ' ...
'Please specify explicitly.'], verbose, 2);
end
else
if hasJsonFile
% add both delays
tRelStartScan = log_files.relative_start_acquisition - val.StartTime;
else
tRelStartScan = log_files.relative_start_acquisition;
end
end

% default columns in text file for phys recordings; overruled by JSON file
% 1 = cardiac, 2 = resp, 3 = trigger
bidsColumnNamesCardiac = {'cardiac', 'trigger'};
idxCol = 1:2; %set default values for columns from BIDS
for iCol = 1:2
if hasJsonFile
idxCurrCol = find(cellfun(@(x) isequal(lower(x), bidsColumnNamesCardiac{iCol}), val.Columns));
if ~isempty(idxCurrCol)
idxCol(iCol) = idxCurrCol;
end
end
end

C = tapas_physio_read_columnar_textfiles(fileNameCardiac, 'BIDS');
c = double(C{idxCol(1)});
iAcqStart = (double(C{idxCol(2)})~=0); % trigger has 1, rest is 0;

%% delete temporary unzipped file
if isZipped
[status,message,messageId] = rmdir(tempFilePath, 's');
% warning if deletion failed
if status == 0
tapas_physio_log(sprintf('%s: %s', messageId, message), verbose, 1)
end
end

%% Read respiratory
%% check for .gz files and unzip to temporary file
[fp, fn, ext] = fileparts(fileName);
[~, ~, ext] = fileparts(fileName{2});

isZipped = strcmpi(ext, '.gz');

if isZipped
fileJson = regexprep(fileName{2}, '\.tsv\.gz', '\.json');
tempFilePath = tempname; % tempname is matlab inbuilt
fileNameRespiration = gunzip(fileName{2}, tempFilePath);
fileNameRespiration = fileNameRespiration{1};
else
fileJson = regexprep(fileName{2}, '\.tsv', '\.json');
end

if hasExplicitJsonFile
fileJson = log_files.scan_timing;
end

hasJsonFile = isfile(fileJson);

if hasJsonFile
val = jsondecode(fileread(fileJson));
else
verbose = tapas_physio_log(...
['No .json file found. Please specify log_files.sampling_interval' ...
' and log_files.relative_start_acquisition explicitly.'], verbose, 1);
end

dtRespiration = log_files.sampling_interval;
if isempty(dtRespiration)
if hasJsonFile
dtRespiration = 1/val.SamplingFrequency;
else
verbose = tapas_physio_log(...
['No .json file found and empty log_files.sampling_interval. ' ...
'Please specify explicitly.'], verbose, 2);
end
end

% default columns in text file for phys recordings; overruled by JSON file
% 1 = resp, 2 = trigger
bidsColumnNamesRespiration = {'respiratory', 'trigger'};
idxCol = 1:2; %set default values for columns from BIDS
for iCol = 1:2
if hasJsonFile
idxCurrCol = find(cellfun(@(x) isequal(lower(x), bidsColumnNamesRespiration{iCol}), val.Columns));
if ~isempty(idxCurrCol)
idxCol(iCol) = idxCurrCol;
end
end
end

C = tapas_physio_read_columnar_textfiles(fileNameRespiration, 'BIDS');
r = double(C{idxCol(1)});

%% Create timing vector from samples
nSamplesCardiac = length(c);
nSamplesRespiration = length(r);

tCardiac = -tRelStartScan + ((0:(nSamplesCardiac-1))*dtCardiac)';
tRespiration = -tRelStartScan + ((0:(nSamplesRespiration-1))*dtRespiration)';

%% Deal with NaNs in c and r timecourse
c(isnan(c)) = interp1(tCardiac(~isnan(c)), c(~isnan(c)), tCardiac(isnan(c)));
r(isnan(r)) = interp1(tRespiration(~isnan(r)), r(~isnan(r)), tRespiration(isnan(r)));

%% occasionally, cardiac time course is instead containing 0/1 cardiac triggers,
% and not raw trace; check this and populate cpulse accordingly
if all(ismember(unique(c), [1 0]))
cpulse = tCardiac(c==1);
else
cpulse = [];
end

%% interpolate to greater precision, if both files exist and
% 2 different sampling rates are given
% interpolate acq_codes and trace with lower sampling rate to higher
%rate

%dtCardiac = tCardiac(2)-tCardiac(1);
%dtRespiration = tRespiration(2) - tRespiration(1);

isHigherSamplingCardiac = dtCardiac < dtRespiration;
if isHigherSamplingCardiac
t = tCardiac;
rInterp = interp1(tRespiration, r, t);
%racq_codesInterp = interp1(tRespiration, racq_codes, t, 'nearest');
%acq_codes = cacq_codes + racq_codesInterp;

if DEBUG
fh = plot_interpolation(tRespiration, r, t, rInterp, ...
{'respiratory', 'cardiac'});
verbose.fig_handles(end+1) = fh;
end
r = rInterp;

else
t = tRespiration;
cInterp = interp1(tCardiac, c, t);
%cacq_codesInterp = interp1(tCardiac, cacq_codes, t, 'nearest');
%acq_codes = racq_codes + cacq_codesInterp;

if DEBUG
fh = plot_interpolation(tCardiac, c, t, cInterp, ...
{'cardiac', 'respiratory'});
verbose.fig_handles(end+1) = fh;
end
c = cInterp;

end

%% Recompute acq_codes as for Siemens (volume on/volume off)
acq_codes = [];

if ~isempty(iAcqStart) % otherwise, nothing to read ...
% iAcqStart is a columns of 0 and 1, 1 for the trigger event of a new
% volume start

% sometimes, trigger is on for several samples; ignore these extended
% phases of "on-triggers" as duplicate values, if trigger distance is
% very different
%
% fraction of mean trigger distance; if trigger time difference below that, will be removed
outlierThreshold = 0.2;

idxAcqStart = find(iAcqStart);
dAcqStart = diff(idxAcqStart);

% + 1 because of diff
iAcqOutlier = 1 + find(dAcqStart < outlierThreshold*mean(dAcqStart));
iAcqStart(idxAcqStart(iAcqOutlier)) = 0;

acq_codes = zeros(nSamplesCardiac,1);
acq_codes(iAcqStart) = 8; % to match Philips etc. format

nAcqs = numel(find(iAcqStart));

if nAcqs >= 1
% report time of acquisition, as defined in SPM
meanTR = mean(diff(t(iAcqStart)));
stdTR = std(diff(t(iAcqStart)));
verbose = tapas_physio_log(...
sprintf('TR = %.3f +/- %.3f s (Estimated mean +/- std time of repetition for one volume; nTriggers = %d)', ...
meanTR, stdTR, nAcqs), verbose, 0);
end
end

%% Plot, if wanted

if DEBUG
stringTitle = 'Read-In: Raw BIDS physlog data (TSV file)';
verbose.fig_handles(end+1) = ...
tapas_physio_plot_raw_physdata_siemens_hcp(t, c, r, acq_codes, ...
stringTitle);
end

end

%% Function for unified physio file (cardiac + resp)
function [c,r,t,cpulse,acq_codes] = tapas_physio_read_physlogfiles_bids_unified(fileName,log_files,hasExplicitJsonFile,DEBUG,verbose)
%% check for .gz files and unzip to temporary file
[~, ~, ext] = fileparts(fileName);

isZipped = strcmpi(ext, '.gz');

Expand All @@ -105,13 +367,7 @@
fileJson = log_files.scan_timing;
end

% 'isfile' is clearer than 'exist', because checks for correct absolute/relative path, not existence anywhere on Matlab path
% but only available in Matlab version >= R2017b
if exist('isfile', 'builtin')
hasJsonFile = isfile(fileJson);
else
hasJsonFile = exist(fileJson, 'file') > 0;
end
hasJsonFile = isfile(fileJson);

if hasJsonFile
val = jsondecode(fileread(fileJson));
Expand Down Expand Up @@ -189,7 +445,7 @@
%
% fraction of mean trigger distance; if trigger time difference below that, will be removed
outlierThreshold = 0.2;

idxAcqStart = find(iAcqStart);
dAcqStart = diff(idxAcqStart);

Expand All @@ -212,7 +468,6 @@
end
end


%% Plot, if wanted

if DEBUG
Expand All @@ -222,6 +477,7 @@
stringTitle);
end


%% occasionally, cardiac time course is instead containing 0/1 cardiac triggers,
% and not raw trace; check this and populate cpulse accordingly
if all(ismember(unique(c), [1 0]))
Expand All @@ -238,3 +494,22 @@
tapas_physio_log(sprintf('%s: %s', messageId, message), verbose, 1)
end
end

end

%% Local function to plot interpolation result
function fh = plot_interpolation(tOrig, yOrig, tInterp, yInterp, ...
stringOrigInterp)
fh = tapas_physio_get_default_fig_params();
stringTitle = sprintf('Read-In: Interpolation of %s signal', stringOrigInterp{1});
set(fh, 'Name', stringTitle);
plot(tOrig, yOrig, 'go--'); hold all;
plot(tInterp, yInterp,'r.');
legend({
sprintf('after interpolation to %s timing', ...
stringOrigInterp{1}), ...
sprintf('original %s time series', stringOrigInterp{2}) });
title(stringTitle);
xlabel('time (seconds');
end

0 comments on commit 2f07036

Please sign in to comment.