GN_abs_solve

 GN_ABS_SOLVER absolute solver using Gauss Newton approximation
 img= gn_abs_solve( inv_model, data1, data2)
 img        => output image (or vector of images)
 inv_model  => inverse model struct
 data1      => EIT data

 Parameters:
   inv_model.parameters.max_iterations = N_max iter

0001 function img= GN_abs_solve( inv_model, data1);
0002 % GN_ABS_SOLVER absolute solver using Gauss Newton approximation
0003 % img= gn_abs_solve( inv_model, data1, data2)
0004 % img        => output image (or vector of images)
0005 % inv_model  => inverse model struct
0006 % data1      => EIT data
0007 %
0008 % Parameters:
0009 %   inv_model.parameters.max_iterations = N_max iter
0010 
0011 % (C) 2010 Andy Adler. License: GPL version 2 or version 3
0012 % $Id: GN_abs_solve.html 2819 2011-09-07 16:43:11Z aadler $
0013 
0014 % Step 1: fit to background
0015 img = homogeneous_estimate( inv_model, data1 );
0016 
0017 hp  = calc_hyperparameter( inv_model );
0018 RtR = calc_RtR_prior( inv_model );
0019 W   = calc_meas_icov( inv_model );
0020 hp2RtR= hp*RtR;
0021 
0022 iters = 1;
0023 try 
0024    iters = inv_model.parameters.max_iterations; 
0025 end
0026 
0027 img0 = img;
0028 for i = 1:iters  
0029   vsim = fwd_solve( img ); 
0030   dv = calc_difference_data( vsim , data1, img.fwd_model);
0031   J = calc_jacobian( img );
0032 
0033   RDx = hp2RtR*(img0.elem_data - img.elem_data);
0034   dx = (J'*W*J + hp2RtR)\(J'*dv + RDx);
0035 
0036   img = line_optimize(img, dx, data1);
0037 end
0038 
0039 % Fit a parabola to the linefit and pick the best point
0040 % This is better than doing an exhaustive search
0041 function  img = line_optimize(imgk, dx, data1);
0042   flist = [ 0.1,  0.5, 1.0];
0043   clim = mean(imgk.elem_data)/10; % prevent zero and negative conductivity
0044   img = imgk;
0045   for i = 1:length(flist);
0046      img.elem_data = imgk.elem_data + flist(i)*dx;
0047      img.elem_data(img.elem_data <= clim ) = clim;
0048      vsim = fwd_solve( img );
0049      dv = calc_difference_data( vsim , data1, img.fwd_model);
0050      mlist(i) = norm(dv);
0051   end
0052   pf = polyfit(flist, mlist, 2);
0053   fmin = -pf(2)/pf(1)/2; % poly minimum
0054   fmin(fmin>1) = 1; fmin(fmin<0) = 0;
0055 
0056   img.elem_data = imgk.elem_data + flist(i)*dx;
0057   img.elem_data(img.elem_data <= clim ) = clim;
0058 
0059 function img = homogeneous_estimate( imdl, data );
0060    img = calc_jacobian_bkgnd( imdl );
0061    vs = fwd_solve(img);
0062 
0063    if isstruct(data)
0064       data = data.meas;
0065    else
0066      meas_select = [];
0067      try
0068         meas_select = imdl.fwd_model.meas_select;
0069      end
0070      if length(data) == length(meas_select)
0071         data = data(meas_select);
0072      end
0073    end
0074 
0075    pf = polyfit(data,vs.meas,1);
0076 
0077    if isfield(img.fwd_model,'coarse2fine');
0078 % TODO: the whole coarse2fine needs work here.
0079 %   what happens if c2f doesn't cover the whole region
0080       nc = size(img.fwd_model.coarse2fine,2);
0081       img.elem_data = mean(img.elem_data)*ones(nc,1)*pf(1);
0082    else
0083       img.elem_data = img.elem_data*pf(1);
0084    end