nodal_solve

 AA_INV_SOLVE inverse solver using approach of Adler&Guardo 1996
 img= nodal_solve( inv_model, data1, data2)
 img        => output image (or vector of images)
 inv_model  => inverse model struct
 data1      => differential data at earlier time
 data2      => differential data at later time

 both data1 and data2 may be matrices (MxT) each of
  M measurements at T times
 if either data1 or data2 is a vector, then it is expanded
  to be the same size matrix

0001 function img= nodal_solve( inv_model, data1, data2)
0002 % AA_INV_SOLVE inverse solver using approach of Adler&Guardo 1996
0003 % img= nodal_solve( inv_model, data1, data2)
0004 % img        => output image (or vector of images)
0005 % inv_model  => inverse model struct
0006 % data1      => differential data at earlier time
0007 % data2      => differential data at later time
0008 %
0009 % both data1 and data2 may be matrices (MxT) each of
0010 %  M measurements at T times
0011 % if either data1 or data2 is a vector, then it is expanded
0012 %  to be the same size matrix
0013 
0014 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0015 % $Id: nodal_solve.html 2819 2011-09-07 16:43:11Z aadler $
0016 
0017 
0018 dv = calc_difference_data( data1, data2, inv_model.fwd_model);
0019 
0020 sol = get_RM( inv_model ) * dv;
0021 
0022 % create a data structure to return
0023 img.name= 'solved by nodal_solve';
0024 img.node_data = sol;
0025 img.fwd_model= inv_model.fwd_model;
0026 
0027 function one_step_inv = get_RM( inv_model );
0028    fwd_model= inv_model.fwd_model;
0029 
0030    % The one_step reconstruction matrix is cached
0031    one_step_inv = eidors_obj('get-cache', inv_model, 'nodal_solve');
0032    if ~isempty(one_step_inv)
0033        eidors_msg('nodal_solve: using cached value', 3);
0034        return;
0035    end
0036 
0037    img_bkgnd= calc_jacobian_bkgnd( inv_model );
0038    J = calc_jacobian( fwd_model, img_bkgnd);
0039 
0040    RtR = calc_RtR_prior( inv_model );
0041    W   = calc_meas_icov( inv_model );
0042    hp  = calc_hyperparameter( inv_model );
0043 
0044    [e2n,Ne, Nn] = elem2node( fwd_model.elems );
0045    if size(J,2) == Ne;
0046       J= J*e2n;
0047    end
0048    if size(RtR,2) == Ne;
0049       RtR = e2n'*RtR*e2n;
0050       RtR = RtR +  1e-4*spdiags(diag(RtR),0,Nn,Nn); %Need just a little regularization
0051    end
0052 
0053    one_step_inv= (J'*W*J +  hp^2*RtR)\J'*W;
0054 
0055    eidors_obj('set-cache', inv_model, 'nodal_solve', one_step_inv);
0056    eidors_msg('nodal_solve: setting cached value', 3);
0057 
0058 function [e2n, Ne, Nn] = elem2node( elems )
0059    [Ne,d] = size(elems);
0060    elemo = (1:Ne)'*ones(1,d);
0061    e2n = sparse(elemo,elems,1/d); 
0062    Nn = size(e2n,2);