readCMRRPhysio.m

function varargout = readCMRRPhysio(varargin)
% -------------------------------------------------------------------------
% readCMRRPhysio.m
% -------------------------------------------------------------------------
% Read physiological log files from CMRR MB sequences (>=R013, >=VD13A)
%   E. Auerbach, CMRR, 2015-2023
%
% Usage #1 (individual .log files):
%    physio = readCMRRPhysio(base_filename, [show_plot]);
% Usage #2 (single encoded DICOM file):
%    physio = readCMRRPhysio(DICOM_filename, [show_plot, [output_path]]);
%
% This function expects to find either some combination of individual log
% files (*_ECG.log, *_RESP.log, *_PULS.log, *_EXT.log, *_Info.log)
% generated by the CMRR C2P sequences >=R013, or a single encoded "_PHYSIO"
% DICOM file generated by the CMRR C2P sequences >=R015.
%
% Inputs:
%    base_filename  = 'Physio_DATE_TIME_UUID'
%    DICOM_filename = 'XXX.dcm'
%    show_plot = 1 to graphically display traces after import (optional)
%    output_path = '/path/to/output/' (optional; path to write .log files)
%
% Returns:
%    Physio traces will be returned for ECG1, ECG2, ECG3, ECG4, RESP, PULS,
%    EXT1, and EXT2 signals. Only active traces (with nonzero values) will
%    be returned.
%
%    If input is DICOM format and the optional output path is specified,
%    the extracted log data will be written to disk (this functionality
%    replaces extractCMRRPhysio.m).
%
%    The unit of time is clock ticks (2.5 ms per tick).
%        physio.UUID: unique identifier string for this measurement
%        physio.SliceMap: [2 x Volumes x Slices] array
%            (1,:,:) = start time stamp of each volume/slice
%            (2,:,:) = finish time stamp of each volume/slice
%        physio.ACQ: [total scan time x 1] array
%            value = 1 if acquisition is active at this time; 0 if not
%        physio.ECG1: [total scan time x 1] array
%        physio.ECG2: [total scan time x 1] array
%        physio.ECG3: [total scan time x 1] array
%        physio.ECG4: [total scan time x 1] array
%            value = ECG signal on this channel
%        physio.RESP: [total scan time x 1] array
%            value = RESP signal on this channel
%        physio.PULS: [total scan time x 1] array
%            value = PULS signal on this channel
%        physio.EXT: [total scan time x 1] array
%            value = 1 if EXT signal detected; 0 if not
%        physio.EXT2: [total scan time x 1] array
%            value = 1 if EXT2 signal detected; 0 if not

VersionString = '2023.11.09';

% this is the file format this function expects; must match log file version
ExpectedVersion = 'EJA_1';

% say hello
fprintf('\nreadCMRRPhysio version %s: E. Auerbach, CMRR\n\n', VersionString);

% check input arguments
show_plot = 0;
outpath = [];
if (nargin < 1) || (nargin > 3)
    usage;
    error('Invalid number of inputs.');
end
fn = varargin{1};
if (nargin == 2)
    test_var = varargin{2};
    if (isnumeric(test_var) || islogical(test_var))
        show_plot = test_var;
    else
        outpath = test_var;
    end
elseif (nargin == 3)
    show_plot = varargin{2};
    outpath = varargin{3};
end
if (~isempty(outpath) && (exist(outpath, 'dir') ~= 7))
    error('Could not locate requested output path: %s', outpath);
end

if (nargout > 1)
    usage;
    error('Invalid number of outputs.');
elseif (~nargout && isempty(outpath) && ~show_plot)
    usage;
    error('Nothing to do!');
end

% first, check if the base is pointing to a DICOM we should extract
if (2 == exist(fn,'file'))
    fprintf('Found %s\n', fn);
    if (~isdicom(fn))
        usage;
        error('Invalid syntax for text log file import!');
    end
    fprintf('Attempting to read CMRR Physio DICOM format file...\n');
    warning('off','images:dicominfo:attrWithSameName');
    dcmInfo = dicominfo(fn);
    if (isempty(dcmInfo)), error('%s could not be read as a valid DICOM format file!', fn); end
    
    % all versions for VE-line store data in Private_7fe1_10xx_Creator
    if (isfield(dcmInfo,'ImageType') && strcmp(dcmInfo.ImageType,'ORIGINAL\PRIMARY\RAWDATA\PHYSIO') ...
        && isfield(dcmInfo,'Private_7fe1_10xx_Creator') && strcmp(deblank(char(dcmInfo.Private_7fe1_10xx_Creator)),'SIEMENS CSA NON-IMAGE'))
        pdata = dcmInfo.Private_7fe1_1010;

    % case 2: up until R017pre6, XA-line ImageType field is wrong, stores in SpectroscopyData
    elseif (isfield(dcmInfo,'ImageType') && strcmp(dcmInfo.ImageType,'ORIGINAL\PRIMARY\RAWDATA\NONE') ... % XA30 bug
        && isfield(dcmInfo,'SpectroscopyData'))
        pdata = typecast(dcmInfo.SpectroscopyData, 'uint8');

    % case 3: R017pre7 and newer, XA-line ImageType field is wrong, but stored in Private_7fe1_10xx_Creator again
    elseif (isfield(dcmInfo,'ImageType') ...
        && (strcmp(dcmInfo.ImageType,'ORIGINAL\PRIMARY\RAWDATA\NONE') || strcmp(dcmInfo.ImageType,'ORIGINAL\PRIMARY\OTHER\NONE')) ...
        && isfield(dcmInfo,'Private_7fe1_10xx_Creator') && strcmp(deblank(char(dcmInfo.Private_7fe1_10xx_Creator)),'SIEMENS MR IMA'))
        pdata = dcmInfo.Private_7fe1_1010;

    % unknown case
    else
        error('could not find physio data in %s (DICOM format)!', fn);
    end
    
    % read in the data
    np = size(pdata,1);
    rows = dcmInfo.AcquisitionNumber;
    columns = np/rows;
    numFiles = columns/1024;
    if (rem(np,rows) || rem(columns,1024)), error('Invalid image size (%dx%d)!', columns, rows); end
    dcmData = reshape(pdata,[],numFiles)';
    % encoded DICOM format: columns = 1024*numFiles
    %                       first row: uint32 datalen, uint32 filenamelen, char[filenamelen] filename
    %                       remaining rows: char[datalen] data
    foundECG  = 0;
    foundRESP = 0;
    foundPULS = 0;
    foundEXT  = 0;
    [~,~,endian] = computer;
    needswap = ~strcmp(endian,'L');
    for idx=1:numFiles
        datalen = typecast(dcmData(idx,1:4),'uint32');
        if needswap, datalen = swapbytes(datalen); end
        filenamelen = typecast(dcmData(idx,5:8),'uint32');
        if needswap, filenamelen = swapbytes(filenamelen); end
        filename = char(dcmData(idx,9:9+filenamelen-1));
        logData = dcmData(idx,1025:1025+datalen-1);
        fprintf('  Decoded: %s\n', filename);
        if (strcmp(filename(end-9+1:end),'_Info.log'))
            fnINFO = logData;
        elseif (strcmp(filename(end-8+1:end),'_ECG.log'))
            fnECG  = logData;
            foundECG = 1;
        elseif (strcmp(filename(end-9+1:end),'_RESP.log'))
            fnRESP = logData;
            foundRESP = 1;
        elseif (strcmp(filename(end-9+1:end),'_PULS.log'))
            fnPULS = logData;
            foundPULS = 1;
        elseif (strcmp(filename(end-8+1:end),'_EXT.log'))
            fnEXT = logData;
            foundEXT = 1;
        end
        if ~isempty(outpath)
            outfn = fullfile(outpath, filename);
            fprintf('  Writing: %s\n', outfn);
            fp = fopen(outfn,'w');
            fwrite(fp, char(logData));
            fclose(fp);
        end
    end
    fprintf('\n');
    
% if we don't have an encoded DICOM, check what text log files we have
else
    if ~isempty(outpath)
        usage;
        error('Invalid syntax for text log file import!');
    end

    fnINFO = [fn '_Info.log'];
    fnECG  = [fn '_ECG.log'];
    fnRESP = [fn '_RESP.log'];
    fnPULS = [fn '_PULS.log'];
    fnEXT  = [fn '_EXT.log'];
    if (2 ~= exist(fnINFO, 'file')), error('%s not found!', fnINFO); end
    foundECG  = (2 == exist(fnECG , 'file'));
    foundRESP = (2 == exist(fnRESP, 'file'));
    foundPULS = (2 == exist(fnPULS, 'file'));
    foundEXT  = (2 == exist(fnEXT , 'file'));
end

if (~foundECG && ~foundRESP && ~foundPULS && ~foundEXT)
    warning('No data files (ECG/RESP/PULS/EXT) found!');
    fprintf('\n');
end

% if we wrote the log files and are not plotting or returning the parsed data, we are done
if (~isempty(outpath) && ~show_plot && ~nargout), return; end

% read in the data
[SliceMap, UUID1, NumSlices, NumVolumes, FirstTime, LastTime, NumEchoes] = readParseFile(fnINFO, 'ACQUISITION_INFO', ExpectedVersion, 0, 0);
if (LastTime <= FirstTime), error('Last timestamp is not greater than first timestamp, aborting...'); end
ActualSamples = LastTime - FirstTime + 1;
ExpectedSamples = ActualSamples + 8; % some padding at the end for worst case EXT sample at last timestamp

if (foundECG)
    [ECG, UUID2] = readParseFile(fnECG, 'ECG', ExpectedVersion, FirstTime, ExpectedSamples);
    if (~strcmp(UUID1, UUID2)), error('UUID mismatch between Info and ECG files!'); end
end

if (foundRESP)
    [RESP, UUID3] = readParseFile(fnRESP, 'RESP', ExpectedVersion, FirstTime, ExpectedSamples);
    if (~strcmp(UUID1, UUID3)), error('UUID mismatch between Info and RESP files!'); end
end

if (foundPULS)
    [PULS, UUID4] = readParseFile(fnPULS, 'PULS', ExpectedVersion, FirstTime, ExpectedSamples);
    if (~strcmp(UUID1, UUID4)), error('UUID mismatch between Info and PULS files!'); end
end

if (foundEXT)
    [EXT, UUID5] = readParseFile(fnEXT, 'EXT', ExpectedVersion, FirstTime, ExpectedSamples);
    if (~strcmp(UUID1, UUID5)), error('UUID mismatch between Info and EXT files!'); end
end

fprintf('Formatting data...\n');
ACQ = zeros(ExpectedSamples,1,'uint16');
for v=1:NumVolumes
    for s=1:NumSlices
        for e=1:NumEchoes
            ACQ(SliceMap(1,v,s,e)+1:SliceMap(2,v,s,e)+1,1) = 1;
        end
    end
end
    
fprintf('\n');
fprintf('Slices in scan:      %d\n', NumSlices);
fprintf('Volumes in scan:     %d\n', NumVolumes);
if (NumEchoes > 1)
    fprintf('Echoes per slc/vol:  %d\n', NumEchoes);
end
fprintf('First timestamp:     %d\n', FirstTime);
fprintf('Last timestamp:      %d\n', LastTime);
fprintf('Total scan duration: %d ticks\n', ActualSamples);
fprintf('Total scan duration: %.4f s\n', double(ActualSamples)*2.5/1000);
fprintf('\n');

% only return active (nonzero) traces
physio.UUID = UUID1;
physio.SliceMap = SliceMap;
physio.ACQ = ACQ;
if ((1 == exist('ECG' , 'var')) && ~isempty(ECG)  && nnz(ECG(:,1))), physio.ECG1 = ECG(:,1); end
if ((1 == exist('ECG' , 'var')) && ~isempty(ECG)  && nnz(ECG(:,2))), physio.ECG2 = ECG(:,2); end
if ((1 == exist('ECG' , 'var')) && ~isempty(ECG)  && nnz(ECG(:,3))), physio.ECG3 = ECG(:,3); end
if ((1 == exist('ECG' , 'var')) && ~isempty(ECG)  && nnz(ECG(:,4))), physio.ECG4 = ECG(:,4); end
if ((1 == exist('RESP', 'var')) && ~isempty(RESP) && nnz(RESP))    , physio.RESP = RESP;     end
if ((1 == exist('PULS', 'var')) && ~isempty(PULS) && nnz(PULS))    , physio.PULS = PULS;     end
if ((1 == exist('EXT' , 'var')) && ~isempty(EXT)  && nnz(EXT(:,1))), physio.EXT  = EXT(:,1); end
if ((1 == exist('EXT' , 'var')) && ~isempty(EXT)  && nnz(EXT(:,2))), physio.EXT2 = EXT(:,2); end

% plot data in a rudimentary way if requested
% if too large, only plot the middle 1k ticks or so
if (show_plot)
    display_max = 1000;
    start_tick = 1;
    end_tick = ActualSamples;
    if (ActualSamples > display_max)
        start_tick = floor(ActualSamples/2) - floor(display_max/2) + 1;
        end_tick = start_tick + display_max - 1;
    end
    figure;
    hold on;
    miny = 50000; maxy = -50000; % actual range is 0..4095
    if (isfield(physio,'ECG1')), [miny, maxy] = plot_trace(physio.ECG1(start_tick:end_tick), miny, maxy, 'y', false); end
    if (isfield(physio,'ECG2')), [miny, maxy] = plot_trace(physio.ECG2(start_tick:end_tick), miny, maxy, 'y', false); end
    if (isfield(physio,'ECG3')), [miny, maxy] = plot_trace(physio.ECG3(start_tick:end_tick), miny, maxy, 'y', false); end
    if (isfield(physio,'ECG4')), [miny, maxy] = plot_trace(physio.ECG4(start_tick:end_tick), miny, maxy, 'y', false); end
    if (isfield(physio,'RESP')), [miny, maxy] = plot_trace(physio.RESP(start_tick:end_tick), miny, maxy, 'm', false); end
    if (isfield(physio,'PULS')), [miny, maxy] = plot_trace(physio.PULS(start_tick:end_tick), miny, maxy, 'r', false); end
    if (isfield(physio,'EXT' )), [miny, maxy] = plot_trace(physio.EXT (start_tick:end_tick), miny, maxy, 'c', true);  end
    if (isfield(physio,'EXT2')), [miny, maxy] = plot_trace(physio.EXT2(start_tick:end_tick), miny, maxy, 'g', true);  end
    [miny, maxy] = plot_trace(physio.ACQ(start_tick:end_tick), miny, maxy, 'k', true);
    axis([1 double(min(display_max, ActualSamples)) miny-maxy*0.05 maxy+maxy*0.05]);
end

if (nargout)
    varargout{1} = physio;
end

%--------------------------------------------------------------------------

function [arr, varargout] = readParseFile(fn, LogDataType, ExpectedVersion, FirstTime, ExpectedSamples)
% read and parse log file

if (isa(fn,'uint8'))
    % if fn is uint8, we read it directly from DICOM
    fprintf('Parsing %s data...\n', LogDataType);
    inData = char(fn);
else
    % otherwise, fn is a filename
    fprintf('Reading %s file...\n', LogDataType);
    fp = fopen(fn);
    inData = fread(fp, Inf, '*char');
    fclose(fp);
end

% echoes parameter was not added until R015a, so prefill a default value
% for compatibility with older data
NumEchoes = uint16(1);
varargout{6} = NumEchoes;

% mgetl returns a cell array where each cell is a line of text
[lines, numlines] = mgetl(inData);

arr = [];

% loop through and clean up
for curline=1:numlines
    line = lines{curline};
    
    % strip whitespace
    if (~isempty(line)), line = strtrim(line); end

    % strip any comments
    if (strfind(line, '#') > 1), line = strtrim(line(1:ctest-1)); end
    
    lines{curline} = line;
end

% now loop through find all the header entries, which since XA61 can also be at the end of the file
for curline=1:numlines
    line = lines{curline};

    if (~isempty(line))
        if (contains(line, '='))
            % this is an assigned value; parse it
            varcell = textscan(line, '%s=%s');
            varname = strtrim(varcell{1});
            value   = strtrim(varcell{2});
            
            if (strcmp(varname, 'UUID')), varargout{1} = value; end
            %if (strcmp(varname, 'ScanDate')), ScanDate = value; end
            if (strcmp(varname, 'LogVersion'))
                if (~strcmp(value, ExpectedVersion))
                    error('File format [%s] not supported by this function (expected [%s]).', value, ExpectedVersion);
                end
            end
            if (strcmp(varname, 'LogDataType'))
                if (~strcmp(value, LogDataType))
                    error('Expected [%s] data, found [%s]? Check filenames?', LogDataType, value);
                end
            end
            if (strcmp(varname, 'SampleTime'))
                if (strcmp(LogDataType, 'ACQUISITION_INFO'))
                    error('Invalid [%s] parameter found.',varname);
                end
                SampleTime = uint16(str2double(value));
            end
            if (strcmp(varname, 'NumSlices'))
                if (~strcmp(LogDataType, 'ACQUISITION_INFO'))
                    error('Invalid [%s] parameter found.',varname);
                end
                NumSlices = uint16(str2double(value));
                varargout{2} = NumSlices;
            end
            if (strcmp(varname, 'NumVolumes'))
                if (~strcmp(LogDataType, 'ACQUISITION_INFO'))
                    error('Invalid [%s] parameter found.',varname);
                end
                NumVolumes = uint16(str2double(value));
                varargout{3} = NumVolumes;
            end
            if (strcmp(varname, 'FirstTime'))
                if (~strcmp(LogDataType, 'ACQUISITION_INFO'))
                    error('Invalid [%s] parameter found.',varname);
                end
                FirstTime = uint32(str2double(value));
                varargout{4} = FirstTime;
            end
            if (strcmp(varname, 'LastTime'))
                if (~strcmp(LogDataType, 'ACQUISITION_INFO'))
                    error('Invalid [%s] parameter found.',varname);
                end
                varargout{5} = uint32(str2double(value));
            end
            if (strcmp(varname, 'NumEchoes'))
                if (~strcmp(LogDataType, 'ACQUISITION_INFO'))
                    error('Invalid [%s] parameter found.',varname);
                end
                NumEchoes = uint16(str2double(value));
                varargout{6} = NumEchoes;
            end
        end
    end
end

% finally find and parse the data
for curline=1:numlines
    line = lines{curline};

    if (~isempty(line))
        if (~contains(line, '='))
            % this must be data; currently it is 3-5 columns so we can
            % parse it easily with textscan
            datacells = textscan(line, '%s %s %s %s %s');
            
            if (~isstrprop(datacells{1}{1}(1), 'digit'))
                % if the first column isn't numeric, it is probably the header
            else
                % store data in output array based on the file type
                if (strcmp(LogDataType, 'ACQUISITION_INFO'))
                    if ( (1 ~= exist('NumVolumes', 'var')) || (NumVolumes < 1) || ...
                            (1 ~= exist('NumSlices' , 'var')) || (NumSlices  < 1) || ...
                            (1 ~= exist('NumEchoes' , 'var')) || (NumEchoes  < 1) )
                        error('Failed reading ACQINFO header!');
                    end
                    if (NumVolumes == 1)
                        % this is probably R016a or earlier diffusion data, where NumVolumes is 1 (incorrect)
                        NumVolumes = (numlines-11)/(NumSlices*NumEchoes);
                        warning('Found NumVolumes=1; correcting to %d for R016a and earlier diffusion data!', NumVolumes);
                    end
                    if (isempty(arr)), arr = zeros(2,NumVolumes,NumSlices,NumEchoes,'uint32'); end
                    curvol    = uint16(str2double(datacells{1}{1})) + 1;
                    curslc    = uint16(str2double(datacells{2}{1})) + 1;
                    curstart  = uint32(str2double(datacells{3}{1}));
                    curfinish = uint32(str2double(datacells{4}{1}));
                    if (size(datacells{5},1))
                        cureco = uint16(str2double(datacells{5}{1})) + 1;
                        if (arr(:,curvol,curslc,cureco)), error('Received duplicate timing data for vol%d slc%d eco%d!', curvol, curslc, cureco); end
                    else
                        cureco = uint16(0) + 1;
                        if (arr(:,curvol,curslc,cureco)), warning('Received duplicate timing data for vol%d slc%d (ignore for pre-R015a multi-echo data)!', curvol, curslc); end
                    end
                    arr(:,curvol,curslc,cureco) = [curstart curfinish]; %#ok<AGROW>
                else
                    curstart   = uint32(str2double(datacells{1}{1})) - FirstTime + 1;
                    curchannel = datacells{2}{1};
                    curvalue   = uint16(str2double(datacells{3}{1}));
                    %curtrigger = datacells{4}{1};
                    
                    if (strcmp(LogDataType, 'ECG'))
                        if (isempty(arr)), arr = zeros(ExpectedSamples,4,'uint16'); end
                        if (strcmp(curchannel, 'ECG1'))
                            chaidx = 1;
                        elseif (strcmp(curchannel, 'ECG2'))
                            chaidx = 2;
                        elseif (strcmp(curchannel, 'ECG3'))
                            chaidx = 3;
                        elseif (strcmp(curchannel, 'ECG4'))
                            chaidx = 4;
                        else
                            error('Invalid ECG channel ID [%s]', curchannel);
                        end
                    elseif (strcmp(LogDataType, 'EXT'))
                        if (isempty(arr)), arr = zeros(ExpectedSamples,2,'uint16'); end
                        if (strcmp(curchannel, 'EXT'))
                            chaidx = 1;
                        elseif (strcmp(curchannel, 'EXT2'))
                            chaidx = 2;
                        else
                            error('Invalid EXT channel ID [%s]', curchannel);
                        end
                    else
                        if (isempty(arr)), arr = zeros(ExpectedSamples,1,'uint16'); end
                        chaidx = 1;
                    end
                    
                    arr(curstart:curstart+uint32(SampleTime-1),chaidx) = curvalue*ones(SampleTime,1,'uint16'); %#ok<AGROW>
                end
            end
        end
    end
end

if (strcmp(LogDataType, 'ACQUISITION_INFO'))
    arr = arr - FirstTime;
end

%--------------------------------------------------------------------------

function [lines, numlines] = mgetl(arr)
% mgetl: parse an entire text file into cell array of strings where each
% cell is one line. recognizes dos and unix file formats.

arr = char(arr);
if (size(arr,1) > size(arr,2)), arr = arr'; end
arr_len = length(arr);
lf = strfind(arr, newline);
numlines = length(lf);
if (lf(numlines) < arr_len)
    numlines = numlines + 1;
    lf(numlines) = arr_len;
end

lines = cell(numlines, 1);

stpos = 1;
for x=1:numlines
    if (x > 1), stpos = lf(x-1)+1; end
    endpos = lf(x) - 1;
    if (lf(x) > 1)
        if (arr(endpos) == char(13))
            endpos = endpos - 1; % strip cr and lf if both present
        end
    end
    if (endpos >= stpos), lines{x} = arr(stpos:endpos); end
end

%--------------------------------------------------------------------------

function [newminy, newmaxy] = plot_trace(trace, oldminy, oldmaxy, color, scale)
% plot trace and keep track of minimum and maximum values

miny = double(min(trace));
maxy = double(max(trace));
newminy = min(oldminy, miny);
newmaxy = max(oldmaxy, maxy);
if (scale && ((miny ~= oldminy) || (maxy ~= oldmaxy)))
    trace = double(trace) .* ((newmaxy-newminy)/(maxy-miny));
    trace = trace - min(trace) + newminy;
end
plot(trace,'Color',color);

%--------------------------------------------------------------------------

function usage
% display command syntax

fprintf(['Usage #1 (individual .log files):' ...
         '\n   physio = readCMRRPhysio(base_filename, [show_plot])' ...
         '\nUsage #2 (single encoded DICOM file):' ...
         '\n   physio = readCMRRPhysio(DICOM_filename, [show_plot, [output_path]])\n\n']);