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)
   hp    => hyperparameter. Ratio between the last kept and the first SV
            in percent (default: .1)

 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 %   hp    => hyperparameter. Ratio between the last kept and the first SV
0012 %            in percent (default: .1)
0013 %
0014 % SEE ALSO: calc_TSVD_RM, calc_hyperpameter, solve_use_matrix
0015 
0016 % (C) 2011 Bartlomiej Grychtol. Licenced under GPL v2 or v3
0017 % $Id: inv_solve_TSVD.m 5431 2017-04-27 15:29:35Z aadler $
0018 
0019 if ischar(inv_model) && strcmp(inv_model,'UNIT_TEST'), do_unit_test, return, end;
0020 
0021 hp = calc_hyperparameter(inv_model);
0022 inv_model.solve_use_matrix.RM = calc_TSVD_RM(inv_model, hp);
0023 sol = solve_use_matrix(inv_model,data1,data2);
0024 
0025 
0026 function do_unit_test
0027     % get a fwd_model
0028     mdl = mk_common_model('c3cr',16);fmdl = mdl.fwd_model; clear mdl;
0029     % homogenous measurement
0030     img  = mk_image(fmdl,1);
0031     vh = fwd_solve(img);
0032     % inhomogeneous measurement
0033     str = sprintf('(x - %f).^2+(y - %f).^2+(z - %f).^2<%f^2',[-0.3 0.2 0 0.1]);
0034     select_fcn = inline(str,'x','y','z');
0035     e = elem_select(img.fwd_model, select_fcn);
0036     ind = find(e);
0037     img.elem_data(ind) = img.elem_data(ind) - 0.25*e(ind);
0038     vi = fwd_solve(img);
0039 
0040     % build an inverse model
0041     imdl.name= 'Test TSVD model';
0042     imdl.type= 'inv_model';
0043     imdl.solve= @inv_solve_TSVD;
0044     imdl.hyperparameter.value = 1;
0045     imdl.jacobian_bkgnd.value = 1;
0046     imdl.reconst_type= 'difference';
0047     imdl.fwd_model = fmdl;
0048     
0049     % solve
0050 subplot(221);
0051     show_slices(img,[inf,inf,0]);
0052 subplot(222);
0053     rimg = inv_solve(imdl, vh, vi);
0054     show_slices(rimg,[inf,inf,0]);
0055 
0056 subplot(223);
0057     imdl.hyperparameter.value = .1;
0058     rimg = inv_solve(imdl, vh, vi);
0059     show_slices(rimg,[inf,inf,0]);