eidors_readdata

PURPOSE ^

EIDORS readdata - read data files from various EIT equipment

SYNOPSIS ^

function [vv, auxdata, stim]= eidors_readdata( fname, format, frame_range, extra )

DESCRIPTION ^

 EIDORS readdata - read data files from various EIT equipment
    manufacturers

 [vv, auxdata, stim ]= eidors_readdata( fname, format, frame_range, extra )
    frame_range is the list of frames to load from the file (ie. 1:100).
    if the frames are beyond the number in the file, no data is returned
    to get all frames, use frame_range= []

 Currently the list of supported file formats is:
    - "EIT" Files, may be one of 
          - Draeger Pulmovista (2008+)
          - GoeIIMF/Carefusion (2010+)
          - Swisstom BB2/Pioneer (2010+)
       The function will attempt to autodetect the format
       
    - MCEIT (GoeIIMF / Viasys) "get" file format 
        format = "GET" or "MCEIT"
        Note that the data is "untwisted" to correspond to a "no_rotate_meas" stim pattern
    - MCEIT (GoeIIMF) "get" file format
        format = "GET-RAW"
        Data in original order, corresponds to "rotate_meas" stim pattern

    - Draeger (Pulmovista) "EIT" file format (2008+)
        format = "DRAEGER-EIT"

    - Draeger "get" file format (- 2007 - format for Draeger equipment)
        format = "DRAEGER-GET"

    - Sentec (Swisstom) EIT equipment "EIT" (Pioneer set and BB2):
           'LQ1' (2010 - 2011)
           'LQ2' (2013 - 2014)
           'LQ4' (2015+)
           'LQ5' (2023+)

    - Dixtal file format, from Dixtal inc, Brazil
        format = 'DIXTAL_encode' extract encoder from provided Dll
           - Output: [encodepage] = eidors_readdata( path,'DIXTAL_encode');
              where path= '/path/to/Criptografa_New.dll' (provided with system)
        format = 'DIXTAL'
           - output: [vv] = eidors_readdata( fname, 'DIXTAL', [], encodepage );
    - New Carefusion "EIT" file format
        format = "GOEIIMF-EIT" or "carefusion"

    - HDF5-based format
        format = 'h5' or 'hdf5' documented at
       Possner, Bulst, Adler "HDF5-based data format for EIT data", Conf. EIT2023
        the frame_range parameter specifies which dataset to load
        e.g. eidors_readdata('montreal_data_1995_breathold.h5','h5','datasetEIT')
       If format 'h5-all' or 'hdf5-all' or dataset = '{:all:}', then all datasets are
        loaded into a cell array

    - Medas Medical 'VTM' format
       For the Ventom device

    - Sciospec text 'EIT' format
        format = 'Sciospec-EIT'

    - Sheffield MK I "RAW" file format
        format = "RAW" or "sheffield"

    - ITS (International Tomography Systems)
        format = "ITS" or "p2k"

    - IIRC (Impedance Imaging Research Center, Korea)
        format = "txt" or "IIRC"

    - University of Cape Town formats
        format = "UCT_SEQ"  UCT sequence file
           - Output: [ stimulation, meas_select]= eidors_readdata(fname, 'UCT_SEQ')
        format = "UCT_CAL"  UCT calibration file
           - Output: [vv, no_cur_caldata_raw ]= eidors_readdata( fname, 'UCT_CAL' )
                 where no_cur_caldata_raw is data captured with no current
        format = "UCT_DATA"  UCT data frame file
           - Output: [vv]= eidors_readdata( fname, 'UCT_DATA' )

 Usage
 [vv, auxdata, stim ]= eidors_readdata( fname, format )
     vv      = measurements - data frames in each column
     auxdata = auxillary data - if provided by system 
     stim    = stimulation structure, to be used with
                fwd_model.stimulation. 
     fname = file name
     stim.framerate = acquisition rate (frames/sec) if available

  if format is unspecified, an attempt to autodetect is made
    auxdata.format is the detected format

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [vv, auxdata, stim]= eidors_readdata( fname, format, frame_range, extra )
0002 % EIDORS readdata - read data files from various EIT equipment
0003 %    manufacturers
0004 %
0005 % [vv, auxdata, stim ]= eidors_readdata( fname, format, frame_range, extra )
0006 %    frame_range is the list of frames to load from the file (ie. 1:100).
0007 %    if the frames are beyond the number in the file, no data is returned
0008 %    to get all frames, use frame_range= []
0009 %
0010 % Currently the list of supported file formats is:
0011 %    - "EIT" Files, may be one of
0012 %          - Draeger Pulmovista (2008+)
0013 %          - GoeIIMF/Carefusion (2010+)
0014 %          - Swisstom BB2/Pioneer (2010+)
0015 %       The function will attempt to autodetect the format
0016 %
0017 %    - MCEIT (GoeIIMF / Viasys) "get" file format
0018 %        format = "GET" or "MCEIT"
0019 %        Note that the data is "untwisted" to correspond to a "no_rotate_meas" stim pattern
0020 %    - MCEIT (GoeIIMF) "get" file format
0021 %        format = "GET-RAW"
0022 %        Data in original order, corresponds to "rotate_meas" stim pattern
0023 %
0024 %    - Draeger (Pulmovista) "EIT" file format (2008+)
0025 %        format = "DRAEGER-EIT"
0026 %
0027 %    - Draeger "get" file format (- 2007 - format for Draeger equipment)
0028 %        format = "DRAEGER-GET"
0029 %
0030 %    - Sentec (Swisstom) EIT equipment "EIT" (Pioneer set and BB2):
0031 %           'LQ1' (2010 - 2011)
0032 %           'LQ2' (2013 - 2014)
0033 %           'LQ4' (2015+)
0034 %           'LQ5' (2023+)
0035 %
0036 %    - Dixtal file format, from Dixtal inc, Brazil
0037 %        format = 'DIXTAL_encode' extract encoder from provided Dll
0038 %           - Output: [encodepage] = eidors_readdata( path,'DIXTAL_encode');
0039 %              where path= '/path/to/Criptografa_New.dll' (provided with system)
0040 %        format = 'DIXTAL'
0041 %           - output: [vv] = eidors_readdata( fname, 'DIXTAL', [], encodepage );
0042 %    - New Carefusion "EIT" file format
0043 %        format = "GOEIIMF-EIT" or "carefusion"
0044 %
0045 %    - HDF5-based format
0046 %        format = 'h5' or 'hdf5' documented at
0047 %       Possner, Bulst, Adler "HDF5-based data format for EIT data", Conf. EIT2023
0048 %        the frame_range parameter specifies which dataset to load
0049 %        e.g. eidors_readdata('montreal_data_1995_breathold.h5','h5','datasetEIT')
0050 %       If format 'h5-all' or 'hdf5-all' or dataset = '{:all:}', then all datasets are
0051 %        loaded into a cell array
0052 %
0053 %    - Medas Medical 'VTM' format
0054 %       For the Ventom device
0055 %
0056 %    - Sciospec text 'EIT' format
0057 %        format = 'Sciospec-EIT'
0058 %
0059 %    - Sheffield MK I "RAW" file format
0060 %        format = "RAW" or "sheffield"
0061 %
0062 %    - ITS (International Tomography Systems)
0063 %        format = "ITS" or "p2k"
0064 %
0065 %    - IIRC (Impedance Imaging Research Center, Korea)
0066 %        format = "txt" or "IIRC"
0067 %
0068 %    - University of Cape Town formats
0069 %        format = "UCT_SEQ"  UCT sequence file
0070 %           - Output: [ stimulation, meas_select]= eidors_readdata(fname, 'UCT_SEQ')
0071 %        format = "UCT_CAL"  UCT calibration file
0072 %           - Output: [vv, no_cur_caldata_raw ]= eidors_readdata( fname, 'UCT_CAL' )
0073 %                 where no_cur_caldata_raw is data captured with no current
0074 %        format = "UCT_DATA"  UCT data frame file
0075 %           - Output: [vv]= eidors_readdata( fname, 'UCT_DATA' )
0076 %
0077 % Usage
0078 % [vv, auxdata, stim ]= eidors_readdata( fname, format )
0079 %     vv      = measurements - data frames in each column
0080 %     auxdata = auxillary data - if provided by system
0081 %     stim    = stimulation structure, to be used with
0082 %                fwd_model.stimulation.
0083 %     fname = file name
0084 %     stim.framerate = acquisition rate (frames/sec) if available
0085 %
0086 %  if format is unspecified, an attempt to autodetect is made
0087 %    auxdata.format is the detected format
0088 
0089 % (C) 2005-09 Andy Adler. License: GPL version 2 or version 3
0090 % $Id: eidors_readdata.m 7131 2024-12-29 16:23:02Z aadler $
0091 
0092 % TODO:
0093 %   - output an eidors data object
0094 %   - test whether file format matches specified stimulation pattern
0095 %   - todo MCEIT provides curr and volt on driven electrodes.
0096 %       how can this be provided to system?
0097 
0098 if nargin>=1 && strcmp(fname,'UNIT_TEST'); do_unit_test; return; end
0099 
0100 check_file_exists(fname);
0101 
0102 if nargin < 2
0103 % unspecified file format, autodetect
0104    dotpos = find(fname == '.');
0105    if isempty( dotpos ) 
0106       error('file format unspecified, can`t autodetect');
0107    else
0108       dotpos= dotpos(end);
0109       format= fname( dotpos+1:end );
0110    end
0111 end
0112 
0113 
0114 auxdata = struct(); % default, can be overriden if the format has it
0115 fmt = pre_proc_spec_fmt( format, fname );
0116 switch fmt
0117    case {'h5','hdf5'};
0118       fmt = 'hdf5';
0119       % Possner, Bulst, Adler "HDF5-based data format for EIT data", Conf. EIT2023
0120       if nargin<3; frame_range=[]; end
0121       [vv, auxdata, stim]= h5_readdata( fname, frame_range );
0122 
0123    case {'h5-all','hdf5-all'};
0124       fmt = 'hdf5';
0125       frame_range = '{:all:}';
0126       [vv, auxdata, stim]= h5_readdata( fname, frame_range );
0127 
0128    case 'mceit';
0129       [vv,curr,volt,auxdata_out,auxtime,rawdata] = mceit_readdata( fname );
0130       auxdata.auxdata = auxdata_out;
0131       auxdata.auxtime = auxtime;
0132       auxdata.curr    = curr;
0133       auxdata.volt    = volt;
0134       auxdata.frame_rate = NaN; % This info is only in the *.prl file
0135 
0136 
0137       if strcmp(lower(format),'get-raw')
0138           vv= rawdata(1:208,:);
0139           stim = mk_stim_patterns(16,1,[0,1],[0,1], {'no_meas_current','rotate_meas'}, 1);
0140       else
0141           stim = basic_stim(16);
0142       end
0143    case 'draeger-get'
0144       vv = draeger_get_readdata( fname );
0145 
0146       stim = basic_stim(16);
0147 
0148    case 'draeger-eit'
0149      [fr] = read_draeger_header( fname );
0150      % Currently Draeger equipment uses this pattern with 5mA injection
0151      stim = mk_stim_patterns(16,1,[0,1],[0,1],{'rotate_meas','no_meas_current'},.005);
0152      [stim(:).framerate] = deal(fr);
0153      [vvOrig, tstamp, event, signals, counter] = read_draeger_file( fname );
0154      % Based on the draeger *get file converter
0155      ft = [.00098242, .00019607];% estimated AA: 2016-04-07
0156      vv = ft(1)*vvOrig(1:208,:) - ft(2)*vvOrig(322+(1:208),:);
0157      auxdata.vv = vvOrig; % the actual voltages
0158      auxdata.tstamp = tstamp; % timestamp in number of days
0159 
0160      T_int = mean(diff(tstamp))*24*3600; % in S
0161      auxdata.frame_rate = fr; % use assigned fr
0162      auxdata.event = event;
0163      auxdata.signals = signals;
0164      % Based on correlating and scaling with outputs of *.get file
0165      % Note that this does not 100% fit but is the best possible guess so far
0166      fc = 194326.3536;
0167      auxdata.curr = vvOrig(224+(1:16), :) / fc;
0168      fv = 0.11771;
0169      auxdata.volt = 1/fv * (vvOrig(256+(1:16), :) - vvOrig(578+(1:16), :));
0170      auxdata.counter = counter;
0171 
0172    case {'raw', 'sheffield'}
0173       fmt = 'sheffield';
0174       vv = sheffield_readdata( fname );
0175 
0176       stim = basic_stim(16);
0177    case {'p2k', 'its'}
0178       fmt = 'its';
0179       vv = its_readdata( fname );
0180 
0181       stim = 'UNKNOWN';
0182    case {'txt','iirc'}
0183       fmt = 'iirc';
0184       vv = iirc_readdata( fname );
0185 
0186       stim = basic_stim(16);
0187    case 'uct_seq'
0188       [vv,auxdata] = UCT_sequence_file( fname );
0189 
0190       stim = 'UNKNOWN';
0191    case 'uct_cal'
0192       [vv,auxdata] = UCT_calibration_file( fname );
0193 
0194       stim = 'UNKNOWN';
0195    case 'uct_data'
0196       [vv] = UCT_LoadDataFrame( fname );
0197 
0198       stim = 'UNKNOWN';
0199    case {'carefusion','goeiimf-eit'}
0200       fmt = 'goeiimf-eit';
0201       [vv, auxdata_out, auxtime] = carefusion_eit_readdata( fname );
0202       auxdata.auxdata = auxdata_out;
0203       auxdata.auxtime = auxtime;
0204   
0205       stim = mk_stim_patterns(16,1,[0,1],[0,1], {'no_meas_current','rotate_meas'}, .005);
0206 
0207    case 'lq1'
0208       [vv] = landquart1_readdata( fname );
0209 
0210       stim = mk_stim_patterns(32,1,[0,5],[0,5],{'no_rotate_meas','no_meas_current'},.005);
0211       
0212    case {'lq2','lq3'}
0213       fmt = 'lq2';
0214       [vv, elecImps, tStampsAbs, tStampsRel] = landquart2_readdata( fname );
0215       auxdata.elec_impedance = elecImps;
0216       auxdata.t_abs = tStampsAbs;
0217       auxdata.t_rel = tStampsRel;
0218       auxdata.frame_rate = 1e6/median(diff(tStampsRel)); %median in case there are jumps
0219 
0220       stim = mk_stim_patterns(32,1,[0,5],[0,5],{'no_rotate_meas','no_meas_current'},.005);
0221      
0222    case 'lq4pre'
0223       [vv] = landquart4pre_readdata( fname );
0224 
0225       stim = mk_stim_patterns(32,1,[0,5],[0,5],{'no_rotate_meas','no_meas_current'},.005);
0226 
0227    case {'lq4','lq5'}
0228       [vv, evtlist, elecImps, tStampsAbs, tStampsRel, n_elec, amp, skip, extra, patPosition] = landquart4_readdata( fname );
0229       auxdata.event = evtlist;
0230       auxdata.elec_impedance = elecImps;
0231       auxdata.t_abs = tStampsAbs;
0232       auxdata.t_rel = tStampsRel;
0233       auxdata.pat_position = patPosition;
0234       for fn = fieldnames(extra)' % copy extra's fields into auxdata
0235          auxdata.(fn{1}) = extra.(fn{1});
0236       end
0237 
0238       [stim, auxdata.meas_sel] = mk_stim_patterns(n_elec,1,[0,skip+1],[0,skip+1],{'no_rotate_meas','no_meas_current'},amp);
0239       
0240    case 'dixtal_encode'
0241       [vv] = dixtal_read_codepage( fname );
0242       stim = 'N/A';
0243 
0244    case 'dixtal'
0245       [vv] = dixtal_read_data( fname, frame_range, extra );
0246       auxdata = vv(1025:end,:);
0247       vv      = vv([1:1024],:);
0248       stim= mk_stim_patterns(32,1,[0,4],[0,4], ... 
0249               {'meas_current','no_rotate_meas'}, 1);
0250 
0251    case 'sciospec-eit'
0252       %% TODO: very early code for 16 electrode
0253       %% Not very robust
0254       auxdata = sciospec_eit_loader(fname);
0255       d=cat(3,auxdata(:).data);
0256       d=d(:,1:16,1:end);
0257       vv=reshape(d,256,[]);
0258 
0259    case 'vtm';
0260      auxdata = read_vtm(fname);
0261 %  Using flip is easier with than the code in read_vtm
0262 %    flip = kron(ones(16,1),kron([-1;1],ones(6,1)));
0263 %    vv = flip.*dd;
0264      vv = auxdata.EITdata192;
0265      stim = mk_stim_patterns(16,1,[0,8],[0,1],{'rotate_meas'});
0266 
0267 % It's also possible to use EIDORSdata192, with no_rotate_meas
0268 %    vv = auxdata.EIDORSdata192;
0269 %    stim = mk_stim_patterns(16,1,[0,8],[0,1],{'no_rotate_meas'});
0270 
0271    otherwise
0272       error('eidors_readdata: file "%s" format unknown', format);
0273 end
0274 if ~iscell(auxdata)
0275    auxdata.format = fmt;
0276    if ~isfield(auxdata,'frame_rate');
0277       auxdata.frame_rate = NaN; % NaN for unknown
0278    end
0279 end
0280 
0281 
0282 function check_file_exists(fname)
0283    if ~exist(fname,'file')
0284       error([fname,' does not exist']);
0285    end
0286 
0287 % HDF5 format.
0288 %  vv is the voltages
0289 %  auxdata is all the data
0290 %  stim is reconstructed stim
0291 % Fname : datasetname
0292 function [vv, auxdata, stim]= h5_readdata( fname, datasetname );
0293    version = h5read(fname, '/VERSION');
0294    if version<2023.4; error('HDF5 version too old'); end
0295    try
0296       info = h5info(fname,'/data');
0297    catch
0298       error('HDF5 files must have /data group');
0299    end
0300 
0301    Data = info.Groups;
0302 
0303    if strcmp(datasetname,'{:all:}')
0304       for ff =  1:length(Data)
0305         [vv{ff},auxdata{ff},stim{ff}] = h5_one_dataset(fname, Data(ff)); 
0306       end
0307       return
0308    end
0309    
0310    if length(Data)==1;
0311       ff=1;
0312    else  %% TODO, figure it out
0313       dsn = {Data(:).Name};
0314       ff = find(strcmp(dsn,['/data/',datasetname]));
0315       if isempty(ff);
0316          error(['Need to specify which dataset to load. [', ...
0317                strjoin(dsn, ' , '),'] are available']); 
0318       end
0319    end
0320    [vv,auxdata,stim] = h5_one_dataset(fname, Data(ff)); 
0321 
0322 function [vv,auxdata,stim] = h5_one_dataset(fname, Data); 
0323    auxdata= struct();
0324    assert( strcmp(Data.Name(1:6), '/data/') );
0325    for i=1:length(Data.Datasets)
0326       name = Data.Datasets(i).Name;
0327       data = h5read(fname, [Data.Name, '/', name]);
0328       outname = ['auxdata.',name];
0329       outname( (outname == ')') | (outname == '(') ) = '_';
0330       eval([outname,'=data;']); % Use eval: we want dots to turn into subfields: aux.Meas.V.Abs
0331    end
0332    auxdata.DataSet = Data.Name(7:end);
0333 
0334    try
0335       vv = auxdata.Meas.V;
0336    catch
0337       error('h5 file has no Meas.V field');
0338    end
0339    if isfield(vv,'Real') 
0340       if isfield(vv,'Imag') 
0341          vv = vv.Real + 1j*vv.Imag;
0342       else
0343          vv = vv.Real;
0344       end
0345    elseif isfield(vv,'Abs') 
0346       vv = vv.Abs;
0347    else
0348       error('no Abs, Real or Imag field');
0349    end
0350 
0351    try
0352       protocol = [Data.Name,'/protocol'];
0353       info = h5info(fname,protocol);
0354    catch
0355       error('HDF5 files must have /data/{Dataset}/protocol group');
0356    end
0357 
0358    % Tested that Datasets are in sort order
0359    Meas=[]; Stim=[];
0360    for i=1:length(info.Datasets)
0361       dsn= info.Datasets(i).Name;
0362       ni=  regexp(dsn,'Stim.I.(?<n>\d+)\(A\)','names');
0363       if ~isempty(ni)
0364          Stim(str2num(ni.n),:) = h5read(fname, [protocol, '/', dsn]);
0365       end
0366       ni=  regexp(dsn,'Meas.V.(?<n>\d+)\(V\)','names');
0367       if ~isempty(ni)
0368          Meas(str2num(ni.n),:) = h5read(fname, [protocol, '/', dsn]);
0369       end
0370    end
0371    if size(Meas,2) ~= size(Stim,2);
0372       error('inconsistent frame size');
0373    end
0374    if size(Meas,1) > size(Stim,1);
0375       Stim(size(Meas,1),1) = 0; % increase size
0376    end
0377    if size(Meas,1) < size(Stim,1);
0378       Meas(size(Stim,1),1) = 0; % increase size
0379    end
0380    Nstims = size(Meas,2);
0381    stim=repmat( ...
0382           struct('stim_pattern',[],'meas_pattern',[]), ...
0383           1,Nstims);
0384    for i=1:Nstims
0385       stim(i) = struct('stim_pattern',Stim(:,i), ...
0386                        'meas_pattern',Meas(:,i).');
0387    end
0388 
0389 function stim = basic_stim(N_el);
0390    stim= mk_stim_patterns(16,1,[0,1],[0,1], ... \
0391           {'no_meas_current','no_rotate_meas'}, 1);
0392 
0393 function fmt = pre_proc_spec_fmt( format, fname );
0394    fmt= lower(format);
0395    if strcmp(fmt,'get')
0396       if is_get_file_a_draeger_file( fname)
0397          fmt= 'draeger-get';
0398       else
0399          fmt= 'mceit';
0400       end
0401    end
0402 
0403    if strcmp(fmt,'get-raw')
0404       fmt= 'mceit';
0405    end
0406    
0407 
0408    if strcmp(fmt,'eit')
0409       draeger =   is_eit_file_a_draeger_file( fname );
0410       swisstom=   is_eit_file_a_swisstom_file( fname );
0411       carefusion= is_eit_file_a_carefusion_file( fname );
0412       if carefusion
0413          eidors_msg('"%s" appears to be in GOEIIMF/Carefusion format',fname,3);
0414          fmt= 'carefusion';
0415       elseif draeger
0416          eidors_msg('"%s" appears to be in Draeger format',fname,3);
0417          fmt= 'draeger-eit';
0418       elseif swisstom
0419          fmt= sprintf('lq%d',swisstom);
0420          if swisstom == 3.5
0421             fmt= 'lq4pre';
0422          end
0423          eidors_msg('"%s" appears to be in %s format',fname,upper(fmt),3);
0424       else
0425          error('EIT file specified, but it doesn''t seem to be a Carefusion file')
0426       end
0427    end
0428 
0429 function df= is_get_file_a_draeger_file( fname)
0430    fid= fopen(fname,'rb');
0431    d= fread(fid,[1 26],'uchar');
0432    fclose(fid);
0433    df = all(d == '---Draeger EIT-Software---');
0434 
0435 function df= is_eit_file_a_draeger_file( fname );
0436    fid= fopen(fname,'rb');
0437    d= fread(fid,[1 80],'uchar');
0438    fclose(fid);
0439    ff = strfind(char(d), '---Draeger EIT-Software');
0440    if ff;
0441       df = 1;
0442       eidors_msg('Draeger format: %s', d(ff(1)+(0:30)),4);
0443    else 
0444       df = 0;
0445    end
0446 
0447 function df= is_eit_file_a_swisstom_file( fname );
0448    fid= fopen(fname,'rb');
0449    d= fread(fid,[1 4],'uchar');
0450    fclose(fid);
0451    
0452 % Note that there are two possible LQ4 formats, either with
0453 %     d=[0,0,0,4] or with d = [4,0,0,0]
0454 %  we call the first one version 3.5
0455    if any(d(4)==[2,3,4]) && all(d([1,2,3]) == 0);
0456       df = d(4);
0457       if d(4)==4; df = 3.5; end
0458 %     eidors_msg('Swisstom format: LQ%d', df, 4);
0459    elseif all(d(1:4) == [4,0,0,0])
0460       df = 4;
0461    elseif all(d(1:4) == [5,0,0,0])
0462       df = 5;
0463    else 
0464       df = 0;
0465    end
0466 
0467 function df = is_eit_file_a_carefusion_file( fname )
0468    fid= fopen(fname,'rb');
0469    d= fread(fid,[1 80],'uchar');
0470    fclose(fid);
0471    df = 0;
0472    if d(1:2) ~= [255,254]; return; end
0473    if d(4:2:end) ~= 0; return; end
0474    str= char(d(3:2:end));
0475    tests1= '<?xml version="1.0" ?>';
0476    tests2= '<EIT_File>';
0477    if ~any(strfind(str,tests1)); return; end
0478    if ~any(strfind(str,tests2)); return; end
0479    
0480    df= 1;
0481    return
0482 
0483 function [vv, auxdata, auxtime] = carefusion_eit_readdata( fname );
0484    fid= fopen(fname,'rb');
0485    d= fread(fid,[1 180],'uchar');
0486    str= char(d(3:2:end));
0487    outv = regexp(str,'<Header>(\d+) bytes</Header>','tokens');
0488    if length(outv) ~=1; 
0489       error('format problem reading carefusion eit files');
0490    end
0491    header = str2num(outv{1}{1});
0492 
0493 % NOTE: we're throwing away the header completely. This isn't
0494 % the 'right' way to read a file.
0495 
0496    fseek(fid, header, -1); % From the beginning
0497    d_EIT= fread(fid,[inf],'float32');
0498 
0499    fseek(fid, header, -1); % From the beginning
0500    d_struct= fread(fid,[inf],'int32');  % In the struct, meaning of every int: 1 Type, 2 Waveform channel No., 3 sequence No., 4 size in byte, 5 count
0501    fclose(fid);
0502    pt_EIT=1;               % pointer of the EIT data
0503    vv=[];
0504    auxdata=[];
0505    
0506    while (pt_EIT<=length(d_struct))
0507       switch d_struct(pt_EIT)
0508          case 3, % type=3,  EIT transimpedance measurement
0509            % d_struct(pt_EIT+2), Sequence No., start from 0
0510            vv(1:256,1+d_struct(pt_EIT+2))= d_EIT(pt_EIT+6:pt_EIT+6+255);
0511            pt_EIT=pt_EIT+6+256;
0512          case 7, %type=7, EIT image vector
0513            pt_EIT=pt_EIT+912+6;        % drop the image
0514          case 8, %type=8, Waveform data
0515             aux_seg_len = d_struct(pt_EIT+3)/4*d_struct(pt_EIT+4);
0516            if (d_struct(pt_EIT+1) == 0)     % waveform channel 0 is the analog input
0517                aux_segment = d_EIT(pt_EIT+6 : pt_EIT+6+aux_seg_len-1);
0518                auxdata = [auxdata; aux_segment];
0519            else
0520                % drop all other waveform channels (see Waves/Waveform in header file for what they are)
0521            end  
0522            pt_EIT=pt_EIT+6+aux_seg_len;           
0523          case 10,% type = 10, unknown type
0524            %drop
0525            pt_EIT=pt_EIT+6+d_struct(pt_EIT+3)/4*d_struct(pt_EIT+4);
0526          otherwise,
0527            eidors_msg('WARNING: unknown type in carefusion file type');
0528            pt_EIT=pt_EIT+6+d_struct(pt_EIT+3)/4*d_struct(pt_EIT+4);
0529        end
0530    end
0531    vv=vv(47+[1:208],:);
0532    auxtime = (cumsum([1 13*ones(1,15)])-1)/208;             % relative time as fractions of the EIT frame: assuming uniform sampling - not sure this is 100% correct(!)
0533    auxtime = reshape(repmat(1:size(vv,2), length(auxtime), 1),[],1) + repmat(auxtime, 1, size(vv,2))';
0534    
0535 
0536 % Read the old <2006 Draeger "get" file format
0537 function [vv,curr,volt,auxdata] = draeger_get_readdata( fname );
0538    fid= fopen(fname,'rb');
0539    emptyctr=0;
0540    while emptyctr<2
0541      line= fgetl(fid);
0542      if isempty(line)
0543         emptyctr= emptyctr+1;
0544      else
0545         eidors_msg('data not understood',0);
0546         emptyctr=0;
0547      end
0548    end
0549    d= fread(fid,inf,'float',0,'ieee-le');
0550    fclose(fid);
0551 
0552    if rem( length(d), 256) ~=0
0553       eidors_msg('File length strange - cropping file',0);
0554       d=d(1:floor(length(d)/256)*256);
0555    end
0556 
0557    dd= reshape( d, 256, length(d)/256);
0558    vv= untwist(dd);
0559 
0560    curr=0.00512*dd(209:224,:);  % Amps
0561    volt=12*dd(225:240,:); %Vrms
0562    auxdata= dd(241:255,:);
0563    auxdata= auxdata(:);
0564     
0565 function jnk
0566 %[adler@adler01 sept07]$ xxd Sch_Pneumoperitoneum_01_001.get | head -30
0567 %0000000: 2d2d 2d44 7261 6567 6572 2045 4954 2d53  ---Draeger EIT-S
0568 %0000010: 6f66 7477 6172 652d 2d2d 0d0a 2d2d 2d50  oftware---..---P
0569 %0000020: 726f 746f 636f 6c20 4461 7461 2d2d 2d2d  rotocol Data----
0570 %0000030: 2d0d 0a0d 0a44 6174 653a 2020 3138 2d30  -....Date:  18-0
0571 %0000040: 322d 3230 3034 0d0a 5469 6d65 3a20 2031  2-2004..Time:  1
0572 %0000050: 333a 3138 2050 4d0d 0a0d 0a46 696c 656e  3:18 PM....Filen
0573 %0000060: 616d 653a 2020 2020 2020 2020 2053 6368  ame:         Sch
0574 %0000070: 5f50 6e65 756d 6f70 6572 6974 6f6e 6575  _Pneumoperitoneu
0575 %0000080: 6d5f 3031 5f30 3031 2e67 6574 0d0a 4453  m_01_001.get..DS
0576 %0000090: 5020 5365 7269 616c 204e 722e 3a20 2020  P Serial Nr.:
0577 %00000a0: 4549 5430 322f 3035 2f30 3030 360d 0a0d  EIT02/05/0006...
0578 %00000b0: 0a46 7265 7175 656e 6379 205b 487a 5d3a  .Frequency [Hz]:
0579 %00000c0: 2020 2020 3937 3635 362e 330d 0a47 6169      97656.3..Gai
0580 %00000d0: 6e3a 2020 2020 2020 2020 2020 2020 2020  n:
0581 %00000e0: 2020 2031 3134 350d 0a41 6d70 6c69 7475     1145..Amplitu
0582 %00000f0: 6465 3a20 2020 2020 2020 2020 2020 2031  de:            1
0583 %0000100: 3030 300d 0a53 616d 706c 6520 5261 7465  000..Sample Rate
0584 %0000110: 205b 6b48 7a5d 3a20 2020 2035 3030 300d   [kHz]:    5000.
0585 %0000120: 0a50 6572 696f 6473 3a20 2020 2020 2020  .Periods:
0586 %0000130: 2020 2020 2020 2020 2032 300d 0a46 7261           20..Fra
0587 %0000140: 6d65 733a 2020 2020 2020 2020 2020 2020  mes:
0588 %0000150: 2020 2020 2031 300d 0a0d 0a0d 0a
0589 %                                       8bc33e       10........>
0590 %0000160: 3f 0548633e bf20933d e192393d a568ea  ?.Hc>. .=..9=.h.
0591 %0000170: 3c 2530f63c 27e6043d 2043ad3c 25ce93  <%0.<'..= C.<%..
0592 %0000180: 3c aebcce3c 8714643d e3533d3e 65b6e1  <...<..d=.S=>e..
0593 %0000190: 3e 7a62103f 81c4143e 8c99813d 35921d  >zb.?...>...=5..
0594 %00001a0: 3d 8e6b0d3d 690cf93c 9910713c 3c9289  =.k.=i..<..q<<..
0595 %00001b0: 3c f6736f3c 2291453d ad1cab3d 386f15  <.so<".E=...=8o.
0596 %00001c0: 3e e82a143f 2a952e3f f568493e f08a8c  >.*.?*..?.hI>...
0597 %00001d0: 3d e43e0e3d 7040253d 19f4af3c 67fd93  =.>.=p@%=...<g..
0598 
0599 function vv = sheffield_readdata( fname );
0600    fid=fopen(fname,'rb','ieee-le');
0601    draw=fread(fid,[104,inf],'float32');
0602    fclose(fid);
0603    ldat = size(draw,2);
0604 
0605    [x,y]= meshgrid( 1:16, 1:16);
0606    idxm= y-x;
0607  % HW gain table
0608    gtbl = [0,849,213,87,45,28,21,19,21,28,45,87,213,849,0];
0609    idxm(idxm<=0)= 1;
0610    idxm= gtbl(idxm);
0611 
0612  % Multiply by gains
0613    draw = draw .* (idxm(idxm>0) * ones(1,ldat)); 
0614 
0615    vv= zeros(16*16, size(draw,2));
0616    vv(find(idxm),:) = draw;
0617 
0618  % Add reciprocal measurements
0619    vv= reshape(vv,[16,16,ldat]);
0620    vv= vv + permute( vv, [2,1,3]);
0621    vv= reshape(vv,[16*16,ldat]);
0622 
0623    
0624 
0625 function [vv,curr,volt,auxdata,auxtime,rawdata] = mceit_readdata( fname );
0626 
0627    fid= fopen(fname,'rb');
0628    d= fread(fid,inf,'float');
0629    fclose(fid);
0630 
0631    if rem( length(d), 256) ~=0
0632       eidors_msg('File length strange - cropping file',0);
0633       d=d(1:floor(length(d)/256)*256);
0634    end
0635 
0636    dd= reshape( d, 256, length(d)/256);
0637    rawdata = dd;
0638    no_reciprocity = (dd(39,1)==0);      %104 measurements per frame
0639    if no_reciprocity
0640        dd=transfer104_208(dd);
0641    end
0642    vv= untwist(dd);
0643 
0644    curr=0.00512*dd(209:224,:);  % Amps
0645    volt=12*dd(225:240,:); %Vrms
0646    auxdata= dd(241:256,:);
0647    auxdata= auxdata(:);
0648    %input impedance=voltage./current-440;        Ohm
0649    
0650    if no_reciprocity
0651      % remove every 14th, 15th and 16th sample as they are always zero
0652      auxdata(14:16:end) = []; auxdata(14:15:end) = []; auxdata(14:14:end) = [];
0653      % only 104 measurements were performed with NON-UNIFORM sampling of AUX in between EIT samples
0654      % Sampling procedure was thought to be: 13 AUX1 13 AUX2 12 AUX3 11 AUX4 10 AUX5 9 ... 2 AUX13 1 [END OF FRAME]
0655 %      auxtime = (cumsum([1 13 12:-1:2]) - 1) / 104;
0656      % However, this seems to be a little different, finally the sampling points were determined
0657      % empirically by measuring a sawtooth signal (linear ramp) and fitting the timings
0658      auxtime = [0.0000,0.1390,0.2535,0.3617,0.4642,0.5486,0.6367,0.7110,0.7780,0.8354,0.8865,0.9304,0.9711];  % mean fit at 25 Hz
0659 %      auxtime = [0.0000,0.1460,0.2789,0.3895,0.4875,0.5788,0.6608,0.7356,0.7986,0.8569,0.9045,0.9454,0.9824,]; % mean fit at 44 Hz
0660    else
0661      % all 208 measurements were performed
0662      auxtime = (cumsum([1 13*ones(1,15)])-1)/208;             % relative time as fractions of the EIT frame
0663    end
0664    % time of AUX signal relative to frame number (starting with 1 - Matlab style)
0665    auxtime = reshape(repmat(1:size(dd,2), length(auxtime), 1),[],1) + repmat(auxtime, 1, size(dd,2))';
0666    
0667 
0668 function array208=transfer104_208(array104),
0669 % The order in the 208-vector is
0670 % Index   Inject    Measure Pair
0671 % 1         1-2        3-4        (alternate pair is 3-4, 1-2, at location 39)
0672 % 2         1-2        4-5
0673 % 3         1-2        5-6
0674 % ¡­
0675 % 13        1-2        15-16
0676 % 14        2-3        4-5
0677 % 15        2-3        5-6
0678 % ¡­
0679 % 26        2-3        16-1
0680 % 27        3-4        5-6
0681 % ¡­
0682 % 39        3-4        1-2
0683 % 40        4-5        6-7
0684 % ¡­
0685 
0686 ind=[39,51,63,75,87,99,111,123,135,147,159,171,183,52,64,76, ...
0687      88,100,112,124,136,148,160,172,184,196,65,77,89,101,113, ...
0688      125,137,149,161,173,185,197,78,90,102,114,126,138,150,162, ...
0689      174,186,198,91,103,115,127,139,151,163,175,187,199,104,116, ...
0690      128,140,152,164,176,188,200,117,129,141,153,165,177,189, ...
0691      201,130,142,154,166,178,190,202,143,155,167,179,191,203, ...
0692      156,168,180,192,204,169,181,193,205,182,194,206,195,207,208];
0693 ro=[1:13, 14:26, 27:38, 40:50, 53:62, 66:74, 79:86, 92:98, ...
0694     105:110,118:122,131:134,144:146,157:158,170:170];
0695 
0696 [x,y]=size(array104);
0697 if x~=256 && y~=256,
0698     eidors_msg(['eidors_readdata: expectingin an input array ', ...
0699                 'of size 208*n']);
0700     return;
0701 elseif y==256,
0702     array104=array104';
0703     y=x;
0704 end
0705 array208=array104;
0706 for i=1:y,
0707     array208(ind,i)=array104(ro,i);    
0708 end
0709    
0710 
0711 % measurements rotate with stimulation, we untwist them
0712 function vv= untwist(dd);
0713    elec=16;
0714    pos_i= [0,1];
0715    ELS= rem(rem(0:elec^2-1,elec) - ...
0716         floor((0:elec^2-1)/elec)+elec,elec)';
0717    ELS=~any(rem( elec+[-1 0 [-1 0]+pos_i*[-1;1] ] ,elec)' ...
0718         *ones(1,elec^2)==ones(4,1)*ELS')';
0719    twist= [               0+(1:13), ...
0720                          13+(1:13), ...
0721            39-(0:-1:0),  26+(1:12), ...
0722            52-(1:-1:0),  39+(1:11), ...
0723            65-(2:-1:0),  52+(1:10), ...
0724            78-(3:-1:0),  65+(1: 9), ...
0725            91-(4:-1:0),  78+(1: 8), ...
0726           104-(5:-1:0),  91+(1: 7), ...
0727           117-(6:-1:0), 104+(1: 6), ...
0728           130-(7:-1:0), 117+(1: 5), ...
0729           143-(8:-1:0), 130+(1: 4), ...
0730           156-(9:-1:0), 143+(1: 3), ...
0731           169-(10:-1:0),156+(1: 2), ...
0732           182-(11:-1:0),169+(1: 1), ...
0733           195-(12:-1:0), ...
0734           208-(12:-1:0) ];
0735     vv= zeros(256,size(dd,2));
0736     vv(ELS,:)= dd(twist,:);
0737    %vv= dd(1:208,:);
0738 
0739 % Read data from p2k files (I T S system)
0740 % FIXME: this code is very rough, it works for
0741 %   only eight ring data records
0742 function  vv = its_readdata( fname ) 
0743    fid= fopen( fname, 'rb', 'ieee-le');
0744    vv=[];
0745 
0746    % don't know how to interpret header
0747    header= fread(fid, 880, 'uchar');
0748    frameno= 0;
0749    rings= 8;
0750    while( ~feof(fid) )
0751        frameno= frameno+1;
0752        % don't know how to interpret frame header
0753        framehdr= fread(fid, 40);
0754        data= fread(fid, 104*rings, 'double');
0755        vv= [vv, data];
0756    end
0757 
0758    if 0 % convert a ring to 208
0759       ringno= 1;
0760       ld= size(vv,2);
0761       vx= [zeros(1,ld);vv( ringno*104 + (-103:0) ,: )];
0762       idx= ptr104_208;
0763       vv= vx(idx+1,:);
0764    end
0765 
0766 % pointer to convert from 104 to 208 meas patterns
0767 function idx= ptr104_208;
0768     idx= zeros(16);
0769     idx(1,:)= [0,0,1:13,0];
0770     ofs= 13;
0771     for i= 2:14
0772         mm= 15-i;
0773         idx(i,:) = [zeros(1,i+1), (1:mm)+ofs ];
0774         ofs= ofs+ mm;
0775     end
0776 
0777     idx= idx + idx';
0778     
0779 function vv = iirc_readdata( fname );
0780     fid= fopen( fname, 'r');
0781     while ~feof(fid)
0782        line = fgetl(fid);
0783        if isempty(line)
0784            continue;
0785        end
0786        
0787        num= regexp(line,'Channel : (\d+)');
0788        if ~isempty(num)
0789            channels= str2num( line(num(1):end ) );
0790            continue;
0791        end
0792        
0793        num= regexp(line,'Frequency : (\d+)kHz');
0794        if ~isempty(num)
0795            freqency= str2num( line(num(1):end ) );
0796            continue;
0797        end
0798 
0799        num= regexp(line,'Scan Method : (\w+)');
0800        if ~isempty(num)
0801            scan_method=  line(num(1):end );
0802            continue;
0803        end
0804 
0805        num= regexp(line,'Study : (\w+)');
0806        if ~isempty(num)
0807            study=  line(num(1):end);
0808            continue;
0809        end
0810            
0811        if strcmp(line,'Data');
0812            data= fscanf(fid,'%f',[4,inf])';
0813            continue;
0814        end
0815     end
0816     vv= data(:,1) + 1i*data(:,2);
0817     if length(vv) ~= channels^2
0818         error('eidors_readdata: data length wrong')
0819     end
0820 
0821 % stimulation is the fwd_model stimulation data structure
0822 % meas_select indicates if data is NOT measures on current electrode
0823 function [stimulations,meas_select] = UCT_sequence_file( fname );
0824    % (c) Tim Long
0825    % 21 January 2005
0826    % University of Cape Town
0827 
0828 
0829    % open the file
0830    fid = fopen(fname, 'rt');
0831 
0832    % check to see if file opened ok
0833    if fid == -1
0834          errordlg('File not found','ERROR')  
0835          return;
0836    end
0837 
0838 
0839    tline = fgetl(fid);             % get the spacer at top of text file
0840 
0841    % the measurement and injection pattern is stored as follows:
0842    % I1V1:  db #$00,#$0F,#$00,#$00,#$10,#$00,#$00,#$21,#$00 ...
0843    % I2V2:  db #$11,#$0F,#$11,#$11,#$10,#$11,#$11,#$21,#$11 ...
0844    % etc
0845    % need to put all the bytes in a vector
0846 
0847    % tokenlist will store the list of bytes as strings
0848    tokenlist = [];
0849 
0850 
0851    tline = fgetl(fid);             % get first line of data
0852 
0853    while length(tline) ~= 0
0854        
0855        % the first few characters in the line are junk
0856        rem = tline(11:end);            % extract only useful data
0857        
0858        % extract each byte
0859        while length(rem) ~=0
0860            [token, rem] = strtok(rem, ',');
0861            tokenlist = [tokenlist; token];
0862        end
0863        
0864        % get the next line in sequence file
0865        tline = fgetl(fid);
0866    end
0867 
0868    fclose(fid);
0869 
0870    % got everything in string form... need to covert to number format
0871 
0872    drive_lay = [];
0873    drive_elec = [];
0874    sense_lay = [];
0875 
0876    injection_no = 1;
0877    % for each injection
0878    for i=1:3:length(tokenlist)
0879        
0880        % get injection layer
0881        tsource_layer = tokenlist(i,3);
0882        tsink_layer = tokenlist(i,4);
0883        source_layer = sscanf(tsource_layer, '%x');
0884        sink_layer = sscanf(tsink_layer, '%x');
0885        
0886        drive_lay = [drive_lay; [source_layer sink_layer]];
0887          
0888        
0889        % get drive pair
0890        tsource_elec = tokenlist(i+1,3);
0891        tsink_elec = tokenlist(i+1,4);
0892        source_elec = sscanf(tsource_elec, '%x');
0893        sink_elec = sscanf(tsink_elec, '%x');
0894        
0895        drive_elec = [drive_elec; [source_elec sink_elec]];
0896        
0897        
0898        % get sense layer pair
0899        tpos_layer = tokenlist(i+2,3);
0900        tneg_layer = tokenlist(i+2,4);
0901        pos_sense_layer = sscanf(tpos_layer, '%x');
0902        neg_sense_layer = sscanf(tneg_layer, '%x');
0903        
0904        sense_lay = [sense_lay; [pos_sense_layer neg_sense_layer]];
0905    end
0906 
0907    n_elec = size(sense_lay,1);
0908    n_inj  = size(drive_lay,1);       % every injection
0909    elecs_per_plane = 16; % FIXED FOR THE UCT DEVICE
0910    raw_index = 0;
0911    meas_select = [];
0912 
0913    for i=1:n_inj      % for every injection
0914        stimulations(i).stimulation= 'mA';
0915        
0916        % find the injection electrodes
0917        e_inj_p = drive_lay(i, 1) * elecs_per_plane + drive_elec(i,1) + 1;
0918        e_inj_n = drive_lay(i, 2) * elecs_per_plane + drive_elec(i,2) + 1;
0919 
0920        % create the stimulation pattern for this injection
0921        inj = zeros(n_elec, 1);
0922        inj(e_inj_p) = 1;
0923        inj(e_inj_n) = -1;
0924        stimulations(i).stim_pattern = inj;
0925      
0926        % the UCT instrument always makes 16 measurements per injection.
0927        % the +ve and -ve electrodes are always adjacent, but might be on
0928        % different planes
0929        meas_pat = [];
0930        for e = 0:15
0931            raw_index = raw_index + 1;
0932            meas = zeros(1, n_elec);   % the measurement electrodes for this sample
0933 
0934            % find the measurement electrodes for this measurement (+ve elec is
0935            % next to -ve electrode)
0936            e_meas_p = sense_lay(i,1) * elecs_per_plane + mod(e+1,elecs_per_plane) + 1;
0937            e_meas_n = sense_lay(i,2) * elecs_per_plane + e + 1;
0938 
0939            % if either of the drive electrodes are equal to any of the sense
0940            % electrodes, we must not include this sample
0941            if any( e_meas_p == [e_inj_p, e_inj_n] ) | ...
0942               any( e_meas_n == [e_inj_p, e_inj_n] ) 
0943                continue;
0944            end
0945 
0946            meas(e_meas_p) = -1;
0947            meas(e_meas_n) = 1;
0948            meas_select = [meas_select;raw_index];
0949            
0950            % add this measurement to the measurement pattern
0951            meas_pat = [meas_pat; meas];
0952 
0953        end     % for each injection there are actually only 13-16 measurements
0954        stimulations(i).meas_pattern = sparse(meas_pat);
0955    end
0956 
0957 function [cur_data,no_cur_data] = UCT_calibration_file( fname );
0958    fid = fopen(fname, 'rb');
0959    mag_num = fread(fid, 1, 'int32');
0960    version = fread(fid, 1, 'int32');
0961 %  comments = fread(fid, 2048, 'char'); MATLAB CHANGED CHAR to 2 bytes
0962    comments = fread(fid, 2048, 'uint8');
0963    comments = char(comments');
0964    no_of_layers = fread(fid, 1, 'int32');
0965    uppa_dac = fread(fid, 1, 'float64');
0966    lowa_dac = fread(fid, 1, 'float64');
0967    raw_frame_size = fread(fid, 1, 'int32');
0968 
0969    no_cur_data = [];
0970    cur_data = [];
0971 
0972    for i=1:no_of_layers
0973        no_cur_data = [no_cur_data;fread(fid, raw_frame_size, 'float64')];
0974        cur_data = [cur_data;fread(fid, raw_frame_size, 'float64')];
0975    end
0976    fclose(fid);
0977 
0978 %  no_cur_data = UCT_ShuffleData( no_cur_data);
0979 %  cur_data =    UCT_ShuffleData( cur_data);
0980 
0981 function [v_raw] = UCT_LoadDataFrame(infilename)
0982 % [v_raw] = LoadDataFrame(infilename)
0983 %
0984 % Loads the data from the tomography file
0985 %
0986 % (c) Tim Long
0987 % 21 January 2005
0988 % University of Cape Town
0989 
0990 
0991 % Open the file
0992 fid = fopen(infilename, 'rb');
0993 
0994 if fid == -1
0995       errordlg('File not found','ERROR')  
0996       return;
0997 end
0998 
0999 %%% new file changes
1000 magic_number = fread(fid, 1, 'uint32');
1001 
1002 if magic_number ~= 2290649224
1003    disp('UCT File: wrong file type'); 
1004 end
1005 version = fread(fid, 1, 'uint32');
1006 foffset = fread(fid, 1, 'uint32');
1007 no_of_layers = fread(fid, 1, 'uint32');
1008 frame_size = fread(fid, 1, 'uint32');
1009 fseek(fid, foffset-8, 'cof');
1010 
1011 %%% end of new file changes
1012 %%% old file stuff
1013 % no_of_layers = fread(fid, 1, 'uint32');
1014 % frame_size = fread(fid, 1, 'uint32');
1015 
1016 frame_no = fread(fid, 1, 'uint32');
1017 
1018 v_raw = [];
1019 while feof(fid)==0
1020     v_raw = [v_raw, fread(fid, frame_size*no_of_layers, 'float64')];
1021     frame_no = fread(fid, 1, 'uint32');
1022 end
1023 
1024 fclose(fid);
1025 % UCT_ShuffleData??
1026 
1027 % Read data from the file format develped by Pascal
1028 % Gaggero at CSEM Landquart in Switzerland.
1029 function [vv] = landquart1_readdata( fname );
1030    FrameSize = 1024;
1031 
1032    enableCounter=1;
1033    counter=0;
1034 
1035    fid=fopen(fname);
1036 
1037    codec=fread(fid,1,'int'); %read codec
1038    if codec~=1; error('Codec unexpected value'); end
1039 
1040    nbFileInHeader=fread(fid,1,'int'); %read codec
1041    
1042    for i=1:nbFileInHeader
1043        lenghtFile = fread(fid,1,'int64'); 
1044        jnk= fread(fid,lenghtFile,'int8');% discard the file
1045    end
1046    
1047    
1048    vv= []; 
1049    while 1; 
1050        type=fread(fid,1,'int'); %check if not End of file
1051        if type == -1; break ; end
1052        
1053        nel= fread(fid,1,'int');
1054        sz= fread(fid,1,'int');
1055        iq=fread(fid,[2,sz/8],'int')';
1056 
1057        vv = [vv, iq*[1;1j]]; %I+ 1j*Q vectors
1058        
1059        for j=1:nel-1
1060            sz= fread(fid,1,'int');
1061            jnk = fread(fid,sz/4,'int');
1062        end
1063        
1064    end
1065    
1066 % Read data from the file format develped by Swisstom, Landquart, Switzerland.
1067 function [vv, elecImps, tStampsAbs, tStampsRel] = landquart2_readdata( fname )
1068    [fid msg]= fopen(fname,'r','ieee-be','UTF-8');
1069    try
1070       format_version = fread(fid,1,'int32','ieee-be');
1071       if format_version ~= 3
1072          error('unsupported file format version');
1073       else
1074          % get expected number of frames and preallocate variables
1075          fseek(fid,16,'cof');
1076          nFrames = fread(fid, 1, 'int32', 'ieee-be');
1077          tStampsAbs = nan(1, nFrames);
1078          tStampsRel = nan(1, nFrames);
1079          viPayload = nan(64, nFrames);
1080          iqPayload = nan(2048, nFrames);
1081          
1082          progress_msg('Reading EIT frame', 0, nFrames);
1083 
1084          header_size = 2264; 
1085          fseek(fid,header_size + 8,'bof');
1086          
1087          %%% get frame size and length of payload at end of frame
1088          frame_length = fread(fid, 1, 'int32', 'ieee-be') + 12;
1089          %%% move back to start of 1. frame
1090          fseek(fid,header_size,'bof');
1091          
1092          %%% Read frames
1093          i = 1;
1094          % check there is 'frame_length' ahead in the file
1095          while fseek(fid, frame_length,'cof') ~= -1
1096             if mod(i, 100) == 0
1097                 progress_msg(i, nFrames);
1098             end
1099             
1100             fseek(fid, -frame_length,'cof');
1101             % read absolute timestamp (milliseconds resolution)
1102             tStampsAbs(i) = fread(fid,1,'int64','ieee-be');
1103             pl = fread(fid,1,'int32','ieee-le');    % payload length (i.e. frame length)
1104             frame_header = fread(fid,15,'int32','ieee-le'); 
1105             % read relative timestamp (microseconds resolution)
1106             tStampsRel(i) = frame_header(5);    
1107             
1108             % drop some unintersting parts
1109             fseek(fid,12, 'cof');
1110             
1111             % get electrode impedance measurement
1112             viPayload(:,i) = fread(fid,64,'int32','ieee-le');
1113             % get effective payload (voltage measurements)
1114             iqPayload(:,i) = fread(fid,2048,'int32','ieee-le');
1115             fseek(fid,header_size + i*frame_length,'bof');
1116             i = i +1;
1117          end
1118       end
1119    catch err
1120       fclose(fid);
1121       rethrow(err);
1122    end
1123    fclose(fid);
1124    
1125    progress_msg(inf);
1126    
1127    i = i-1;
1128    if i ~= nFrames       
1129       % remove data which were preallocated but not read
1130       if i < nFrames       
1131          tStampsAbs(i+1:end) = [];
1132          tStampsRel(i+1:end) = [];
1133          viPayload(:,i+1:end) = [];
1134          iqPayload(:,i+1:end) = [];
1135       end
1136       eidors_msg('"%s": expected %.0f frames but read %.0f',fname, nFrames, i ,3);
1137    end
1138    
1139    % convert into matlab date format, e.g. read via datestr(tStampAbs(1))
1140    tStampsAbs = tStampsAbs/(24*3600*1E3) + datenum(1970, 1, 1, 1, 0, 0);
1141    % strictly speaking we would need to correct for DST, but as this is only
1142    % supported by ML R2014b or higher we avoid it to be backwards compatible
1143 %    tStampsAbs = tStampsAbs + isdst(datetime(datevec(tStampsAbs(1)), 'TimeZone', 'local'))/24;   % add one hour if DST
1144    
1145    % this is just a simple guess
1146    amplitudeFactor = 2.048 / (2^20 * 360 * 1000);
1147    vv = amplitudeFactor * (iqPayload(1:2:end,:) + 1i*iqPayload(2:2:end,:));
1148    
1149    elecImps = viPayload(1:2:end,:) + 1i*viPayload(2:2:end,:);
1150 
1151 % Read data from the file format develped by Swisstom, Landquart, Switzerland.
1152 
1153 % notes from Beat about format_version==5
1154 % There are only small changes, main structure is the same, but with these changes:
1155 % - header-size: 2672
1156 % - each frame starts with 16 bytes: timestamp of packet reception (8 bytes), event Type (4 bytes), payload size (4 bytes).
1157 %
1158 % See attached screenshot below: Red corner marks the start of the first frame. With the red arrow, the event type is marked, and the red cursor is on the payload size. After those 16 bytes, the eit-frame follows as it was for version 4.
1159 %
1160 % Beside different header size, the following loop is needed while reading in a frame:
1161 % check event Type:
1162 %                 If type != 0, then ignore payload and move by payload size to the next frame.
1163 %                 Otherwise, frame-payload starts
1164 function [vv, evtlist, elecImps, tStampsAbs, tStampsRel, n_elec, amp, skip, extra, patPosition] = landquart4_readdata( fname )
1165    evtlist = [];
1166    [fid msg]= fopen(fname,'r','ieee-le','UTF-8');
1167    try
1168       format_version = fread(fid,1,'int32','ieee-le');
1169       switch format_version
1170          case {4,5}; % OK, we can process
1171          otherwise;
1172             error('unsupported file format version');
1173       end
1174          header_size = fread(fid,1,'int32', 'ieee-le');   
1175          eit_frame_offset = 328; % 60 + 12 + 256 = 328
1176          iq_payload = 2048;
1177          vi_payload = 64;
1178          pat_position = 3; 
1179          
1180          % get expected number of frames and preallocate variables
1181          fseek(fid,16,'cof');
1182          nFrames = fread(fid, 1, 'int32', 'ieee-le');
1183          tStampsAbs = nan(1, nFrames);
1184          tStampsRel = nan(1, nFrames);
1185          viPayload = nan(vi_payload, nFrames);
1186          iqPayload = nan(iq_payload, nFrames);
1187          patPosition = nan(pat_position, nFrames);
1188 
1189          progress_msg('Reading EIT frame', 0, nFrames);
1190          fseek(fid,424,'bof');
1191          x = fread(fid, 2, 'int16', 'ieee-le');
1192          extra.filename = fread(fid,[1,100], 'int16=>char', 'ieee-le');
1193          extra.conditions = fread(fid,[1,300], 'int16=>char', 'ieee-le');
1194          extra.comments = fread(fid,[1,600], 'int16=>char', 'ieee-le');
1195 
1196          %fseek(fid,2428,'bof');
1197          extra.frame_rate = fread(fid, 1, 'float', 'ieee-le');
1198          amp = fread(fid, 1, 'float', 'ieee-le');
1199          extra.freq = fread(fid, 1, 'float', 'ieee-le');
1200          extra.settle = fread(fid, 1, 'float', 'ieee-le');
1201 
1202          %fseek(fid,2444,'bof');
1203          skip = fread(fid, 1, 'int16', 'ieee-le');
1204          x = fread(fid, 1, 'int16', 'ieee-le');
1205          n_elec = fread(fid, 1, 'int16', 'ieee-le');
1206 
1207          fseek(fid,header_size,'bof');
1208          
1209          %%% Read frames
1210          i = 1;
1211          evti = 1;
1212          % Check there is at least 'eit_frame_offset' ahead in file
1213          while fseek(fid, eit_frame_offset,'cof') ~= -1
1214             if mod(i, 100) == 0
1215                 progress_msg(i, nFrames);
1216             end
1217             
1218             fseek(fid, -eit_frame_offset,'cof');
1219             % drop frame header and some payload:
1220             % read absolute timestamp (milliseconds resolution)
1221             tStampsAbs(i) = fread(fid,1,'int64','ieee-le'); 
1222             ft = fread(fid,1,'int32','ieee-le'); 
1223             pl = fread(fid,1,'int32','ieee-le');
1224             if ft == 1
1225                % event
1226                evtlist(evti).timestamp = tStampsAbs(i) ;
1227                evtlist(evti).frame = i ;
1228                evti = evti + 1;
1229                if pl > 0
1230                   evtlist(evti).eventId = fread(fid, 1, 'int32', 'ieee-le');
1231                 end
1232             elseif ft == 0
1233                 frame_header = fread(fid,15,'int32','ieee-le'); 
1234                 % read relative timestamp (microseconds resolution)
1235                 tStampsRel(i) = frame_header(5);   
1236                
1237                 % get patient position
1238                 patPosition(:, i) = fread(fid,pat_position,'int32','ieee-le');                
1239                 % get electrode impedance measurement
1240                 viPayload(:,i) = fread(fid,vi_payload,'int32','ieee-le');
1241                 % get effective payload (voltage measurements)
1242                 iqPayload(:,i) = fread(fid,iq_payload,'int32','ieee-le');
1243                 fseek(fid,pl-4*iq_payload-eit_frame_offset,'cof');
1244                 i = i+1;
1245             elseif pl > 0
1246                fseek(fid,pl,'cof'); 
1247             else
1248                % nothing to do
1249             end
1250       end
1251    catch err
1252       fclose(fid);
1253       rethrow(err);
1254    end
1255    fclose(fid);
1256    
1257    progress_msg(inf);
1258    
1259    i = i-1;
1260    if i ~= nFrames       
1261       % remove data which were preallocated but not read
1262       if i < nFrames       
1263          tStampsAbs(i+1:end) = [];
1264          tStampsRel(i+1:end) = [];
1265          viPayload(:,i+1:end) = [];
1266          iqPayload(:,i+1:end) = [];
1267          patPosition(:,i+1:end) = [];
1268       end
1269       eidors_msg('"%s": expected %.0f frames but read %.0f',fname, nFrames, i ,3);
1270    end
1271    
1272    % convert into matlab date format, e.g. read via datestr(tStampAbs(1))
1273    tStampsAbs = tStampsAbs/(24*3600*1E3) + datenum(1, 1, 1, 0, 0, 0);
1274    % strictly speaking we would need to correct for DST, but as this is only
1275    % supported by ML R2014b or higher we avoid it to be backwards compatible
1276 %    tStampsAbs = tStampsAbs + isdst(datetime(datevec(tStampsAbs(1)), 'TimeZone', 'local'))/24;   % add one hour if DST
1277 
1278    % this is just a simple guess
1279    amplitudeFactor = 2.048 / (2^20 * 360 * 1000);
1280    vv = amplitudeFactor * (iqPayload(1:2:end,:) + 1i*iqPayload(2:2:end,:));
1281    
1282    %% elecImps depends on a factor which varies with the SBC version,
1283    %%   however the following value is close for most versions
1284    impedanceFactor =  2.048 / (2^12 * 0.173 * 0.003) / 2^15; % = 0.9633 / 2^15;
1285    elecImps = impedanceFactor * (viPayload(1:2:end,:) + 1i*viPayload(2:2:end,:));
1286 
1287 function [vv] = landquart4pre_readdata( fname )
1288    [fid msg]= fopen(fname,'r','ieee-le','UTF-8');
1289    try
1290       format_version = fread(fid,1,'int32','ieee-le');
1291       if format_version ~= 4
1292          error('unsupported file format version');
1293       else
1294          header_size = fread(fid,1,'int32', 'ieee-le'); 
1295          
1296          fseek(fid,header_size + 12,'bof');
1297          %%% get frame size and length of payload at end of frame
1298          frame_length = fread(fid, 1, 'int32', 'ieee-le') + 16;
1299          %%% move back to start of 1. frame
1300          fseek(fid,header_size,'bof');
1301          
1302          %%% Read frames
1303          i = 1;
1304          while fseek(fid, 1,'cof') ~= -1
1305             fseek(fid, -1,'cof');
1306             % drop frame header and some payload:
1307             % (16 + 60) + 12 + 256 = 344
1308             fseek(fid,344, 'cof');
1309             iqPayload(:,i) = fread(fid,2048,'int32','ieee-le');
1310             fseek(fid,header_size + i*frame_length,'bof');
1311             i = i +1;
1312          end
1313 
1314       end
1315    catch err
1316       fclose(fid);
1317       rethrow(err);
1318    end
1319    fclose(fid);
1320 
1321    % this is just a simple guess
1322    amplitudeFactor = 2.048 / (2^20 * 360 * 1000);
1323    vv = amplitudeFactor * (iqPayload(1:2:end,:) + 1i*iqPayload(2:2:end,:));
1324    
1325 %  Output: [encodepage] = eidors_readdata( path,'DIX_encode');
1326 %   where path= '/path/to/Criptografa_New.dll' (provided with system)
1327 function [vv] = dixtal_read_codepage( fname );
1328    %Fname must be the Criptografa_New dll
1329    fid = fopen(fname,'rb');
1330    b1234= fread(fid, [1,4], 'uint8');
1331    if b1234~= [77, 90, 144, 0];
1332       error('This does not appear to be the correct Dll');
1333    end
1334    fseek(fid, hex2dec('e00'), 'bof');
1335    encodepage1 = fread(fid, hex2dec('1000'), 'uint8');
1336    fclose(fid);
1337    encodepage1 = flipud(reshape(encodepage1,4,[]));
1338    vv = encodepage1(:);
1339 
1340 %        format = 'DIXTAL'
1341 %           - output: [vv] = eidors_readdata( fname, 'DIXTAL', encodepage );
1342 function [vv] = dixtal_read_data( file_name, frame_range, encodepage1 );
1343 
1344    % Get encodepage from the file
1345    fid=fopen(file_name,'rb');
1346    n1 = fgetl(fid);
1347    n2 = fgetl(fid);
1348    n3 = fgetl(fid); 
1349    nframes = str2num( n3 );
1350    
1351    fseek(fid, hex2dec('800'), 'bof');
1352    encodepage2 = fread(fid, hex2dec('1000'), 'uint8');
1353    fclose(fid);
1354 
1355    encodepage = bitxor(encodepage1, encodepage2);
1356 
1357    % Read the file
1358    dplen = hex2dec('1090');
1359    start = hex2dec('2800');
1360    fid=fopen(file_name,'rb');
1361    if isempty(frame_range)
1362       frame_range = 0:nframes;
1363    else
1364       frame_range = frame_range - 1;
1365       frame_range(frame_range>nframes) = [];
1366    end      
1367    vv= zeros(dplen/4, length(frame_range));
1368    k= 1;
1369    for i = frame_range
1370       status= fseek(fid, start + i*dplen, 'bof');
1371       if status~=0; break; end
1372       datapage = fread(fid, dplen, 'uint8');
1373       if length(datapage)== 0; 
1374          vv= vv(:,1:k-1); break; % frame number inside file is wrong
1375       end
1376       vv(:,k) = proc_dixtal_data( datapage, encodepage );
1377       k=k+1;
1378    end
1379    fclose(fid);
1380 
1381 
1382 function  data = proc_dixtal_data( b, encodepage );
1383 
1384    cryptolen = 1024*4;
1385    b(1:cryptolen) = bitxor(b(1:cryptolen), encodepage);
1386 
1387    b1 = b(1:4:end,:); %b1 = bitand(b1,127);
1388    b2 = b(2:4:end,:);
1389    b3 = b(3:4:end,:);
1390    b4 = b(4:4:end,:);
1391    bb  = [b1,b2,b3,b4]*(256.^[3;2;1;0]);
1392 
1393    sgnbit = bitand( bb,  2^31);
1394    sgnbit = +1*(sgnbit == 0) - 1*(sgnbit == 2^31);
1395 
1396    expbit = bitand( bb, 255*2^23 ) /2^23; 
1397    expbit = 2.^(expbit - 127);
1398 
1399    fracbt = bitand( bb, 2^23-1)/2^23 + 1;
1400    data = sgnbit .* expbit .* fracbt;
1401 
1402 function [fr] = read_draeger_header( filename );
1403 %READ_DRAEGER_HEADER   Opens and reads the header portion of the Draeger .eit
1404 %file. Current parameter returned is frame rate.
1405 %
1406 % function [fr,s] = ReadDraegerHeader(filename)
1407 % fr        double      scalar      frame rate in hertz
1408 
1409 % Determine file version
1410 K0 = 2048;   % bytes to read in char class
1411 K1='Framerate [Hz]:'; % Draeger frame rate line in header
1412 K2=15;                % Length of K1 string
1413 K3=13;                % Following F1 Field for delimiter (look for ^M)
1414 
1415 % Open file for reading in little-endian form
1416 fid = fopen(filename,'r','l');
1417 if fid == -1
1418     error('Error read_draeger_header: can''t open file');
1419 end
1420 header = fread(fid,K0,'*char')';
1421 index = strfind(header,K1);
1422 if ~isempty(index)
1423     [tok,rem]= strtok(header(index+K2:end),K3);
1424     fr = str2num(tok); 
1425 else
1426     error('Error read_draeger_header: frame rate unspecified');
1427 end
1428 
1429 if isempty(fr)
1430     error('Error read_draeger_header: Frame rate could not be read');
1431 end
1432 
1433 fclose(fid);
1434 
1435 
1436 function [vd, tstamp, event, signals, counter] = read_draeger_file(filename)
1437 %READDRAEGERFILE   Opens and reads a Draeger .eit file. Returns an array
1438 %containing the voltage data arranged per frame.
1439 %
1440 % Input:
1441 % filename  char        1xN
1442 %
1443 % Output:
1444 % vd        double      MxN         M = volt measurements per frame (mV)
1445 %                                   N = number of frames
1446 
1447 
1448    ff = dir(filename);
1449    filelen = ff.bytes;
1450 
1451    % Open file for reading in little-endian form
1452    fid = fopen(filename,'r','l');
1453    if fid == -1
1454        error('Error read_draeger_file: file could not be opened');
1455    end
1456 
1457    % Determine file version
1458    K0 = 128;   % bytes to read in char class
1459 
1460    header = fread(fid,K0,'*char')';
1461 if 0 % THis analysis doesn't seem useful
1462    if ~isempty(strfind(header,'V3.2'))
1463        version = '3.2';
1464    elseif ~isempty(strfind(header,'V4.01'))
1465        version = '4.01'; 
1466    elseif ~isempty(strfind(header,'V1.10')) % 2014!!
1467        version = '1.10'; 
1468    else
1469        error('Error read_draeger_file: unknown file version');
1470    end
1471 end
1472 
1473    % Read Format version
1474    fseek(fid,0,-1); % beginning of file
1475    hdr = fread(fid, 8, 'uint8');
1476    offset1 =  (256.^(0:3))*hdr(5:8); % "good" data starts at offset #2
1477    switch sprintf('%02X-',header(1:4));
1478       case '1F-00-00-00-';
1479          type='Draeger v31';  SPC = 4112;
1480       case '20-00-00-00-';
1481          type='Draeger v32';  SPC = 5200;
1482       case '33-00-00-00-';
1483          type='Draeger v51';  SPC = 5495;
1484       otherwise;
1485          error('File "%s" is format version %d, which we don''t know how to read', ...
1486                hdr(1))
1487    end
1488 
1489    len = floor( (filelen-offset1)/SPC );
1490    vd= zeros(600,len);
1491    tstamp = zeros(1, len, 'double'); 
1492    counter = zeros(1, len, 'uint16');     
1493    event = repmat(char(0), 30, len);
1494    rawsigs = zeros(67, len, 'single');
1495    ss= offset1 + (0:len)*SPC;
1496 
1497    progress_msg('Reading EIT frame', 0, length(ss)-1);
1498    for k= 1:length(ss)-1;
1499        if mod(k, 100) == 0
1500            progress_msg(k, length(ss)-1);
1501        end
1502        % for newest Dräger files: one frame has 5495 bytes
1503        % - [0001-0008] 8 bytes of whatever in first frame then zero
1504        % - [0009-5168] 5160 bytes of signals: as 645 float64 (double)
1505        % - [5169-5436] 268 bytes of "medibus" signals: as 67 float32 (single)
1506        % - [5437-5467] 30 bytes of "event" string: as 30 chars
1507        % - [5468-5491] ??? some bytes (mostly int16 and zeros) of we do not know what
1508        % - [5492-5493] 2 bytes of some counter signal which gets reset from time to time: as uint16
1509        % - [5494-5495] 2 bytes of zeros
1510        
1511       if fseek(fid,ss(k)+8,-1)<0; break; end
1512       % timestamp in number of days
1513       tstamp(k) = fread(fid,1,'float64', 0, 'ieee-le'); 
1514       
1515       if fseek(fid,ss(k)+16,-1)<0; break; end
1516       vd(:,k)=fread(fid,600,'double', 0, 'ieee-le');
1517       
1518       start = 16+600*8;
1519       if fseek(fid,ss(k)+start,-1)<0; break; end
1520       gugus(:,k) = fread(fid, floor((5169-start)/8), 'double', 'ieee-le');
1521       
1522       if fseek(fid,ss(k)+5169,-1)<0; break; end
1523 %       rawsigs(:,k) = fread(fid, 52, 'single', 'ieee-le');
1524       rawsigs(:,k) = fread(fid, 67, 'single', 'ieee-le');
1525       
1526       if fseek(fid,ss(k)+5437,-1)<0; break; end
1527       event(:, k) = fread(fid, 30, 'char', 'ieee-le');
1528             
1529 %       fseek(fid,ss(k)+5437+30,-1);
1530 %       dada(:,k) = fread(fid, floor(SPC-(5437+30)/8), 'double', 'ieee-le');
1531 
1532       if fseek(fid,ss(k)+5491,-1)<0; break; end
1533       counter(k) = fread(fid, 1, 'uint16', 'ieee-le');      
1534    end
1535    
1536    % convert raw signals into a structure
1537    % WARNING: this list is far from correct or complete
1538    % TODO: verify and fix this!
1539    SignalNames = {'RTD_Paw_mbar_01', 'RTD_Paw_mbar_02', 'RTD_Flow_L_min_03', 'RTD_Flow_L_min_04', 'RTD_Vol__mL_05', 'RTD_Vol__mL_06', 'unknown_07', 'unknown_08', 'unknown_09', 'unknown_10', 'unknown_11', 'Tispon_sec_12', 'Pmin_mbar_13', 'P0_1_mbar_14', 'Pmean_mbar_15', 'unknown_16', 'PEEP_mbar_17', 'unknown_18', 'unknown_19', 'unknown_20', 'unknown_21', 'PIP_mbar_22', 'unknown_23', 'unknown_24', 'unknown_25', 'unknown_26', 'VT_mL_27', 'unknown_28', 'unknown_29', 'unknown_30', 'unknown_31', 'VT_mL_32', 'VTispon_mL_33', 'NIF_mbar_34', 'MVleak_L_min_35', 'unknown_36', 'RRspon_1_min_37', 'unknown_38', 'MVspon_L_min_39', 'unknown_40', 'MVspon_L_min_41', 'unknown_42', 'RSB_1_LXMin_43', 'RRspon_1_min_44', 'unknown_45', 'unknown_46', 'unknown_47', 'unknown_48', 'unknown_49', 'etCO2___50', 'etCO2___51', 'unknown_52', 'unknown_53', 'unknown_54', 'unknown_55', 'unknown_56', 'unknown_57', 'unknown_58', 'unknown_59', 'unknown_60', 'unknown_61', 'unknown_62', 'unknown_63', 'unknown_64', 'unknown_65', 'unknown_66', 'unknown_67'};
1540 %    SignalNames = {'RTD_Paw_mbar_01', 'RTD_Paw_mbar_02', 'unknown_03', 'RTD_Flow_L_min_04', 'RTD_Vol__mL_05', 'RTD_Vol__mL_06', 'Cdyn_mL_mbar_07', 'unknown_08', 'unknown_09', 'R_mbar_L_s_10', 'unknown_11', 'unknown_12', 'unknown_13', 'unknown_14', 'unknown_15', 'unknown_16', 'unknown_17', 'unknown_18', 'unknown_19', 'unknown_20', 'unknown_21', 'unknown_22', 'unknown_23', 'unknown_24', 'unknown_25', 'unknown_26', 'unknown_27', 'unknown_28', 'unknown_29', 'unknown_30', 'unknown_31', 'VT_mL_32', 'unknown_33', 'unknown_34', 'MVleak_L_min_35', 'unknown_36', 'unknown_37', 'unknown_38', 'unknown_39', 'unknown_40', 'MV_L_min_41', 'unknown_42', 'unknown_43', 'RR_1_min_44', 'unknown_45', 'unknown_46', 'unknown_47', 'unknown_48', 'unknown_49', 'etCO2___50', 'etCO2_kPa_51', 'etCO2_mmHg_52', 'unknown_53', 'unknown_54', 'unknown_55', 'unknown_56', 'unknown_57', 'unknown_58', 'unknown_59', 'unknown_60', 'unknown_61', 'unknown_62', 'unknown_63', 'unknown_64', 'unknown_65', 'unknown_66', 'unknown_67'};
1541     for iS = 1 : size(rawsigs,1) 
1542 %         signals.(['sig_',num2str(iS)]) = rawsigs(iS,:);
1543         SigTmp = rawsigs(iS,:);
1544         % set invalid values to NAN (on time signals < -1E38, others -1000)
1545         SigTmp((SigTmp == -1000) | (SigTmp < -1E38)) = nan;
1546         signals.(SignalNames{iS}) = SigTmp;
1547     end
1548    
1549    % End of function
1550    progress_msg(inf);
1551    fclose(fid);
1552 
1553 % Call with the name of the folder. It looks for all eit files
1554 function s=sciospec_eit_loader(fname);
1555   if isdir(fname);
1556     files = dir([fname,'/*.eit']);
1557     nfiles= length(files);
1558     progress_msg('Sciospec Load:', 0,nfiles);
1559     for i=nfiles:-1:1
1560        progress_msg(nfiles-i+1,nfiles);
1561        s(i) = sciospec_frame([fname,'/',files(i).name]);
1562     end
1563     progress_msg(inf); 
1564   else
1565     error 'not dir'
1566   end
1567     
1568 
1569 
1570 function s=sciospec_frame(fname);
1571   fid = fopen(fname,'r');
1572   indata = 0;
1573   s.inj=[];
1574   s.data=[];
1575   while(~feof(fid))
1576     line = fgetl(fid);
1577     if indata
1578       num = textscan(line, '%f', 'Delimiter',' ');
1579       num = cell2mat(num)';
1580       if rem(indata,2)==1; %odd is inj
1581          s.inj(end+1,:) = num;
1582       else
1583          s.data(end+1,:) = num(1:2:end) + 1j*num(2:2:end);
1584       end
1585       indata = indata+1;
1586     end
1587     if length(line)>20 && strcmp(line(1:20),'MeasurementChannels:');
1588        num = textscan(line(21:end), '%f', 'Delimiter',',');
1589        s.MeasurementChannels=cell2mat(num);
1590     end
1591     if length(line)>51 && strcmp(line(1:51),'MeasurementChannelsIndependentFromInjectionPattern:');
1592        num = textscan(line(52:end), '%f', 'Delimiter',',');
1593        s.MeasurementChannelsIndependentFromInjectionPattern=cell2mat(num);
1594        indata = 1;
1595     end
1596   end
1597   fclose(fid);
1598 
1599 % Contribution by Zhanqi Zhao <zhanqizhao@gzhmu.edu.cn>
1600 function [vtmData] = read_vtm(full_dfile_name)
1601 %READ_vtm
1602 %   ´Ë´¦ÏÔʾÏêϸ˵Ã÷
1603     [~,~,file_ext]= fileparts(full_dfile_name);
1604     if  strcmp(file_ext,'.vtm') %.bin File, image data
1605         
1606         len_header  = 333  ;  % 2+204+2+125 = 333
1607         len_data    = 530  ;  % 2+528
1608         len_Medibus = 242  ;  % 2+60*4
1609         len_imginfo = 102  ;  % 2+25*4
1610         len_frame = len_data + len_Medibus +len_imginfo;
1611         try
1612             file_inf = dir(full_dfile_name);
1613             len = file_inf.bytes;
1614             framenums=(len-len_header-2)/len_frame;
1615             framenums=int32(framenums);
1616         end
1617         fid=fopen(full_dfile_name,'r');
1618         Header = fread(fid,len_header,'uint8'); 
1619         if ~strcmp(sprintf('%02X',Header(1:2)),'AA55')
1620             error('δʶ±ðµÄÎļþ');
1621         end
1622 %%
1623             systemdata1 = native2unicode(Header(3:206))';
1624             str1 = regexpi(systemdata1, '([^\t\f\n\r]*)', 'match');
1625 %             systemparm = {'DeviceModel','Frequency','Framerate','Amplitude','Time'};
1626             for i = 1:length(str1)
1627                 parm = regexpi(str1{1,i},':','split');
1628                 parmname = char(parm{1,1});                
1629                 if length(parm)>2
1630                     parmvalue = parm{1,2};
1631                     for j = 3:length(parm)
1632                         parmvalue = [parmvalue,':',parm{1,j}];
1633                     end
1634                 elseif length(parm)<2
1635                     break;
1636                 else
1637                     parmvalue = parm{1,2};
1638                 end
1639                 Systemparm.(parmname) = parmvalue;
1640             end
1641 
1642             patientdata1 = Header(209:209+124)';
1643             Patient.PatientID   = native2unicode(patientdata1(1:20), 'UTF-8');
1644             Patient.Name        = native2unicode(patientdata1(21:50), 'UTF-8');
1645             Patient.Age         = native2unicode(patientdata1(51:53), 'UTF-8');
1646             Patient.Gender      = native2unicode(patientdata1(54:55), 'UTF-8');
1647             Patient.Departments = native2unicode(patientdata1(56:85), 'UTF-8');
1648             Patient.Notes       = native2unicode(patientdata1(86:125), 'UTF-8');
1649 
1650 %%
1651 %         fseek(fid,len_header,-1); % beginning of file
1652 %         dd = fread(fid,[len_frame framenums],'uint8',0,'ieee-le');
1653 %         dd = reshape(d,[],framenums);
1654         dd= zeros(530,framenums); 
1655         ss= len_header + (0:framenums-1)*len_frame;
1656         for k= 1:length(ss);
1657            if fseek(fid,ss(k),-1)<0; break; end
1658            dd(:,k)=fread(fid,530,'uint8', 0, 'ieee-le');
1659         end
1660         ETDdata = dd(3:3+527,:);% 528ԭʼÊý¾Ý
1661         contact0=ETDdata((9+192)*2+1:((9+192)*2+32),:);
1662         for row=1:16
1663             RawDataContact(row,:)=contact0((row-1)*2+2,:)*256+contact0((row-1)*2+1,:); 
1664         end
1665         temp=RawDataContact(1,:);
1666         
1667         for row=1:15
1668             RawDataContact(row,:)=min(RawDataContact(row,:),RawDataContact(row+1,:));
1669         end
1670         RawDataContact(16,:)=min(RawDataContact(16,:),temp);
1671         
1672 %         fseek(fid,len_header,-1); % beginning of file
1673 %         dd = fread(fid,[len_frame/2 framenums],'uint16',0,'ieee-le');
1674 % %         dd = reshape(d,[],framenums);
1675         dd= zeros(192,framenums);
1676         di= zeros(192,framenums);
1677         ss= len_header + 20  + (0:framenums-1)*len_frame;
1678         for k= 1:length(ss);
1679            if fseek(fid,ss(k),-1)<0; break; end
1680 %          di(:,k)=fread(fid,192,'int16', 0, 'ieee-le');
1681            dd(:,k)=fread(fid,192,'uint16', 0, 'ieee-le');
1682         end
1683 %       di(:,end) = []; % Last seems bad
1684         EITdata192 = dd;
1685         EITdata192=  -EITdata192;
1686         for i=1:16
1687            EITdata192(12*i-5:12*i,:)=  -EITdata192(12*i-5:12*i,:);
1688         end
1689         move192_index = [0,0,1,2,3,4,5,6,6,6,7,8,9,10,11,12];
1690         columnIdx192 = bsxfun(@circshift,reshape(1:192,12,16),move192_index);
1691         columnIdx192 = columnIdx192(:);
1692         EIDORSdata192 = EITdata192(columnIdx192,:);% ²»°üº¬ÎÞЧµç¼«µÄ¾ø¶Ôµç¼«Î»ÖÃÅÅÐò
1693         vtmData.dd = dd;
1694  %%
1695 %         fseek(fid,len_header,-1); % beginning of file
1696 %         d = fread(fid,len-335,'float32',0,'ieee-le');%°´ÕÕfloatÀàÐͶÁÈ¡Êý¾ÝΪµ½ÁÐÏòÁ¿d
1697 %         dd = reshape(d,[],framenums);
1698         dd= zeros(60,framenums);
1699         ss= len_header + len_data + 2  + (0:framenums-1)*len_frame;
1700         for k= 1:length(ss);
1701            if fseek(fid,ss(k),-1)<0; break; end
1702            dd(1:52,k)=fread(fid,52,'float32', 0, 'ieee-le');
1703         end
1704         Medibus = dd;
1705 %%
1706         dd= zeros(25,framenums);
1707         ss= len_header + len_data + len_Medibus + 2 + (0:framenums-1)*len_frame;
1708         for k= 1:length(ss);
1709            if fseek(fid,ss(k),-1)<0; break; end
1710            dd(:,k)=fread(fid,25,'float32', 0, 'ieee-le');
1711         end
1712         Imginfo = dd;
1713         fclose(fid);
1714 %%
1715         vtmData.Patient = Patient;
1716         vtmData.Systemparm = Systemparm;        
1717         vtmData.ETDdata = ETDdata;
1718         vtmData.EITdata192 = EITdata192;
1719         vtmData.EIDORSdata192 = EIDORSdata192;
1720         vtmData.Medibus = Medibus;
1721         vtmData.Imginfo = Imginfo;
1722         vtmData.RawDataContact = RawDataContact/58;
1723         % vtmData.RawDataContact = RawDataContact/55.2688;
1724         vtmData.frame_rate = 19.7636;
1725     else
1726         error('File type cannot be recognised');
1727     end
1728 
1729 
1730 function do_unit_test
1731 
1732    % LQ2 test
1733    pth = get_contrib_data('human-ventilation-2014','human-ventilation-2014.eit'); 
1734    [dd,auxdata,stim]= eidors_readdata(pth);
1735    unit_test_cmp('LQ2',dd(100,100), -8.820502983940973e-04 - 8.022670735677084e-04i,1e-16);
1736    unit_test_cmp('LQ2:elecZ', auxdata.elec_impedance(3), -6.883231000000000e+06 - 1.805360000000000e+05i,1e-10);
1737    unit_test_cmp('LQ2:t_rel', auxdata.t_rel(2), 38442185);
1738    unit_test_cmp('LQ2:frame_rate', auxdata.frame_rate, 28.041838422927,1e-10);
1739 
1740    pth = get_contrib_data('et-healthy-breathing','ET_Test_File_Healthy_Lung_01.eit');
1741    [dd,auxdata,stim]= eidors_readdata(pth);
1742    unit_test_cmp('Draeger.eit:frame_rate', auxdata.frame_rate, 20,1e-10);
1743 
1744    pth = get_contrib_data('et-healthy-breathing','ET_Test_File_Healthy_Lung_01_01.get');
1745    [dd,auxdata,stim]= eidors_readdata(pth);
1746    unit_test_cmp('Draeger.get:frame_rate', isnan(auxdata.frame_rate), true);
1747 
1748    pth = write_hdf5_sample;
1749    [dd,~,stim]= eidors_readdata(pth,'h5','datasetEIT');
1750    unit_test_cmp('hdf5:size', size(dd), [208, 34]);
1751    unit_test_cmp('hdf5:stim', size(stim), [1, 208]);
1752    unit_test_cmp('hdf5:data', dd(1:6), [13441   10288   24740   19477  15029 9780]);
1753 
1754    [dd,~,stim]= eidors_readdata(pth,'h5-all');
1755    [dd,~,stim]= eidors_readdata(pth,'hdf5-all');
1756    unit_test_cmp('hdf5:size:1', size(dd{1}), [544,1]);
1757    unit_test_cmp('hdf5:size:2', size(dd{2}), [208,34]);
1758    unit_test_cmp('hdf5:data:1', dd{1}(1:4), [ -3674 -3610 -3871 -4177]');
1759    unit_test_cmp('hdf5:data:2', dd{2}(1:4), [13441   10288   24740   19477]);
1760 
1761    % Midas Med test
1762    pth = get_contrib_data('midas-medical-sample2024','EITData2024_10_31_08_43_09__S01.vtm');
1763    [dd,auxdata,stim]= eidors_readdata(pth);
1764    fmdl = mk_library_model('adult_male_16el');
1765    fmdl.stimulation = stim;
1766    fmdl.frame_rate = auxdata.frame_rate;
1767 %  vs = fwd_solve(mk_image(fmdl,1));
1768 
1769    imdl = select_imdl(fmdl, {'GREIT:NF=0.5 32x32'});
1770    imgr = inv_solve(imdl, mean(dd,2),dd);
1771    imgr.elem_data = subtract_rank(imgr.elem_data,0.95);
1772    [ROIs,imgR] = define_ROIs(imdl,[-2,0,2],[2,0,-2],struct('normalize',0));
1773    sig = -ROIs*imgr.elem_data;
1774    lt = {'LineWidth',2};
1775    subplot(211);
1776    plot(imgr.time,sig',lt{:}); box off; xlim tight;
1777     legend('RV','LV','RD','LD');
1778    vertlines(50 + [1,50]*.5);
1779    subplot(212);
1780    imgr.get_img_data.frame_select = 500 + (1:50)*5;
1781    imgr.calc_colours.ref_level = 0;
1782    imgr.show_slices.img_cols = 17;
1783    imgr.calc_colours.greylev =  .01;
1784    imgr.calc_colours.backgnd =[1,1,1];
1785    show_slices(imgr)

Generated on Sun 29-Dec-2024 11:41:59 by m2html © 2005