inv_solve_TSVD

 INV_SOLVE_TSVD: inverse solver based on truncatated SVD
   img= inv_solve_TSVD(inv_model, data1, data2)
 
   img        => output image (or vector of images)
   data1      => differential data at earlier time
   data2      => differential data at later time
   inv_model  => inverse model struct. Requires:
                 inv_model.hyperparameter.value OR
                 inv_model.hyperparameter.func (and possibly others)

 SEE ALSO: calc_TSVD_RM, calc_hyperpameter, solve_use_matrix

0001 function sol = inv_solve_TSVD(inv_model, data1, data2)
0002 % INV_SOLVE_TSVD: inverse solver based on truncatated SVD
0003 %   img= inv_solve_TSVD(inv_model, data1, data2)
0004 %
0005 %   img        => output image (or vector of images)
0006 %   data1      => differential data at earlier time
0007 %   data2      => differential data at later time
0008 %   inv_model  => inverse model struct. Requires:
0009 %                 inv_model.hyperparameter.value OR
0010 %                 inv_model.hyperparameter.func (and possibly others)
0011 %
0012 % SEE ALSO: calc_TSVD_RM, calc_hyperpameter, solve_use_matrix
0013 
0014 % (C) 2011 Bartlomiej Grychtol. Licenced under GPL v2 or v3
0015 % $Id: inv_solve_TSVD.html 2819 2011-09-07 16:43:11Z aadler $
0016 
0017 if isstr(inv_model) && strcmp(inv_model,'UNIT_TEST'), do_unit_test, return, end;
0018 
0019 hp = calc_hyperparameter(inv_model);
0020 inv_model.solve_use_matrix.RM = calc_TSVD_RM(inv_model, hp);
0021 sol = solve_use_matrix(inv_model,data1,data2);
0022 
0023 
0024 function do_unit_test
0025     % get a fwd_model
0026     mdl = mk_common_model('c3cr',16);fmdl = mdl.fwd_model; clear mdl;
0027     % homogenous measurement
0028     img  = mk_image(fmdl,1);
0029     vh = fwd_solve(img);
0030     % inhomogeneous measurement
0031     str = sprintf('(x-%f).^2+(y-%f).^2+(z-%f).^2<%f^2',[-0.3 0.2 0 0.1]);
0032     select_fcn = inline(str,'x','y','z');
0033     e = elem_select(img.fwd_model, select_fcn);
0034     ind = find(e);
0035     img.elem_data(ind) = img.elem_data(ind) - 0.25*e(ind);
0036     vi = fwd_solve(img);
0037 
0038     % build an inverse model
0039     imdl.name= 'Test TSVD model';
0040     imdl.type= 'inv_model';
0041     imdl.solve= @inv_solve_TSVD;
0042     imdl.hyperparameter.value = 1;
0043     imdl.jacobian_bkgnd.value = 1;
0044     imdl.reconst_type= 'difference';
0045     imdl.fwd_model = fmdl;
0046     
0047     % solve
0048     rimg = inv_solve(imdl, vh, vi);
0049     show_slices(rimg,[inf,inf,0]);