calc_R_prior

PURPOSE ^

R = calc_R_prior( inv_model, varargin )

SYNOPSIS ^

function R_prior = calc_R_prior( inv_model, varargin )

DESCRIPTION ^

 R = calc_R_prior( inv_model, varargin )
 CALC_R_PRIOR: calculate regularization matrix R
   The image prior is matrix n_elem x ??? 
 
 calc_R_prior can be called as
    R_prior= calc_R_prior( inv_model, ... )

 and will call the function inv_model.R_prior
 parameters to R_prior should be passed in the field
 inv_model.R_prior_function_name.parameters
 
 If inv_model.R_prior is a matrix, calc_R_prior will return that matrix,
 possibly correcting for coarse2fine

 R_prior      calculated regularization prior R
 inv_model    is an inv_model structure

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function R_prior = calc_R_prior( inv_model, varargin )
0002 % R = calc_R_prior( inv_model, varargin )
0003 % CALC_R_PRIOR: calculate regularization matrix R
0004 %   The image prior is matrix n_elem x ???
0005 %
0006 % calc_R_prior can be called as
0007 %    R_prior= calc_R_prior( inv_model, ... )
0008 %
0009 % and will call the function inv_model.R_prior
0010 % parameters to R_prior should be passed in the field
0011 % inv_model.R_prior_function_name.parameters
0012 %
0013 % If inv_model.R_prior is a matrix, calc_R_prior will return that matrix,
0014 % possibly correcting for coarse2fine
0015 %
0016 % R_prior      calculated regularization prior R
0017 % inv_model    is an inv_model structure
0018 
0019 % (C) 2005-2008 Andy Adler. License: GPL version 2 or version 3
0020 % $Id: calc_R_prior.m 5555 2017-06-18 17:07:12Z aadler $
0021 
0022 if ischar(inv_model) && strcmp(inv_model,'UNIT_TEST'); do_unit_test; return; end
0023 
0024 inv_model = rec_or_fwd_model( inv_model);
0025 
0026 
0027 if isfield(inv_model,'R_prior')
0028    if isnumeric(inv_model.R_prior)
0029       R_prior = inv_model.R_prior;
0030    else
0031       R_prior= eidors_cache( inv_model.R_prior, inv_model );
0032    end
0033 elseif isfield(inv_model,'RtR_prior')
0034    R_prior = eidors_cache(@calc_from_RtR_prior, inv_model,'calc_R_prior');
0035 else
0036    error('calc_R_prior: neither R_prior or RtR_prior func provided');
0037 end
0038 
0039 if isfield(inv_model.fwd_model,'coarse2fine')
0040    c2f= inv_model.fwd_model.coarse2fine;
0041    if size(R_prior,1)==size(c2f,1)
0042 %     we need to take into account coarse2fine - using a reasonable tol
0043       f2c = c2f';         
0044       R_prior = f2c*R_prior*c2f;
0045    end
0046 end
0047 
0048 function R_prior = calc_from_RtR_prior(inv_model)
0049    % The user has provided an RtR prior. We can use this to
0050    % get R =RtR^(1/2). Not that this is non unique
0051    if isnumeric(inv_model.RtR_prior)
0052       RtR_prior = inv_model.RtR_prior;
0053    else
0054       RtR_prior= eidors_cache( inv_model.RtR_prior, inv_model );
0055    end
0056    
0057    [L,D,P] = ldl(RtR_prior);
0058    R_prior = sqrt(D)*L'*P';
0059    [L,D,P] = ldl(RtR_prior);
0060 
0061 function inv_model = rec_or_fwd_model( inv_model);
0062 
0063    if isfield(inv_model,'rec_model');
0064       use_rec_model = 1;
0065       try if inv_model.prior_use_fwd_not_rec== 1;
0066          use_rec_model = 0;
0067       end; end
0068 
0069       if use_rec_model
0070           % copy the normalize flag from the fwd_model to prevent warnings
0071          inv_model.rec_model = mdl_normalize(inv_model.rec_model, ...
0072                                        mdl_normalize(inv_model.fwd_model));
0073          inv_model.fwd_model= inv_model.rec_model;
0074          inv_model= rmfield(inv_model,'rec_model');
0075       end
0076    end
0077 
0078 function do_unit_test
0079 
0080    priors={'prior_tikhonov',
0081            'prior_laplace',
0082            'prior_gaussian_HPF',
0083            'prior_movement'};
0084    models= {'a2c0','f2c2','n3r2'};
0085 
0086    for j= 1:length(models);
0087       imdl = mk_common_model(models{j},16);
0088       for i = 1:length(priors);
0089          imdl.RtR_prior = feval(priors{i},imdl);
0090 
0091          R   = calc_R_prior(imdl);
0092          RtR = calc_RtR_prior(imdl);  
0093          unit_test_cmp(['P:',priors{i}], R'*R, RtR, 1e-10);
0094       end
0095    end
0096 
0097

Generated on Wed 21-Jun-2017 09:29:07 by m2html © 2005