aa_inv_total_var

 AA_INV_TOTAL_VARIANCE inverse solver for difference EIT
 img= aa_inv_total_var( inv_model, data1, data2)

 img        => output image
 inv_model  => inverse model struct
 data1      => differential data at earlier time
 data2      => differential data at later time

 ALGORITHM: for image at iteration k x{k}
 Iterate:
   x{k+1} = (H'*W*H + l^2*R'*WB( x{k} )*R )\H'*W*y
 For:
   WB( x ) = diag( 0.5 ./ sqrt( (l*R*x).^2 + Beta ) )
 With Beta > 0 to approximate the L1 norm

 Comment on hyperparameter, we have l^2 in the iteration
 formula and l in WB, so overall hyperparameter is l. This
 makes the units work for the L1 norm.

0001 function img= aa_inv_total_var( inv_model, data1, data2)
0002 % AA_INV_TOTAL_VARIANCE inverse solver for difference EIT
0003 % img= aa_inv_total_var( inv_model, data1, data2)
0004 %
0005 % img        => output image
0006 % inv_model  => inverse model struct
0007 % data1      => differential data at earlier time
0008 % data2      => differential data at later time
0009 %
0010 % ALGORITHM: for image at iteration k x{k}
0011 % Iterate:
0012 %   x{k+1} = (H'*W*H + l^2*R'*WB( x{k} )*R )\H'*W*y
0013 % For:
0014 %   WB( x ) = diag( 0.5 ./ sqrt( (l*R*x).^2 + Beta ) )
0015 % With Beta > 0 to approximate the L1 norm
0016 %
0017 % Comment on hyperparameter, we have l^2 in the iteration
0018 % formula and l in WB, so overall hyperparameter is l. This
0019 % makes the units work for the L1 norm.
0020  
0021 % REFS:
0022 % Andrea Borsic, Ph.D. Thesis, Oxford Brookes, UK
0023 % See Section 7.2.4 Lagged Diffusivity
0024 % http://www.sc-aip.com/thesis/Andrea%20Borsic%20PhD.pdf
0025 %
0026 % W. Clem Karl, "Regularization in Image Restoration
0027 % and Reconstruction," in Handbook of Image and Video
0028 % Processing, A. Bovic, Ed. San Diego, CA: Academic Press,
0029 % 2000, ch. 3.6, pp. 141-160.
0030 
0031 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0032 % $Id: aa_inv_total_var.html 2819 2011-09-07 16:43:11Z aadler $
0033 
0034 fwd_model= inv_model.fwd_model;
0035 pp= aa_fwd_parameters( fwd_model );
0036 
0037 img_bkgnd= calc_jacobian_bkgnd( inv_model );
0038 J = calc_jacobian( fwd_model, img_bkgnd);
0039 
0040 R = calc_R_prior( inv_model );
0041 W = calc_meas_icov( inv_model );
0042 hp= calc_hyperparameter( inv_model );
0043 
0044 
0045 l_data1= length(data1); l1_0 = l_data1 ~=0;
0046 l_data2= length(data2); l2_0 = l_data2 ~=0;
0047 l_data= max( l_data1, l_data2 );
0048 
0049 dva = calc_difference_data( data1, data2, inv_model.fwd_model);
0050 
0051 imax= 3 ; % we need to do this avoid matlab version issues
0052 if isfield(inv_model, 'parameters')
0053    if isfield(inv_model.parameters, 'max_iterations')
0054       imax= inv_model.parameters.max_iterations;
0055    end
0056 end
0057 etol= 1e-3;
0058 n_img= size(dva,2);
0059 
0060 sol = zeros( size(J,2), n_img );
0061 for i=1:n_img
0062    sol(:,i) = tv_inv( J, W, hp*R, dva(:,i), imax, etol );
0063 end
0064 
0065 % create a data structure to return
0066 img.name= 'solved by aa_inv_conj_grad';
0067 img.elem_data = sol;
0068 img.inv_model= inv_model;
0069 img.fwd_model= fwd_model;
0070 
0071 function x= tv_inv( H, W, R, y, imax, etol )
0072    n= size(R,1);
0073    x= (H'*W*H + R'*R)\H'*W*y; %take reasonable first step
0074    Beta = 1e-4;
0075    for iter= 2:imax
0076       eidors_msg('AA_INV_TV: iteration %d',iter,3);
0077       WB = spdiags( 0.5 ./ sqrt( (R*x).^2 + Beta ), 0, n, n );
0078       x = (H'*W*H + R'*WB*R )\H'*W*y;
0079    end