EIDORS: Electrical Impedance Tomography and Diffuse Optical Tomography Reconstruction Software

Hosted by

EIDORS data structures for forward models

2D foward models

% Create 2D model $Id: tutorial022a.m 3850 2013-04-16 18:13:39Z aadler $
clear mdl elec stim

nn= 84;     % number of nodes
ww=4;       % width = 4
conduc= 1;  % conductivity in Ohm-meters
mdl= eidors_obj('fwd_model','2D rectangle');
mdl= mdl_normalize( mdl, 0);
mdl.nodes = [floor( (0:nn-1)/ww );rem(0:nn-1,ww)]';
mdl.elems = delaunayn(mdl.nodes);
mdl.gnd_node = 1;
elec(1).nodes= [1:ww];      elec(1).z_contact= 0.01;
elec(2).nodes= nn-[0:ww-1]; elec(2).z_contact= 0.01;
stim.stim_pattern= [-1;1];
stim.meas_pattern= [-1,1];
mdl.stimulation= stim;
mdl.electrode= elec;
mdl.solve = @fwd_solve_1st_order;
mdl.system_mat = @system_mat_1st_order;

show_fem(mdl); axis('equal'); set(gca,'Ylim',[-.5,ww-.5]);
print_convert tutorial022a.png '-density 75'

Figure: this shows a single line representing the resistor in 3D. Not very interesting, but this is a simple example.

Defining stimulation patterns

% Create stimulation patterns and solve fwd_model
% $Id: tutorial020b.m 3850 2013-04-16 18:13:39Z aadler $

% Define stimulation patterns
for i=1:3
    r_mdl.stimulation(i).stimulation= 'Amps';
    r_mdl.stimulation(i).stim_pattern= ( 0.001*i );
    r_mdl.stimulation(i).meas_pattern= 1; % measure electrode 1
end

r_mdl.solve=      @tutorial020_f_solve;

% Define an 'image'
resistor = eidors_obj('image', 'resistor');
resistor.elem_data= 1000;
resistor.fwd_model= r_mdl;

% Calculate data for 1k resistor
data_1k0 =fwd_solve( resistor );

% Now change resistor to be 1.2k
resistor.elem_data= 1200;
data_1k2 =fwd_solve( resistor );

>> disp(data_1k0)
    name: 'resistor model data'
    meas: [3x1 double]
    type: 'data'

>> disp(data_1k0.meas')
    1.0200    2.0400    3.0600

Forward solver function

The forward solver used here is tutorial020_f_solve, shown below:

function data =tutorial020_f_solve( img )
% Forward Model for a resistor
% 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)

  R= img.elem_data;
  stim = img.fwd_model.stimulation;

  n_stim= length( stim );
  V= zeros(n_stim, 1);

  for i=1:n_stim
    I        = stim(i).stim_pattern;
    meas_pat = stim(i).meas_pattern;

    zc       = img.fwd_model.electrode( find(I) ).z_contact;

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

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

Image reconstruction

% Solve resistor model
% $Id: tutorial020c.m 3127 2012-06-08 16:19:25Z bgrychtol $

% Now we complete the fwd_model
r_mdl.jacobian= @jacobian_perturb;

% Now create an inverse model
i_mdl= eidors_obj('inv_model','resistor inverse');
i_mdl.fwd_model= r_mdl;
i_mdl.jacobian_bkgnd.value= 1000;

% regulatization not needed for this problem
i_mdl.RtR_prior= @prior_tikhonov;
i_mdl.hyperparameter.value= 0;

i_mdl.reconst_type= 'difference';
i_mdl.solve= @inv_solve_diff_GN_one_step;

% Reconstruct resistor change
reconst= inv_solve(i_mdl, data_1k0, data_1k2);

>>reconst = 
         name: 'solved by aa_inv_solve'
    elem_data: 200.0000
    inv_model: [1x1 struct]
    fwd_model: [1x1 struct]
         type: 'image'