inv_solve_TV_irls

 INV_SOLVE_TV_IRLS Iteratively Reweighted Least Squares inverse solver
 img= ab_inv_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

 Parameters:
  max_iters =  inv_model.parameters.max_iteration (default 10)
      Max number of iterations before stopping
  min change = inv_model.parameters.min_change   (default 0)
      Min Change in objective fcn (norm(y-Jx)^2 + hp*TV(x)) before stopping
  beta      =  inv_model.inv_solve_TV_irls.beta   (default 1e-3)
 beta is the parameter that smooths the TV functional

0001 function img=inv_solve_TV_irls( inv_model, data1, data2)
0002 % INV_SOLVE_TV_IRLS Iteratively Reweighted Least Squares inverse solver
0003 % img= ab_inv_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 % Parameters:
0010 %  max_iters =  inv_model.parameters.max_iteration (default 10)
0011 %      Max number of iterations before stopping
0012 %  min change = inv_model.parameters.min_change   (default 0)
0013 %      Min Change in objective fcn (norm(y-Jx)^2 + hp*TV(x)) before stopping
0014 %  beta      =  inv_model.inv_solve_TV_irls.beta   (default 1e-3)
0015 % beta is the parameter that smooths the TV functional
0016 
0017 % (C) 2008 Andrea Borsic. License: GPL version 2 or version 3
0018 % $Id: inv_solve_TV_irls.m 4026 2013-05-18 09:22:32Z bgrychtol $
0019 
0020 try    max_iter = inv_model.parameters.max_iterations;
0021 catch  max_iter = 10;
0022 end
0023 
0024 try    min_change = inv_model.parameters.min_change;
0025 catch  min_change = 0;
0026 end
0027 
0028 try    beta = inv_model.inv_solve_TV_irls.beta; 
0029 catch  beta = 1e-3;
0030 end
0031 
0032 
0033 fwd_model= inv_model.fwd_model;
0034 d=calc_difference_data( data1, data2, inv_model.fwd_model);
0035 
0036 L=calc_R_prior( inv_model );
0037 
0038 img_bkgnd=calc_jacobian_bkgnd( inv_model );
0039 J=calc_jacobian( img_bkgnd);
0040 
0041 alpha=calc_hyperparameter( inv_model );
0042 
0043 
0044 delta_sigma = zeros(size(J,2),1); % we start from no initial difference
0045 
0046 Obj_Fcn = inf; %initial value
0047 
0048 for k=1:max_iter
0049  
0050     dv =  J*delta_sigma - d;
0051 
0052 % STOP IF OBJECTIVE FCN IS NOT CHANGING by more than min_change
0053     Obj_Fcnk = norm(dv)^2 + alpha*sum(abs( delta_sigma ));
0054     delta_ObjFcn = abs(Obj_Fcnk/Obj_Fcn - 1);
0055 %   fprintf('%d %g %g %g\n',k, Obj_Fcnk, Obj_Fcn, delta_ObjFcn);
0056     if delta_ObjFcn < min_change; 
0057        eidors_msg('Lagged_diff: Breaking at iteration %d',k,2);
0058        break;
0059     end
0060     Obj_Fcn = Obj_Fcnk;
0061     
0062     E= sqrt((L*delta_sigma).^2+beta);
0063     inv_E= spdiags( 1./E, 0, length(E), length(E));
0064 
0065     phi1=J'*dv+ alpha*L'*inv_E*L*delta_sigma;
0066     phi2=J'*J + alpha*L'*inv_E*L;
0067 
0068     upd=-(phi2)\phi1;
0069     
0070     delta_sigma=delta_sigma+upd;
0071 
0072     
0073 end % for k
0074 
0075 % create a data structure to return
0076 img.name = 'inv_solve_TV_irls';
0077 img.elem_data = delta_sigma;
0078 img.fwd_model = fwd_model;