EIDORS: Electrical Impedance Tomography and Diffuse Optical Tomography Reconstruction Software

EIDORS (mirror)
− Image Reconst
− Data Structures
− Applications
− FEM Modelling
− Old tutorials
Contrib Data
Browse Docs
Browse SVN

Mailing list


Hosted by
SourceForge.net Logo


GN reconstruction vs Backprojection

The most commonly used EIT reconstruction algorithm is the backprojection algorithm of Barber and Brown, as implemented in the Sheffield Mk I EIT system.

This tutorial explores the concept of reconstruction from electrical field projections, but does not implement the original backprojection algorithm.

The reconstruction matrix output by the original algorithm has been made available as part of EIDORS, here

Calculate the nodal voltage field

%$Id: backproj_solve01.m 4067 2013-05-26 20:17:46Z bgrychtol $
imdl= mk_common_model('c2c',16);

fmdl= imdl.fwd_model;
    ee= fmdl.elems';
    xx= reshape( fmdl.nodes(ee,1),3, []); 
    yy= reshape( fmdl.nodes(ee,2),3, []); 

for i=1:4
   if     i==1; stim= [0,1];
   elseif i==2; stim= [0,2];
   elseif i==3; stim= [0,4];
   elseif i==4; stim= [0,8];

   fmdl.stimulation = mk_stim_patterns(16,1,stim,[0 1], {}, 1);
   img= mk_image(fmdl, 1);
   node_v= calc_all_node_voltages( img );
   zz{i}= reshape(     node_v(ee,4),3, []); 

for i=1:4
   subplot(2,4,i); cla
   patch(xx,yy,zz{i},zz{i}); view(0, 4); axis off

print_convert backproj_solve01a.png

for i=1:4
   subplot(2,4,i); cla
   patch(xx,yy,zz{i},zz{i}); view(0,34); axis off
print_convert backproj_solve01b.png

Figure: Nodal voltages in a mesh with different stimulation patterns. From Left to right: Adjacent stimulation ([0 1]), 45° stimulation ([0 2]), 90° stimulation ([0 4]), 180° stimulation ([0 8])

Calculate Equipotential lines

% $Id: backproj_solve02.m 3879 2013-04-18 09:26:22Z aadler $

clf; img = {};
% Solve voltage for 3 different models
for idx=1:3
  if     idx==1; mdltype= 'd2C2';
  elseif idx==2; mdltype= 'd2d3c';
  elseif idx==3; mdltype= 'd2T2';

  pat = 4; % Stimulation pattern to show

  imdl= mk_common_model(mdltype,16);
  img{idx} = mk_image(imdl); 
  stim = mk_stim_patterns(16,1,'{ad}','{mono}',{'meas_current','rotate_meas'},-1);
  img{idx}.fwd_model.stimulation = stim(pat);
  img{idx}.fwd_solve.get_all_meas = 1;

% Show raw voltage pattern
for idx=1:3
  vh = fwd_solve(img{idx});
  imgn = rmfield(img{idx},'elem_data');
  imgn.node_data= vh.volt;
  subplot(2,3,idx); show_fem(imgn);
print_convert backproj_solve02a.png

% Calculate Equipotential lines
for idx=1:3
  vh = fwd_solve(img{idx});
  imgn = rmfield(img{idx},'elem_data');

  imgn.node_data= zeros(size(vh.volt,1),1);
  for v = 2:16
     imgn.node_data(vh.volt > vh.meas(v)) = v;

  subplot(2,3,idx); show_fem(imgn);
print_convert backproj_solve02b.png

Figure: Equipotential lines for adjacent stimulation From Left to right: 1024 element circular mesh, 6746 element electrode refined element circular mesh, 1024 element mesh of human upper thorax

Create 2D Model and simulate measurements

Here, we reuse the model from the this tutorial, below:
% Compare 2D algorithms
% $Id: tutorial120a.m 3273 2012-06-30 18:00:35Z aadler $

imb=  mk_common_model('c2c',16);

e= size(imb.fwd_model.elems,1);
bkgnd= 1;

% Solve Homogeneous model
img= mk_image(imb.fwd_model, bkgnd);
vh= fwd_solve( img );

% Add Two triangular elements
img.elem_data([25,37,49:50,65:66,81:83,101:103,121:124])=bkgnd * 2;
img.elem_data([95,98:100,79,80,76,63,64,60,48,45,36,33,22])=bkgnd * 2;
vi= fwd_solve( img );

% Add -12dB SNR
vi_n= vi; 
nampl= std(vi.meas - vh.meas)*10^(-18/20);
vi_n.meas = vi.meas + nampl *randn(size(vi.meas));

axis square; axis off
print_convert('tutorial120a.png', '-density 60')

Figure: Sample image to test 2D image reconstruction algorithms

Image reconstruction with GN (NOSER) and Backprojection solvers

% $Id: backproj_solve03.m 4839 2015-03-30 07:44:50Z aadler $

tutorial120a; % get the model from a related tutorial

% Gauss Newton Solver
inv_GN= eidors_obj('inv_model','GN_solver','fwd_model', img.fwd_model);
inv_GN.reconst_type= 'difference';
inv_GN.solve= @inv_solve_diff_GN_one_step;
inv_GN.RtR_prior= @prior_noser;
inv_GN.jacobian_bkgnd.value= 1;
inv_GN.hyperparameter.value= 0.2;

imgr= inv_solve(inv_GN, vh,vi);
subplot(131); show_fem(imgr);
axis equal; axis off

% Backprojection Solver
inv_BP= eidors_obj('inv_model','BP_solver','fwd_model', img.fwd_model);
inv_BP.reconst_type= 'difference';
inv_BP.solve= @inv_solve_backproj;
inv_BP.inv_solve_backproj.type= 'naive';

imgr= inv_solve(inv_BP, vh,vi);
subplot(132); show_fem(imgr);
axis equal; axis off

inv_BP.inv_solve_backproj.type= 'simple_filter';

imgr= inv_solve(inv_BP, vh,vi);
subplot(133); show_fem(imgr);
axis equal; axis off

print_convert inv_solve_backproj03a.png

Figure: Reconstructions using: Left: GN Solver with NOSER Prior Middle: Naive Backprojection Right: Backprojection with a simple filter

GN vs Sheffield Backprojection

There are several different versions of the backprojection algorithm in existence. The matrix available with EIDORS and shown here is the version distributed with the Sheffield Mk I system, and is very similar to the algorithm distributed with the Göttingen Goe MF II EIT system. Almost all clinical and experimental publications which mention "backprojection" use this version of the algorithm.

% Sheffield MKI backprojection $Id: backproj_solve04.m 4839 2015-03-30 07:44:50Z aadler $

% Gauss Newton Solvers
inv_mdl{1} = inv_GN;
inv_mdl{1}.hyperparameter.value= 0.2;
inv_mdl{2} = inv_GN;
inv_mdl{2}.hyperparameter.value= 0.5;
% Sheffield Backprojection solver
inv_mdl{3} = mk_common_gridmdl('backproj');

for loop=1:3
   imgr= inv_solve(inv_mdl{loop}, vh,vi);
   subplot(1,3,loop); show_fem(imgr);
   axis equal; axis off

print_convert inv_solve_backproj04a.png

Figure: Reconstructions using: Left: GN Solver with NOSER Prior (small hyperparameter) Middle:GN Solver with NOSER Prior (larger hyperparameter) Right: Sheffield Mk I Backprojection

Last Modified: $Date: 2017-02-28 13:12:08 -0500 (Tue, 28 Feb 2017) $ by $Author: aadler $