EIDORS: Electrical Impedance Tomography and Diffuse Optical Tomography Reconstruction Software

Hosted by

Compare 2D image reconstructions

% 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

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

% Create Inverse Model
inv2d= eidors_obj('inv_model', 'EIT inverse');
inv2d.reconst_type= 'difference';
inv2d.jacobian_bkgnd.value= 1;

% This is not an inverse crime; inv_mdl != fwd_mdl
imb=  mk_common_model('b2c',16);
inv2d.fwd_model= imb.fwd_model;
inv2d.fwd_model.np_fwd_solve.perm_sym= '{y}';

% Guass-Newton solvers
inv2d.solve=       @np_inv_solve;

% Tikhonov prior
inv2d.hyperparameter.value = 1e-3;
inv2d.RtR_prior=   @tikhonov_image_prior;
imgr(1)= inv_solve( inv2d, vh, vi);
imgn(1)= inv_solve( inv2d, vh, vi_n);

% NOSER prior
inv2d.hyperparameter.value = 1e-2;
inv2d.RtR_prior=   @noser_image_prior;
imgr(2)= inv_solve( inv2d, vh, vi);
imgn(2)= inv_solve( inv2d, vh, vi_n);

% Laplace image prior
inv2d.hyperparameter.value = 1e-3;
inv2d.RtR_prior=   @laplace_image_prior;
imgr(3)= inv_solve( inv2d, vh, vi);
imgn(3)= inv_solve( inv2d, vh, vi_n);

% Automatic hyperparameter selection
inv2d.hyperparameter = rmfield(inv2d.hyperparameter,'value');
inv2d.hyperparameter.func = @choose_noise_figure;
inv2d.hyperparameter.noise_figure= 0.5;
inv2d.hyperparameter.tgt_elems= 1:4;
inv2d.RtR_prior=   @gaussian_HPF_prior;
inv2d.solve=       @aa_inv_solve;
imgr(4)= inv_solve( inv2d, vh, vi);
imgn(4)= inv_solve( inv2d, vh, vi_n);
inv2d.hyperparameter = rmfield(inv2d.hyperparameter,'func');

% Total variation using PDIPM
inv2d.hyperparameter.value = 1e-5;
inv2d.solve=       @ab_tv_diff_solve;
inv2d.R_prior=     @ab_calc_tv_prior;
inv2d.parameters.max_iterations= 10;
inv2d.parameters.term_tolerance= 1e-3;

imgr(5)= inv_solve( inv2d, vh, vi);
imgn(5)= inv_solve( inv2d, vh, vi_n);

% Output image
show_slices(imgr, [inf,inf,0,1,1]);
print -r125 -dpng tutorial120b.png;
show_slices(imgn, [inf,inf,0,1,1]);
print -r125 -dpng tutorial120c.png;

Figure: Images reconstructed with data without noise, From Left to Right: 1) One step Gauss-Newton reconstruction (Tikhonov prior) 2) One step Gauss-Newton reconstruction (NOSER prior) 3) One step Gauss-Newton reconstruction (Laplace filter prior) 5): One step Gauss-Newton reconstruction (automatic hyperparameter selection) 5): Total Variation reconstruction

Figure: Images reconstructed with data with added 12dB SNR. From Left to Right: 1) One step Gauss-Newton reconstruction (Tikhonov prior) 2) One step Gauss-Newton reconstruction (NOSER prior) 3) One step Gauss-Newton reconstruction (Laplace filter prior) 5): One step Gauss-Newton reconstruction (automatic hyperparameter selection) 5): Total Variation reconstruction