Eidors-logo    

EIDORS: Electrical Impedance Tomography and Diffuse Optical Tomography Reconstruction Software

EIDORS (mirror)
Main
Documentation
Tutorials
Install
Contrib Data
GREIT
Browse Docs
Browse SVN

News
Mailing list
(archive)
FAQ
Developer
− Examples
− Structure
− Objects
− Caching
                       

 

Hosted by
SourceForge.net Logo

 

EIDORS: Programming / Examples

Simple Example

Consider the simplest possible EIT system. We have a resistor, and we want to know its value. We therefore attach one electrode to each terminal, and apply several different test currents (I1, I2, I3) and measure the voltages: (V1, V2, V3) From these measurements the a least-squares estimate of the resistance is calculated. This example is available here

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

Forward Solver
    % 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 Solver
    % 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 Calculator
    % 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 (vi−vh) or normalized difference
    • fwd_model.meas_select     measurement reduction (when not measuring on injection electrodes, while given full data set) (vi−vh)/vh?

      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

It is essential that EIDORS be able to cache values that are reused. The design tries to make this as clean as possible, so that the long calculation steps can be sped up without resorting to convoluted code. The design is as follows:
  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.

Last Modified: $Date: 2011-07-15 06:34:46 -0400 (Fri, 15 Jul 2011) $ by $Author: aadler $