EIDORS: Electrical Impedance Tomography and Diffuse Optical Tomography Reconstruction Software

Hosted by

GN reconstruction vs Backprojection

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 1535 2008-07-26 15:36:27Z aadler $
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, []); 

homog= ones(size(fmdl.elems,1),1);
img= eidors_obj('image','','elem_data',homog);


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];
   end

   fmdl.stimulation = mk_stim_patterns(16,1,stim,[0 1], {}, 1);

   img.fwd_model= fmdl;

   node_v= calc_all_node_voltages( img );

   % show voltages at measurement 2
   zz= reshape(     node_v(ee,4),3, []); 
   subplot(2,4,i); cla
   patch(xx,yy,zz,zz); view(0, 4); axis off
   subplot(2,4,i+4); cla
   patch(xx,yy,zz,zz); view(0,34); axis off
end

print -r125 -dpng backproj_solve01a.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 1535 2008-07-26 15:36:27Z aadler $

for idx=1:3
  if     idx==1; mdltype= 'b2c';
  elseif idx==2; mdltype= 'd2c2';
  elseif idx==3; mdltype= 'd2t2';
  end

   imdl= mk_common_model(mdltype,16);
   fmdl= imdl.fwd_model;
   params= aa_fwd_parameters( fmdl );

   homog= ones(size(fmdl.elems,1),1);
   img= eidors_obj('image','','elem_data',homog);
   img.fwd_model= fmdl;

   node_v= calc_all_node_voltages( img );
   elem_v= reshape(node_v( fmdl.elems',:),3,[],16);
   elem_v= squeeze(mean(elem_v,1));
   meas_v= params.N2E * node_v;

   sel= 4;
   ed= elem_v(:,sel);
   img.elem_data= ed;
   subplot(2,3,idx);
   show_fem(img);

   ej = zeros(size(ed));
   for i=1:16
     ej= ej + ( ed < meas_v(i,sel) );
   end
   img.elem_data= ej;
   subplot(2,3,idx+3);
   show_fem(img);

end

print -r125 -dpng backproj_solve02a.png;

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

Create 2D Model and simulate measurements

% Compare 2D algorithms
% $Id: tutorial120a.m 1535 2008-07-26 15:36:27Z aadler $

imb=  mk_common_model('c2c',16);

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

% Solve Homogeneous model
img= eidors_obj('image','');
img.elem_data= bkgnd * ones(e,1);
img.fwd_model= imb.fwd_model;
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));

subplot(221)
show_fem(img);
axis square; axis off
print -r100 -dpng tutorial120a.png;

Figure: Sample image to test 2D image reconstruction algorithms

Image reconstruction with GN (NOSER) and Backprojection solvers

% $Id: backproj_solve03.m 1535 2008-07-26 15:36:27Z aadler $

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

imgr= inv_solve(inv_GN, vh,vi);
imgr.calc_colours.ref_level=0;
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= @backproj_solve;
inv_BP.backproj_solve.type= 'naive';

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

inv_BP.backproj_solve.type= 'simple_filter';

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

print -r125 -dpng backproj_solve03a.png;

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

GN vs Sheffield Backprojection

% Sheffield MKI backprojection $Id: backproj_solve04.m 1535 2008-07-26 15:36:27Z aadler $

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

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

print -r125 -dpng backproj_solve04a.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