calc_GREIT_RM

 CALCULATE GREIT reconstruction matrix
   RM= calc_GREIT_RM(vh,vi, xyc, radius, weight, normalize)
 
 Input:
   vh     = homogeneous (reference) training measurements 
   vi     = inhomogeneous training measurements 
   xyc    = x,y position of targets 2xN
   radius = requested resolution matrix
      if scalar - apply resolution to all images:  recommend 0.25 for 16 elecs
   weight = weighting matrix
      if scalar   = weighting of noise vs signal
      if 32^2 x N = weighting of each image output
   options.normalize = 
            0 -> regular difference EIT
            1 -> normalized difference EIT
   options.meshsz = [xmin xmax ymin ymax] size of mesh
      defaults to [-1 1 -1 1]
   options.imgsz  = [xsz ysz] size of the reconstructed image in pixels
   options.noise_covar [optional]
      covariance matrix of data noise
   options.desired_solution_fn
      specify a function to calculate the desired image. 
      It must have the signature:
      D = my_function( xyc, radius, options);
      uses eidors_defualt('get','GREIT_desired_img') if not specified
 

 (C) 2009 Andy Adler. Licenced under GPL v2 or v3
 $Id: calc_GREIT_RM.m 5662 2017-12-03 12:58:59Z bgrychtol $

0001 function [RM, PJt, M] = calc_GREIT_RM(vh,vi, xyc, radius, weight, options)
0002 % CALCULATE GREIT reconstruction matrix
0003 %   RM= calc_GREIT_RM(vh,vi, xyc, radius, weight, normalize)
0004 %
0005 % Input:
0006 %   vh     = homogeneous (reference) training measurements
0007 %   vi     = inhomogeneous training measurements
0008 %   xyc    = x,y position of targets 2xN
0009 %   radius = requested resolution matrix
0010 %      if scalar - apply resolution to all images:  recommend 0.25 for 16 elecs
0011 %   weight = weighting matrix
0012 %      if scalar   = weighting of noise vs signal
0013 %      if 32^2 x N = weighting of each image output
0014 %   options.normalize =
0015 %            0 -> regular difference EIT
0016 %            1 -> normalized difference EIT
0017 %   options.meshsz = [xmin xmax ymin ymax] size of mesh
0018 %      defaults to [-1 1 -1 1]
0019 %   options.imgsz  = [xsz ysz] size of the reconstructed image in pixels
0020 %   options.noise_covar [optional]
0021 %      covariance matrix of data noise
0022 %   options.desired_solution_fn
0023 %      specify a function to calculate the desired image.
0024 %      It must have the signature:
0025 %      D = my_function( xyc, radius, options);
0026 %      uses eidors_defualt('get','GREIT_desired_img') if not specified
0027 %
0028 %
0029 % (C) 2009 Andy Adler. Licenced under GPL v2 or v3
0030 % $Id: calc_GREIT_RM.m 5662 2017-12-03 12:58:59Z bgrychtol $
0031 
0032    if ~isstruct(options)
0033        options.normalize = options;
0034    end
0035    opt = parse_options(options);
0036 
0037    if opt.normalize
0038       Y = vi./(vh*ones(1,size(vi,2))) - 1; 
0039    else
0040       Y = vi - (vh*ones(1,size(vi,2)));
0041    end
0042    try 
0043        f = opt.desired_solution_fn;
0044    catch
0045        f = eidors_default('get','GREIT_desired_img');
0046    end
0047 
0048    D = feval(f, xyc, radius, opt);
0049    
0050    % PJt is expensive and doesn't change when optimising NF
0051    copt.cache_obj = {vi,vh,xyc,radius,opt};
0052    copt.fstr = 'calc_GREIT_RM_PJt';
0053    PJt = eidors_cache(@calc_PJt,{Y,D},copt);
0054    
0055    if size(weight)==[1,1] % Can't use isscalar for compatibility with M6.5
0056        [RM, M] = calc_RM(PJt, Y, weight, opt);
0057    else
0058        error('not coded yet');
0059    end
0060 
0061 function [RM, M] = calc_RM(PJt, Y, noiselev, opt)
0062 
0063    noiselev = noiselev * mean(abs(Y(:)));
0064    % Desired soln for noise is 0
0065    N_meas = size(Y,1);
0066 
0067    % This implements RM = D*Y'/(J*Sx*J + Sn);
0068    Sn = speye(N_meas) .* opt.noise_covar; % Noise covariance
0069 %    PJt= D*Y';
0070    M  = (Y*Y' + noiselev^2*Sn);
0071    RM =  left_divide(M',PJt')';    %PJt/M;
0072    % This implements RM = D*Y'/(Y*Y');
0073    if 0
0074       Y = [Y, noiselev*eye(N_meas)];
0075       D = [D,          zeros(size(D,1),N_meas)];
0076 
0077       RMold = D*Y'/(Y*Y');
0078       if norm(RM-RMold,'fro')/norm(RM,'fro') > 1e-10; warning('not OK'); end
0079    end
0080 
0081 function PJt = calc_PJt(Y,D)
0082       PJt = D*Y';
0083       
0084 
0085 function opt = parse_options(opt)
0086        if ~isfield(opt, 'normalize'), opt.normalize = 1; end
0087        if ~isfield(opt, 'meshsz'),    opt.meshsz = [-1 1 -1 1]; end
0088        if ~isfield(opt, 'imgsz'),     opt.imgsz = [32 32]; end
0089        if ~isfield(opt, 'noise_covar'),
0090                      opt.noise_covar = 1;
0091        end
0092 %   options.data_covariance [optional]