EIDORS: Electrical Impedance Tomography and Diffuse Optical Tomography Reconstruction Software

Hosted by

Compare 3D image reconstructions

The first step is to create a simulation model

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

imb=  mk_common_model('n3r2',16);
e= size(imb.fwd_model.elems,1);
bkgnd= 1;

% Homogenous Data
img= eidors_obj('image','Demo Image');
img.elem_data= bkgnd*ones(e,1);
img.fwd_model= imb.fwd_model;
vh= fwd_solve( img );

% Inhomogenous Data - Load from file 'datacom'
load datacom A B;
img.elem_data(A)= bkgnd*1.2;
img.elem_data(B)= bkgnd*0.8;
clear A B;
vi= fwd_solve( img );

% Add 12dB noise
vi_n= vi; 
vi_n.meas = vi.meas + std(vi.meas - vh.meas)/10^(12/20) ...
                     *randn(size(vi.meas));
sig= sqrt(norm(vi.meas - vh.meas));

subplot(121);
show_fem(img); axis square;

subplot(122);
show_fem(img); axis square;
crop_model([],  inline('y<0','x','y','z'))
view(-51,14);
print -r100 -dpng tutorial130a.png;

Figure: Simulation image for sample data (two different views)

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

% Create Inverse Model
inv3d= eidors_obj('inv_model', 'EIT inverse');
inv3d.reconst_type= 'difference';
inv3d.jacobian_bkgnd.value = 1;
inv3d.fwd_model= imb.fwd_model;
inv3d.fwd_model.np_fwd_solve.perm_sym= '{y}';

% Nick Polydorides' Gauss-Newton Solver
inv3d.hyperparameter.value = 1e-3;
inv3d.solve=       @np_inv_solve;

% Tikhonov prior
inv3d.R_prior=     @tikhonov_image_prior;
imgr(1)= inv_solve( inv3d, vh, vi);
imgn(1)= inv_solve( inv3d, vh, vi_n);

% Nick Polydorides' Prior (Laplace)
inv3d.R_prior=     @np_calc_image_prior;
inv3d.np_calc_image_prior.parameters= [3 1]; %  deg=1, w=1
imgr(2)= inv_solve( inv3d, vh, vi);
imgn(2)= inv_solve( inv3d, vh, vi_n);

% Andrea Borsic's PDIPM TV solver
inv3d.ab_calc_tv_prior.alpha2 = 1e-5;
inv3d.parameters.max_iterations= 20;
inv3d.parameters.term_tolerance= 1e-3;
inv3d.R_prior=     @ab_calc_tv_prior;
inv3d.solve=       @ab_tv_diff_solve;

imgr(3)= inv_solve( inv3d, vh, vi);
imgn(3)= inv_solve( inv3d, vh, vi_n);

% Output image
posn= [inf,inf,2.5,1,1;inf,inf,1.5,1,2;inf,inf,0.5,1,3];
clf;
show_slices(imgr, posn);
print -r100 -dpng tutorial130b.png;
show_slices(imgn, posn);
print -r100 -dpng tutorial130c.png;

Figure: Images reconstructed with data without noise. Slices are shown at heights of (top to bottom): 1) 2.5, 2) 1.5, 3) 0.5. From Left to Right: 1) One step Gauss-Newton reconstruction (Tikhonov prior) 2) One step Gauss-Newton reconstruction (Laplace filter prior) 3): Total Variation reconstruction

Figure: Images reconstructed with data with 15dB SNR. Slices are shown at heights of (top to bottom): 1) 2.5, 2) 1.5, 3) 0.5. From Left to Right: 1) One step Gauss-Newton reconstruction (Tikhonov prior) 2) One step Gauss-Newton reconstruction (Laplace filter prior) 3): Total Variation reconstruction