diff --git a/Contributor-License-Agreement.md b/Contributor-License-Agreement.md index 80d82149..618c2be1 100644 --- a/Contributor-License-Agreement.md +++ b/Contributor-License-Agreement.md @@ -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") diff --git a/PhysIO/code/readin/tapas_physio_read_physlogfiles_bids.m b/PhysIO/code/readin/tapas_physio_read_physlogfiles_bids.m index d40cdf37..45840fd0 100644 --- a/PhysIO/code/readin/tapas_physio_read_physlogfiles_bids.m +++ b/PhysIO/code/readin/tapas_physio_read_physlogfiles_bids.m @@ -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 @@ -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 @@ -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'); @@ -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)); @@ -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); @@ -212,7 +468,6 @@ end end - %% Plot, if wanted if DEBUG @@ -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])) @@ -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 +