Eidors-logo    

EIDORS: Electrical Impedance Tomography and Diffuse Optical Tomography Reconstruction Software

EIDORS (mirror)
Main
Documentation
Tutorials
− Image Reconst
− Data Structures
− Applications
− FEM Modelling
− GREIT
− Old tutorials
Workshop
Download
Contrib Data
GREIT
Browse Docs
Browse SVN

News
Mailing list
(archive)
FAQ
Developer
                       

 

Hosted by
SourceForge.net Logo

 

Absolute Imaging: Non-Linear Gauss Newton Image Reconstruction

This tutorial shows use of EIDORS for static image reconstructions. First, use a large (20951 element) model to simulate target inhomogeneities. The mat file with the FEM models may be downloaded here (tutorial151_model.mat).
% Create model for Absolute Imaging
% Based on example from Stephen Murphy
% $Id: tutorial151a.m 2784 2011-07-14 21:25:40Z bgrychtol $

load tutorial151_model.mat

% Create Simulation Model
sim_img= eidors_obj('image', 'Simulation Image');
sim_img.fwd_model= simmdl;

backgnd= .02; 
sim_img.elem_data= backgnd*ones(size(simmdl.elems,1),1);
sim_img.elem_data(target1)= backgnd*2;
sim_img.elem_data(target2)= backgnd/4;

clf;
axes('position',[.1,.1,.65,.8]);
   show_fem(sim_img); view(-33,20);
axes('position',[.8,.1,.15,.8]);
   show_slices(sim_img,5);

print_convert tutorial151a.png '-density 75'

Figure: Left: Simulation Model (20951 elements), showing a non-conductive (blue) and conductive (red) target Right: Simulation Model viewed as equispaced horizontal slices
Simulation measured voltages using this model. We simulate both homogeneous measurements and target measurements in order to do both difference and absolute imaging.
% Simulate data for difference and absolute imaging
% $Id: tutorial151b.m 2784 2011-07-14 21:25:40Z bgrychtol $

% Create homogeneous simulation Model
backgnd= .02; 
sim_img.elem_data= backgnd*ones(size(simmdl.elems,1),1);
v_homg= fwd_solve(sim_img);

% Create target simulation Model
sim_img.elem_data(target1)= backgnd*2;
sim_img.elem_data(target2)= backgnd/4;
v_targ= fwd_solve(sim_img);

clf; subplot(211);
plot([v_homg.meas, v_targ.meas]);
print_convert tutorial151b.png '-density 75'

Figure: Green Homogeneous measurements, Blue Target measurements,
Using a smaller model for the inverse, we can reconstruct this pattern with difference imaging.
% Difference imaging result
% $Id: tutorial151c.m 4839 2015-03-30 07:44:50Z aadler $

% imdl is loaded from file tutorial151_model.mat
imdl.reconst_type= 'difference';
imdl.solve=        @inv_solve_diff_GN_one_step;
imdl.jacobian_bkgnd.value= backgnd;
imdl.hyperparameter.value= .1;
img_diff= inv_solve(imdl, v_homg, v_targ);

clf;
axes('position',[.1,.1,.65,.8]);
   show_fem(img_diff); view(-33,20);
axes('position',[.8,.1,.15,.8]);
   show_slices(img_diff,5);

print_convert tutorial151c.png '-density 75'

Figure: Left: Reconstructed Image (2266 elements), showing a non-conductive (blue) and conductive (red) target Right: Simulation Model viewed as equispaced horizontal slices
Absolute Imaging: First, the download the following function file tutorial151_nonlinearGN.m and save it to a directory in your path, or copy-and-paste the file here:
function img= tutorial151_nonlinearGN( inv_model, data )
% TUTORIAL151_NONLINEARGN Non-Linear EIT Inverse
% img        => output image (or vector of images)
% inv_model  => inverse model struct
% data       => measurement data
% $Id: tutorial151_nonlinearGN.m 3661 2012-11-20 21:18:01Z bgrychtol $

fwd_model= inv_model.fwd_model;

img= eidors_obj('image','Solved by tutorial151_nonlinearGN');
img.fwd_model= fwd_model;
sol= inv_model.tutorial151_nonlinearGN.init_backgnd * ...
     ones(size(fwd_model.elems,1),1);

RtR = calc_RtR_prior( inv_model );
hp2  = calc_hyperparameter( inv_model )^2;

factor= 0; norm_d_data= inf;
for iter= 1:inv_model.parameters.max_iterations
   img.elem_data= sol;
   simdata= fwd_solve( img );

   d_data= data - simdata.meas;
   prev_norm_d_data= norm_d_data; norm_d_data= norm(d_data); 
   eidors_msg('tutorial151_nonlinearGN: iter=%d err=%f factor=%f', ...
           iter,norm_d_data, factor, 2);

   if prev_norm_d_data - norm_d_data < inv_model.parameters.term_tolerance
      break;
   end
   
   J = calc_jacobian( img);
   delta_sol = (J.'*J + hp2*RtR)\ (J.' * d_data);
   factor= linesearch(img, data, sol, delta_sol, norm_d_data);
   sol = sol + factor*delta_sol;
end

   img.elem_data= sol;

% Test several different factors to see which minimizes best
function factor= linesearch(img, data, sol, delta_sol, norm_d_data);
   facts= [0, logspace(-3,0,19)];  
   norms= norm_d_data;
   for f= 2:length(facts);
      img.elem_data= sol + facts(f)*delta_sol;
      simdata= fwd_solve( img );
      norms(f)= norm(data - simdata.meas);
   end
   ff= find(norms==min(norms));
   factor= facts(ff(end));
Using this function, the following code will do absolute imaging
% Absolute imaging result
% $Id: tutorial151d.m 4839 2015-03-30 07:44:50Z aadler $

% imdl is loaded from file tutorial151_model.mat
imdl.reconst_type= 'static';
imdl.solve=        @tutorial151_nonlinearGN;
imdl.parameters.term_tolerance= 1e-5;
imdl.parameters.max_iterations= 10;
% special parameter for this model
imdl.tutorial151_nonlinearGN.init_backgnd= backgnd;

imdl.hyperparameter.value= 2;
img_diff= inv_solve(imdl, v_targ);

clf;
axes('position',[.1,.1,.65,.8]);
   show_fem(img_diff); view(-33,20);
axes('position',[.8,.1,.15,.8]);
   show_slices(img_diff,5);

print_convert tutorial151d.png '-density 75'
Running this code gives the following output:
>> tutorial151d;
EIDORS:[ inv_solve:nonlinear Gauss Newton for absolute imaging demo ]
EIDORS:[ tutorial151_nonlinearGN: iter=1 err=52.820965 factor=0.000000 ]
EIDORS:[ tutorial151_nonlinearGN: iter=2 err=11.468709 factor=0.014678 ]
EIDORS:[ tutorial151_nonlinearGN: iter=3 err=7.208743 factor=0.014678 ]
EIDORS:[ tutorial151_nonlinearGN: iter=4 err=2.284559 factor=0.031623 ]
EIDORS:[ tutorial151_nonlinearGN: iter=5 err=0.850762 factor=0.014678 ]
EIDORS:[ tutorial151_nonlinearGN: iter=6 err=0.599933 factor=0.021544 ]
EIDORS:[ tutorial151_nonlinearGN: iter=7 err=0.519560 factor=0.031623 ]
EIDORS:[ tutorial151_nonlinearGN: iter=8 err=0.341341 factor=0.068129 ]
EIDORS:[ tutorial151_nonlinearGN: iter=9 err=0.307938 factor=0.014678 ]
EIDORS:[ tutorial151_nonlinearGN: iter=10 err=0.262837 factor=0.068129 ]

Figure: Left: Reconstructed absolute Image (2266 elements), showing a non-conductive (blue) and conductive (red) target Right: Simulation Model viewed as equispaced horizontal slices

Last Modified: $Date: 2017-03-01 08:44:21 -0500 (Wed, 01 Mar 2017) $ by $Author: aadler $