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.m 3961 2013-04-22 09:48:11Z aadler $
0013 
0014 warning('EIDORS:deprecated','NP_INV_SOLVE is deprecated as of 07-Jun-2012. Use INV_SOLVE_DIFF_GN_ONE_STEP instead.');
0015 
0016 [maxiter, tol] = get_parameters(inv_model);
0017   
0018 dv = calc_difference_data( data1, data2, inv_model.fwd_model);
0019 
0020 sol = one_step_inv_matrix(inv_model) * dv;
0021 
0022 if maxiter>1
0023    RtR = calc_RtR_prior( inv_model );
0024    hp= calc_hyperparameter( inv_model );
0025 
0026    for iter=2:maxiter
0027       dv_sim= forward_solve_diff(inv_model, sol);
0028       eidors_msg('iter=%d, norm(err)= %f', iter, norm(dv_sim - dv),3);
0029       if norm(dv_sim - dv)<tol; break; end  % test tolerance
0030 
0031       img_bkgnd= calc_jacobian_bkgnd( inv_model );
0032       img_bkgnd.elem_data = img_bkgnd.elem_data + sol;
0033       J = calc_jacobian(img_bkgnd);
0034 
0035       sol_upd= (J'*J +  hp^2*RtR)\(J' * (dv - dv_sim));
0036       sol = sol + sol_upd;
0037    end
0038 end
0039 
0040 % create a data structure to return
0041 img.name= 'solved by np_inv_solve';
0042 img.elem_data = sol;
0043 img.fwd_model= inv_model.fwd_model;
0044 
0045 function one_step_inv = one_step_inv_matrix(inv_model)
0046 % The one_step reconstruction matrix is cached
0047    one_step_inv = eidors_obj('get-cache', inv_model, 'np_2003_one_step_inv');
0048    if ~isempty(one_step_inv)
0049        eidors_msg('np_inv_solve: using cached value', 3);
0050    else
0051        img_bkgnd= calc_jacobian_bkgnd( inv_model );
0052        J = calc_jacobian( img_bkgnd);
0053 
0054        RtR = calc_RtR_prior( inv_model );
0055        hp= calc_hyperparameter( inv_model );
0056 
0057        % Calculating a linear inverse solution
0058        one_step_inv= (J'*J +  hp^2*RtR)\J';
0059 
0060        eidors_obj('set-cache', inv_model, 'np_2003_one_step_inv', one_step_inv);
0061        eidors_msg('np_inv_solve: setting cached value', 3);
0062    end
0063 
0064 function dv_sim= forward_solve_diff(inv_model, sol) 
0065    img= calc_jacobian_bkgnd( inv_model );
0066    v_homg = fwd_solve(img);
0067    img.elem_data = img.elem_data + sol;
0068    v_solv = fwd_solve(img);
0069    dv_sim= calc_difference_data( v_homg, v_solv, inv_model.fwd_model);
0070 
0071 function [maxiter, tol] = get_parameters(inv_model);
0072 
0073    try
0074      maxiter= inv_model.parameters.max_iterations;
0075    catch
0076      maxiter= 1;
0077    end
0078 
0079    try
0080      tol = inv_model.parameters.term_tolerance;
0081    catch
0082      tol= 1e-3;
0083    end