np_inv_solve

 NP_INV_SOLVE inverse solver for Nick Polydorides EIDORS3D code
 img= np_inv_solve( inv_model, data1, data2)
 img        => output image
 inv_model  => inverse model struct
 data1      => differential data at earlier time
 data2      => differential data at later time
 inv_model.parameters.max_iterations (default 1);
 inv_model.parameters.term_tolerance (default 1e-3);

0001 function img= np_inv_solve( inv_model, data1, data2)
0002 % NP_INV_SOLVE inverse solver for Nick Polydorides EIDORS3D code
0003 % img= np_inv_solve( inv_model, data1, data2)
0004 % img        => output image
0005 % inv_model  => inverse model struct
0006 % data1      => differential data at earlier time
0007 % data2      => differential data at later time
0008 % inv_model.parameters.max_iterations (default 1);
0009 % inv_model.parameters.term_tolerance (default 1e-3);
0010 
0011 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0012 % $Id: np_inv_solve.html 2819 2011-09-07 16:43:11Z aadler $
0013 
0014 [maxiter, tol] = get_parameters(inv_model);
0015   
0016 dv = calc_difference_data( data1, data2, inv_model.fwd_model);
0017 
0018 sol = one_step_inv_matrix(inv_model) * dv;
0019 
0020 if maxiter>1
0021    RtR = calc_RtR_prior( inv_model );
0022    hp= calc_hyperparameter( inv_model );
0023 
0024    for iter=2:maxiter
0025       dv_sim= forward_solve_diff(inv_model, sol);
0026       eidors_msg('iter=%d, norm(err)= %f', iter, norm(dv_sim - dv),3);
0027       if norm(dv_sim - dv)<tol; break; end  % test tolerance
0028 
0029       img_bkgnd= calc_jacobian_bkgnd( inv_model );
0030       img_bkgnd.elem_data = img_bkgnd.elem_data + sol;
0031       J = calc_jacobian( inv_model.fwd_model, img_bkgnd);
0032 
0033       sol_upd= (J'*J +  hp^2*RtR)\(J' * (dv - dv_sim));
0034       sol = sol + sol_upd;
0035    end
0036 end
0037 
0038 % create a data structure to return
0039 img.name= 'solved by np_inv_solve';
0040 img.elem_data = sol;
0041 img.fwd_model= inv_model.fwd_model;
0042 
0043 function one_step_inv = one_step_inv_matrix(inv_model)
0044 % The one_step reconstruction matrix is cached
0045    one_step_inv = eidors_obj('get-cache', inv_model, 'np_2003_one_step_inv');
0046    if ~isempty(one_step_inv)
0047        eidors_msg('np_inv_solve: using cached value', 3);
0048    else
0049        img_bkgnd= calc_jacobian_bkgnd( inv_model );
0050        J = calc_jacobian( inv_model.fwd_model, img_bkgnd);
0051 
0052        RtR = calc_RtR_prior( inv_model );
0053        hp= calc_hyperparameter( inv_model );
0054 
0055        % Calculating a linear inverse solution
0056        one_step_inv= (J'*J +  hp^2*RtR)\J';
0057 
0058        eidors_obj('set-cache', inv_model, 'np_2003_one_step_inv', one_step_inv);
0059        eidors_msg('np_inv_solve: setting cached value', 3);
0060    end
0061 
0062 function dv_sim= forward_solve_diff(inv_model, sol) 
0063    img= calc_jacobian_bkgnd( inv_model );
0064    v_homg = fwd_solve(img);
0065    img.elem_data = img.elem_data + sol;
0066    v_solv = fwd_solve(img);
0067    dv_sim= calc_difference_data( v_homg, v_solv, inv_model.fwd_model);
0068 
0069 function [maxiter, tol] = get_parameters(inv_model);
0070 
0071    try
0072      maxiter= inv_model.parameters.max_iterations;
0073    catch
0074      maxiter= 1;
0075    end
0076 
0077    try
0078      tol = inv_model.parameters.term_tolerance;
0079    catch
0080      tol= 1e-3;
0081    end