EIDORS: Electrical Impedance Tomography and Diffuse Optical Tomography Reconstruction Software

Hosted by

Compare selection of hyperparameters for Total Variation

Simluation object

% TV: Create object $Id: TV_hyperparams01.m 1535 2008-07-26 15:36:27Z aadler $

imb=  mk_common_model('c2c2',16);
fmdl= imb.fwd_model;
sigma= ones(size(fmdl.elems,1),1);
img= eidors_obj('image','','elem_data',sigma,'fwd_model',fmdl);

vh= fwd_solve( img );

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

subplot(221)
show_fem(img);
axis equal; axis off;
print -r100 -dpng TV_hyperparams01a.png;

Figure: Shape that is to be reconstructed (generated on 576 element mesh)

Create TV reconstruction model

% TV: Reconstruction model $Id: TV_hyperparams02.m 1535 2008-07-26 15:36:27Z aadler $

maxit=20;  % max number of iterations
imdl=mk_common_model('b2c0',16);

invtv= eidors_obj('inv_model', 'EIT inverse');
invtv.reconst_type= 'difference';
invtv.jacobian_bkgnd.value= 1;

invtv.fwd_model=                  imdl.fwd_model;
invtv.solve=                      @ab_tv_diff_solve;
invtv.R_prior=                    @ab_calc_tv_prior;
invtv.parameters.term_tolerance=  1e-6;
invtv.parameters.keep_iterations= 1;
invtv.parameters.max_iterations=  maxit;

subplot(221)
show_fem(invtv.fwd_model);
axis equal; axis off;
print -r100 -dpng TV_hyperparams02a.png;

Figure: Reconstruction model (256 elements)

Reconstruction with no noise

name_base= 'tv_hp_00_img-a1=%3.1f-a2=%3.1f.jpg';
% Iterate over params $Id: TV_hyperparams03.m 1535 2008-07-26 15:36:27Z aadler $
name_base = 'tv_hp_00_img-a1=%3.1f-a2=%3.1f.jpg';
alpha1list= [1.5:0.5:4.5];
alpha2list= [3.5:1.0:9.5];
for alpha2= alpha2list
   for alpha1= alpha1list
      invtv.hyperparameter.value =   10^-alpha2;
      invtv.ab_tv_diff_solve.alpha1= 10^-alpha1;
      imgs= inv_solve(invtv,vh,vi);
      imgs.elem_data= imgs.elem_data(:,[1,2,5,maxit]);
      imgs.calc_colours.window_range=.2;
      imgs.calc_colours.clim=1.2;
      out_img= show_slices(imgs);
      imwrite(out_img, colormap, sprintf(name_base,alpha1,alpha2) );
    end
end

To display these results, create an html table with matlab

fid= fopen('TV-params-NSR=0.html','w');
% Generate HTML frame to view $Id: TV_hyperparams04.m 1535 2008-07-26 15:36:27Z aadler $

a=sprintf('%calpha;',38); % alpha
m=sprintf('%cminus;',38); % alpha
s=sprintf('%c',60); % less than 
e=sprintf('%c',62); % greater than
tr= [s,'TR',e]; etr= [s,'/TR',e];
th= [s,'TH',e]; eth= [s,'/TH',e];
td= [s,'TD',e]; etd= [s,'/TD',e];
sub=[s,'SUB',e];esub=[s,'/SUB',e];
sup=[s,'SUP',e];esup=[s,'/SUP',e];


fprintf(fid,[s,'TABLE',e,tr,th]);
for alpha1= alpha1list
   fprintf(fid,[th,a,sub,'1',esub,'=10',sup,m,'%3.1f',esup],alpha1);
end
fprintf(fid,'\n');

for alpha2= alpha2list
   fprintf(fid,['\n',tr,'\n',th,a,sub,'2',esub,'=10', ...
                sup,m,'%3.1f',esup,'\n'],alpha2);
   for alpha1= alpha1list
      fprintf(fid,[td,s,'img src="%s"',e,'\n'], sprintf(name_base,alpha1,alpha2) );
   end
end

fprintf(fid,['\n',s,'/TABLE',e,'\n']);
fclose(fid);

Reconstruction Results (No Noise) (GN param = α₁, TV param= α₂)

	α1=10−1.5	α1=10−2.0	α1=10−2.5	α1=10−3.0	α1=10−3.5	α1=10−4.0	α1=10−4.5
α2=10−3.5
α2=10−4.5
α2=10−5.5
α2=10−6.5
α2=10−7.5
α2=10−8.5
α2=10−9.5

Reconstruction with 20db SNR noise

Add Noise

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

sig= sqrt(norm(vi.meas - vh.meas));
m= size(vi.meas,1);
vi.meas = vi.meas + .01*sig*randn(m,1);

Calculations

name_base= 'tv_hp_20_img-a1=%3.1f-a2=%3.1f.jpg';
% Iterate over params $Id: TV_hyperparams03.m 1535 2008-07-26 15:36:27Z aadler $
name_base = 'tv_hp_00_img-a1=%3.1f-a2=%3.1f.jpg';
alpha1list= [1.5:0.5:4.5];
alpha2list= [3.5:1.0:9.5];
for alpha2= alpha2list
   for alpha1= alpha1list
      invtv.hyperparameter.value =   10^-alpha2;
      invtv.ab_tv_diff_solve.alpha1= 10^-alpha1;
      imgs= inv_solve(invtv,vh,vi);
      imgs.elem_data= imgs.elem_data(:,[1,2,5,maxit]);
      imgs.calc_colours.window_range=.2;
      imgs.calc_colours.clim=1.2;
      out_img= show_slices(imgs);
      imwrite(out_img, colormap, sprintf(name_base,alpha1,alpha2) );
    end
end

To display these results, create an html table with matlab

fid= fopen('TV-params-NSR=.01.html','w');
% Generate HTML frame to view $Id: TV_hyperparams04.m 1535 2008-07-26 15:36:27Z aadler $

a=sprintf('%calpha;',38); % alpha
m=sprintf('%cminus;',38); % alpha
s=sprintf('%c',60); % less than 
e=sprintf('%c',62); % greater than
tr= [s,'TR',e]; etr= [s,'/TR',e];
th= [s,'TH',e]; eth= [s,'/TH',e];
td= [s,'TD',e]; etd= [s,'/TD',e];
sub=[s,'SUB',e];esub=[s,'/SUB',e];
sup=[s,'SUP',e];esup=[s,'/SUP',e];


fprintf(fid,[s,'TABLE',e,tr,th]);
for alpha1= alpha1list
   fprintf(fid,[th,a,sub,'1',esub,'=10',sup,m,'%3.1f',esup],alpha1);
end
fprintf(fid,'\n');

for alpha2= alpha2list
   fprintf(fid,['\n',tr,'\n',th,a,sub,'2',esub,'=10', ...
                sup,m,'%3.1f',esup,'\n'],alpha2);
   for alpha1= alpha1list
      fprintf(fid,[td,s,'img src="%s"',e,'\n'], sprintf(name_base,alpha1,alpha2) );
   end
end

fprintf(fid,['\n',s,'/TABLE',e,'\n']);
fclose(fid);

Reconstruction Results (20dB SNR) (GN param = α₁, TV param= α₂)

	α1=10−1.5	α1=10−2.0	α1=10−2.5	α1=10−3.0	α1=10−3.5	α1=10−4.0	α1=10−4.5
α2=10−3.5
α2=10−4.5
α2=10−5.5
α2=10−6.5
α2=10−7.5
α2=10−8.5
α2=10−9.5