calc_JiRtRJt

 [JiRtRJt,iRtRJt] = calc_JiRtRJt( inv_model, varargin )
 CALC_iRtR_PRIOR: calculate regularization matrix J*inv(R'*R)*J'
   This is a model of the covariance of image elements
   The image prior is matrix n_elem x n_elem 
 
 calc_JiRtRJt can be called as
    JiRtRJt= calc_JiRtRJt( inv_model, ... )

 and will call the function inv_model.JiRtRJt
 parameters to JiRtRJt should be passed in the field
 inv_model.JiRtRJt_function_name.parameters

 if JiRtRJt is a matrix, than calc_JiRtRJt will return that matrix

 JiRtRJt         calculated regularization prior R
 inv_model       is an inv_model structure
 inv_model.JiRtRJt_func function to make calculation

 TODO: think about how to implement this!!

0001 function [JiRtRJt,iRtRJt] = calc_JiRtRJt( inv_model, varargin )
0002 % [JiRtRJt,iRtRJt] = calc_JiRtRJt( inv_model, varargin )
0003 % CALC_iRtR_PRIOR: calculate regularization matrix J*inv(R'*R)*J'
0004 %   This is a model of the covariance of image elements
0005 %   The image prior is matrix n_elem x n_elem
0006 %
0007 % calc_JiRtRJt can be called as
0008 %    JiRtRJt= calc_JiRtRJt( inv_model, ... )
0009 %
0010 % and will call the function inv_model.JiRtRJt
0011 % parameters to JiRtRJt should be passed in the field
0012 % inv_model.JiRtRJt_function_name.parameters
0013 %
0014 % if JiRtRJt is a matrix, than calc_JiRtRJt will return that matrix
0015 %
0016 % JiRtRJt         calculated regularization prior R
0017 % inv_model       is an inv_model structure
0018 % inv_model.JiRtRJt_func function to make calculation
0019 %
0020 % TODO: think about how to implement this!!
0021 
0022 % (C) 2006 Andy Adler. License: GPL version 2 or version 3
0023 % $Id: calc_JiRtRJt.m 4832 2015-03-29 21:13:53Z bgrychtol-ipa $
0024 
0025 if isfield(inv_model,'JiRtRJt') 
0026    if isnumeric(inv_model.JiRtRJt)
0027       JiRtRJt = inv_model.JiRtRJt;
0028    else
0029       try inv_model.JiRtRJt = str2func(inv_model.JiRtRJt); end
0030       JiRtRJt = eidors_cache(inv_model.JiRtRJt,{inv_model});
0031    end
0032    return
0033 end
0034 
0035 % try to calculate from scratch
0036 JiRtRJt = eidors_cache(@calc_JiRtRJt_from_scratch,{inv_model},'calc_JiRtRJt');
0037 
0038 
0039 function JiRtRJt = calc_JiRtRJt_from_scratch(inv_model)
0040    eidors_msg( ...
0041       'calc_JiRtRJt: trying to calculate JiRtRJt from RtR_prior',2);
0042    RtR_prior = calc_RtR_prior(inv_model);
0043    % regularize slightly so inverse exists
0044    % This is very slow, but I can't think of a better idea
0045    RtR_p_reg = spdiags( spdiags(RtR_prior,0)*1.00001, 0, RtR_prior);
0046 
0047    img_bkgnd= calc_jacobian_bkgnd( inv_model );
0048    J = calc_jacobian( img_bkgnd);
0049 
0050    JiRtRJt= J*(RtR_p_reg\J');