EIDORS: Electrical Impedance Tomography and Diffuse Optical Tomography Reconstruction Software

Hosted by

EIDORS: Programming / Examples

Simple Example

Step 1: Define functions for fwd_solve, inv_solve, and calc_jacobian

    % Forward Model:
    % For each stimulation there is I1 into Node1
    %  Node2 is connected to gnd with Zcontact
    %
    % each stim has one measurement pattern where
    %  Vmeas= Meas_pat * Node1
    %       = Meas_pat * I1 * ( Zcontact*2 + R )
    %
    % Thus
    %  V= IR    => [V1;V2;V3] = [I1;I2*I3]*(R + 2*Zcontact)
    function data =f_solve( f_mdl, img )
      R= img.elem_data;

      n_stim= length( f_mdl.stimulation );
      V= zeros(n_stim, 1);

      for i=1:n_stim
        I        = f_mdl.stimulation(i).stim_pattern / 1000; % mA
        meas_pat = f_mdl.stimulation(i).meas_pattern;

        stim_elec= find( I );
        zc       = f_mdl.electrode( stim_elec ).z_contact;

        V(i)= meas_pat * I * ( R + 2*zc);
      end

      data.name='resistor model data';
      data.meas= V;

    % Inverse Model: R= inv(J'*J)*J'*V
    %    This corresponds to the least squares solution
    function img= i_solve( i_mdl, data )
      % Normally the Jacobian depends on an image. Create a dummy one here
      i_img= eidors_obj('image','Unused');
      f_mdl= i_mdl.fwd_model;
      J = calc_jacobian( f_mdl, i_img);

      img.name= 'solved by i_solve';
      img.elem_data= (J'*J)\J'* data.meas;
      img.inv_model= i_mdl;
    

    % Jacobian:      J= dV/dR =I = [I1; I2; I3]
    function J= c_jacobian( f_mdl, img)
      n_stim= length( f_mdl.stimulation );
      J= zeros(n_stim, 1);
      for i=1:n_stim
        J(i)     = f_mdl.stimulation(i).stim_pattern / 1000; % mA
      end
    

Step 2: Create Eidors based code using functions

    %
    % create FEM model structure
    %
    % Fwd model:
    %  Two nodes are in space at [1,1,1] and [2,2,2]
    %  The resistor is connected between them

    r_mdl.name = 'demo resistor model';
    r_mdl.nodes= [1,1,1;  2,2,2];
    r_mdl.elems= [1;2];
    r_mdl.solve=      @f_solve;
    r_mdl.jacobian=   @c_jacobian;

    %
    % create FEM model electrode definitions
    %

    r_mdl.electrode(1).z_contact= 10; % ohms
    r_mdl.electrode(1).nodes=     1;
    r_mdl.gnd_node= 2;

    %
    % create stimulation and measurement patterns
    % patterns are 0.010,0.020,0.030 mA

    for i=1:3
        r_mdl.stimulation(i).stimulation= 'mA';
        r_mdl.stimulation(i).stim_pattern= ( 0.010*i );
        r_mdl.stimulation(i).meas_pattern= 1; % measure electrode 1
    end

    r_mdl= eidors_obj('fwd_model', r_mdl);

    %
    % simulate data for medium with R=1 kohms
    % This medium is called an 'image'
    %

    img_1k = eidors_obj('image', 'homogeneous image', ...
                         'elem_data', 1e3, ...
                         'fwd_model', r_mdl );

    data_1k =fwd_solve( r_mdl, img_1k );

    %
    % add noise to simulated data
    %

    data_noise= eidors_obj('data', 'noisy data', ...
                           'meas', data_1k.meas + 1e-3*randn(3,1));

    %
    % create inverse model
    %

    % create an inv_model structure of name 'demo_inv'
    r_inv.name=  'Resistor Model inverse';
    r_inv.solve= @i_solve;
    r_inv.reconst_type= 'static';
    r_inv.fwd_model= r_mdl;
    r_inv= eidors_obj('inv_model', r_inv);

    %
    % solve inverse model');
    %

    R= inv_solve( r_inv, data_1k );
    fprintf('R calculated with clean data= %5.3f kOhm\n', R.elem_data / 1000 );

    R= inv_solve( r_inv, data_noise );
    fprintf('R calculated with noisy data= %5.3f kOhm\n', R.elem_data / 1000 );
    

Step 3: Run code

    >> resistor_model
        EIDORS:[ fwd_solve: setting cached value ]
        EIDORS:[ inv_solve ]
        EIDORS:[ calc_jacobian: setting cached value ]
    Value calculated with clean data= 1.020 kOhm
        EIDORS:[ inv_solve ]
        EIDORS:[ calc_jacobian: setting cached value ]
    Value calculated with noisy data= 1.068 kOhm
    

EIDORS Objects Reference

• data
  A data object represents one set of measurement data. It is a collection of all measurements for each stimulation pattern. While not simultaneous, we conceptually represent this as representing the conductivity distribution at an instant. It is invisaged that data converter software be written to load from the various hardware systems into this format.
  Properties:
  • data.name     string name of data (or empty string)
  • data.meas     vector (Num_meas × 1) actual measured data, ordered as measurements for each stimulation pattern
  • data.time     scalar measurement start time in seconds since the epoch, 0=unknown. In Matlab this is given by time.
  • data.comment     string comments on this measurement
  • data.meas_config     EIDORS configuration structure meas_config common to all measurements with the same configuration

    Methods: None

• meas_config
  this structure represents a given measurement configuration. This would be created by a data reading function as it loads the data from a file (or gets the input directly from the hardware)
  Properties:
  • meas_config.name     string name of measurement configuration (or empty string)
  • meas_config.units     string measurement units (ie. volts)
  • meas_config.subject     string description of subject (or empty string)
  • meas_config.electrodes     vector (Num_elec × 1) vector of descripions of electrodes, where each electrode has
    • electrode(index).position     vector (x,y,z) position of centre of electrode in units of mm
    • electrode(index).diameter     scalar diameter of electrode in mm
    • electrode(index).diameter     scalar diameter of electrode in mm
    • electrode(index).z_contact     contact impedance (in Ω) may be complex
  • meas_config.stimulation     vector (Num_stim × 1) stimulation patterns stim_model

    Methods: None

• fwd_model
  
  Properties:
  • fwd_model.name     Model name (if known)
  • fwd_model.nodes     position of FEM nodes (Nodes×Dims)
  • fwd_model.elems     definition of FEM elements (Elems×Dims+1)
    Currently defined only for simplex element shapes (i.e. each element had dimentions + 1 nodes)
  • fwd_model.boundary     nodes of element faces on the medium surface
  • fwd_model.gnd_node     Number of node connected to ground
  • fwd_model.electrode     Vector (Num_elecs ×1) of electrode models (elec_model)
  • fwd_model.stimulation     Vector (Num_Stim ×1) of stimulation patterns (stim_model) (current in EIT)
  • fwd_model.dimention     2D, 3D, etc.
  • fwd_model.normalize_measurements     Do we do difference (v_i−v_h) or normalized difference
  • fwd_model.meas_select     measurement reduction (when not measuring on injection electrodes, while given full data set) (v_i−v_h)/v_h?

    Methods:

  • fwd_model.solve     Calculate data object:
    usage: data = fwd_solve( fwd_model, image )
  • fwd_model.jacobian     usage: data = jacobian( fwd_model, image )
  • fwd_model.image_prior     usage: data = image_prior( fwd_model, image )
  • fwd_model.data_prior     usage: data = data_prior( fwd_model, image )
  Questions:
  − No definition of conductivities in model?
  − We need an approach to cache 'hard to calculate' parameters for some of these methods

• elec_model
  
  Properties:
  • elec_model.name     Electrode name (optional)
  • elec_model.z_contact     contact impedance (in Ω) may be complex
  • elec_model.nodes     nodes to which this electrode is attached

    Methods: None

• stim_model
  
  model of a stimulation pattern and accociated measurements
  Properties:
  • stim_model.name     Stimulation name (optional)
  • stim_model.stimulation     Quantity stimulated (mA) in EIT, light in DOT
  • stim_model.stim_pattern     Quantity of stimulation on each electrode (Num_elecs×1) in units of .stimulation
    
  • stim_model.meas_pattern     Measurements pattern for this stimulation (Num_meas×Num_elecs). This is a sparse matrix of the contribution of each electrode to the measurement.

    Methods: None

• inv_model
  
  Note all properties are required for all inv_models
  Properties:
  • inv_model.name     Model name (if known)
  • inv_model.hyperparameter.func     function to call to set hyperparameter value
  • inv_model.hyperparameter.value     specified value of hyperparameter (if inv_model.hyperparameter.func doesn't exist)
  • inv_model.image_prior.func     function to calculate image prior
  • inv_model.image_prior.parameters     parameters to inv_model.image_prior.func
  • inv_model.term_tolerance     termination tolerance (or array of parameters)
  • inv_model.iterations    
  • inv_model.type     'differential' or 'static'
  • inv_model.fwd_model     pointer to fwd_model

    Methods:

  • inv_model.solve     Calculate image object:
    usage: image = inv_solve( model_static, data )
    or image = inv_solve( model_diff, data_1, data_2 )
  Questions:
  −

• image
  
  Properties:
  • image.name     name of image (optional)
  • image.elem_data     data for each element
  • image.type     real, complex, tensor?
  • image.fwd_model     pointer to fwd_model
  • image.inv_model     pointer to inv_model

    Methods: None

  Questions:
  −

Caching

1. Caching should be 'natural'

   This part of the 'overloaded' accessor functions, so for example,

   calc_image_prior used to be
   

          image_prior= feval( inv_model.image_prior.func, inv_model);
   

   now it is (using the eidors_obj function):

          image_prior = eidors_obj('cache', inv_model, 'image_prior');
   
          if isempty(image_prior)
              image_prior= feval( inv_model.image_prior.func, inv_model);
              eidors_obj('cache', inv_model, 'image_prior', image_prior);
          end
   

   so this means that the function 'pointer' in inv_model.image_prior.func = 'np_calc_image_prior' doesn't need to know anything about possible caching.

2. Cached values should not appear when the underlying model has changed.

   This is ensured by creating an 'object repository' using the eidors_obj function. eidors objects now must be constructed using this function, either as

          demo_inv.name= 'Nick Polydorides EIT inverse';
          demo_inv.solve=       'np_inv_solve';
          demo_inv.hyperparameter= 1e-8;
          demo_inv.image_prior.func= 'np_calc_image_prior';
          demo_inv= eidors_obj('inv_model', demo_inv);
   

   or as

          demo_inv= eidors_obj( ...
               'inv_model', 'Nick Polydorides EIT inverse',...
               'solve',          'np_inv_solve', ...
               'hyperparameter', 1e-8, ...
               'func',           'np_calc_image_prior');
   

   whenever an 'object' is modified, such as

          eidors_obj('set', demo_inv, 'solve', 'NEW_SOLVER_CODE' );
   

   then all cached values are flushed.