compare_2d_algs

 Compare different 2D reconstructions
 [imgr, img]= compare_2d_algs(option);

 imgr - reconstructed image (256 elements)
 img  - original image      (576 elements)

 option -> select algorithm
 OPTION   SOLVER               PRIOR             HP
   1   inv_solve_diff_GN_one_step       prior_laplace   1e-3
   2   np_inv_solve       prior_laplace   1e-3
   3   inv_solve_diff_GN_one_step       prior_gaussian_HPF   NF=2
   3.1 inv_solve_diff_GN_one_step       prior_noser     NF=2
   4   inv_solve_TV_pdipm   prior_TV      1e-4
   5   aa_inv_total_var   prior_laplace   1e-4 (not the usual prior)
   6   aa_inv_total_var   prior_TV      1e-4
   7   aa_inv_conj_grad   prior_TV      ??? 
   8   inv_solve_TV_irls  prior_TV      ???

  OPTION = 1dd => do OPTION=dd with normalize_measurements

 shape
   0  two triangles
   1  round

0001 function [imgr, img]= compare_2d_algs(option,shape);
0002 % Compare different 2D reconstructions
0003 % [imgr, img]= compare_2d_algs(option);
0004 %
0005 % imgr - reconstructed image (256 elements)
0006 % img  - original image      (576 elements)
0007 %
0008 % option -> select algorithm
0009 % OPTION   SOLVER               PRIOR             HP
0010 %   1   inv_solve_diff_GN_one_step       prior_laplace   1e-3
0011 %   2   np_inv_solve       prior_laplace   1e-3
0012 %   3   inv_solve_diff_GN_one_step       prior_gaussian_HPF   NF=2
0013 %   3.1 inv_solve_diff_GN_one_step       prior_noser     NF=2
0014 %   4   inv_solve_TV_pdipm   prior_TV      1e-4
0015 %   5   aa_inv_total_var   prior_laplace   1e-4 (not the usual prior)
0016 %   6   aa_inv_total_var   prior_TV      1e-4
0017 %   7   aa_inv_conj_grad   prior_TV      ???
0018 %   8   inv_solve_TV_irls  prior_TV      ???
0019 %
0020 %  OPTION = 1dd => do OPTION=dd with normalize_measurements
0021 %
0022 % shape
0023 %   0  two triangles
0024 %   1  round
0025 
0026 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0027 % $Id: compare_2d_algs.m 3499 2012-07-04 21:00:30Z bgrychtol $
0028 
0029 if nargin<2
0030     shape=0;
0031 end
0032 
0033 calc_colours('ref_level','auto');
0034 
0035 imb=  mk_common_model('c2c',16);
0036 e= size(imb.fwd_model.elems,1);
0037 sigma= ones(e,1);
0038 img= eidors_obj('image','');
0039 img.elem_data= sigma;
0040 img.fwd_model= imb.fwd_model;
0041 vh= fwd_solve( img );
0042 
0043 if shape==0
0044     sigma([25,37,49:50,65:66,81:83,101:103,121:124])=2;
0045     sigma([95,98:100,79,80,76,63,64,60,48,45,36,33,22])=2;
0046 elseif shape==1
0047     sigma([65,81,82,101,102,122])=2; 
0048 elseif shape==2
0049     sigma(1:4)=2;
0050 else
0051     error('shape not defined (%d)',shape);
0052 end
0053     
0054 img.elem_data= sigma;
0055 vi= fwd_solve( img );
0056 
0057 sig= sqrt(norm(vi.meas - vh.meas));
0058 m= size(vi.meas,1);
0059 vi.meas = vi.meas + .0001*sig*randn(m,1);
0060 figure(2); img.elem_data= img.elem_data - 1; show_slices(img); figure(1);
0061 
0062 %show_slices(img);
0063 imb=  mk_common_model('b2c2',16);
0064 inv2d= eidors_obj('inv_model', 'EIT inverse');
0065 inv2d.reconst_type= 'difference';
0066 inv2d.jacobian_bkgnd.value= 1;
0067 inv2d.fwd_model= imb.fwd_model;
0068 inv2d.fwd_model.np_fwd_solve.perm_sym= '{y}';
0069 inv2d.parameters.term_tolerance= 1e-4;
0070 
0071 if option>100
0072    inv2d.fwd_model = mdl_normalize(inv2d.fwd_model,1);
0073    option=option-100;
0074 end
0075 
0076 switch option
0077    case 1,
0078      inv2d.hyperparameter.value = 1e-3;
0079      inv2d.solve=       'inv_solve_diff_GN_one_step';
0080      inv2d.RtR_prior=   'prior_laplace';
0081 
0082    case 2,
0083      inv2d.hyperparameter.value = 1e-3;
0084      inv2d.RtR_prior=   'prior_laplace';
0085      inv2d.solve=       'np_inv_solve';
0086 
0087    case 3,
0088      inv2d.hyperparameter.func = @choose_noise_figure;
0089      inv2d.hyperparameter.noise_figure= 2;
0090      inv2d.hyperparameter.tgt_elems= 1:4;
0091      inv2d.RtR_prior=   'prior_gaussian_HPF';
0092      inv2d.solve=       'inv_solve_diff_GN_one_step';
0093 
0094    case 3.1,
0095      inv2d.hyperparameter.func = @choose_noise_figure;
0096      inv2d.hyperparameter.noise_figure= 2;
0097      inv2d.hyperparameter.tgt_elems= 1:4;
0098      inv2d.RtR_prior=   @prior_laplace;
0099      inv2d.solve=       'inv_solve_diff_GN_one_step';
0100 
0101    case 4,
0102      inv2d.hyperparameter.value = 1e-6;
0103      inv2d.parameters.max_iterations= 10;
0104      inv2d.R_prior=     'prior_TV';
0105      inv2d.solve=       'inv_solve_TV_pdipm';
0106      inv2d.parameters.keep_iterations=1;
0107 
0108    case 5,
0109      inv2d.hyperparameter.value = 1e-4;
0110      inv2d.solve=       'aa_inv_total_var';
0111      inv2d.R_prior=     'prior_laplace';
0112      inv2d.parameters.max_iterations= 10;
0113 
0114    case 6,
0115      subplot(141); show_slices(img);
0116      inv2d.hyperparameter.value = 1e-4;
0117      inv2d.solve=       'aa_inv_total_var';
0118      inv2d.R_prior=     'prior_TV';
0119      inv2d.parameters.max_iterations= 1;
0120      inv2d.parameters.keep_iterations=1;
0121      subplot(142); show_slices( inv_solve( inv2d, vh, vi) );
0122      inv2d.parameters.max_iterations= 2;                
0123      subplot(143); show_slices( inv_solve( inv2d, vh, vi) );
0124      inv2d.parameters.max_iterations= 5;                
0125      subplot(144); show_slices( inv_solve( inv2d, vh, vi) );
0126      return;
0127 
0128    case 7,
0129      inv2d.hyperparameter.value = 1e-2;
0130      inv2d.parameters.max_iterations = 1e3;
0131      inv2d.parameters.term_tolerance = 1e-3;
0132      inv2d.solve=          'aa_inv_conj_grad';
0133      inv2d.R_prior=        'prior_TV';
0134 
0135      
0136    case 8,
0137      inv2d.hyperparameter.value = 1e-5;
0138      inv2d.parameters.max_iterations= 20;
0139      inv2d.R_prior=     'prior_TV';
0140      inv2d.solve=       'inv_solve_TV_irls';
0141      inv2d.parameters.keep_iterations=1;
0142      
0143    otherwise,
0144      error('action unknown');
0145 end
0146 
0147 %
0148 % Step 3: Reconst and show image
0149 %
0150 imgr= inv_solve( inv2d, vh, vi);
0151 
0152 figure(1); show_slices(imgr);