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
      - use to specify electrode errors
      - use to add movement (Jm * Jm')
   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 6242 2022-03-22 17:08:28Z aadler $

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 %      - use to specify electrode errors
0023 %      - use to add movement (Jm * Jm')
0024 %   options.desired_solution_fn
0025 %      specify a function to calculate the desired image.
0026 %      It must have the signature:
0027 %      D = my_function( xyc, radius, options);
0028 %      uses eidors_defualt('get','GREIT_desired_img') if not specified
0029 %
0030 %
0031 % (C) 2009 Andy Adler. Licenced under GPL v2 or v3
0032 % $Id: calc_GREIT_RM.m 6242 2022-03-22 17:08:28Z aadler $
0033 
0034    if ~isstruct(options)
0035        options.normalize = options;
0036    end
0037    opt = parse_options(options);
0038 
0039    if opt.normalize
0040       Y = vi./(vh*ones(1,size(vi,2))) - 1; 
0041    else
0042       Y = vi - (vh*ones(1,size(vi,2)));
0043    end
0044    try 
0045        f = opt.desired_solution_fn;
0046    catch
0047        f = eidors_default('get','GREIT_desired_img');
0048    end
0049 
0050    D = feval(f, xyc, radius, opt);
0051    
0052    % PJt is expensive and doesn't change when optimising NF
0053    copt.cache_obj = {vi,vh,xyc,radius,opt};
0054    copt.fstr = 'calc_GREIT_RM_PJt';
0055    PJt = eidors_cache(@calc_PJt,{Y,D},copt);
0056    
0057    if size(weight)==[1,1] % Can't use isscalar for compatibility with M6.5
0058        [RM, M] = calc_RM(PJt, Y, weight, opt);
0059    else
0060        error('not coded yet');
0061    end
0062 
0063 function [RM, M] = calc_RM(PJt, Y, noiselev, opt)
0064 
0065    noiselev = noiselev * mean(abs(Y(:)));
0066    % Desired soln for noise is 0
0067    N_meas = size(Y,1);
0068 
0069    % This implements RM = D*Y'/(J*Sx*J + Sn);
0070    Sn = speye(N_meas) * opt.noise_covar; % Noise covariance
0071 %    PJt= D*Y';
0072    M  = (Y*Y' + noiselev^2*Sn);
0073    RM =  left_divide(M',PJt')';    %PJt/M;
0074    % This implements RM = D*Y'/(Y*Y');
0075    if 0
0076       Y = [Y, noiselev*eye(N_meas)];
0077       D = [D,          zeros(size(D,1),N_meas)];
0078 
0079       RMold = D*Y'/(Y*Y');
0080       if norm(RM-RMold,'fro')/norm(RM,'fro') > 1e-10; warning('not OK'); end
0081    end
0082 
0083 function PJt = calc_PJt(Y,D)
0084       PJt = D*Y';
0085       
0086 
0087 function opt = parse_options(opt)
0088        if ~isfield(opt, 'normalize'), opt.normalize = 1; end
0089        if ~isfield(opt, 'meshsz'),    opt.meshsz = [-1 1 -1 1]; end
0090        if ~isfield(opt, 'imgsz'),     opt.imgsz = [32 32]; end
0091        if ~isfield(opt, 'noise_covar'),
0092                      opt.noise_covar = 1;
0093        end
0094 %   options.data_covariance [optional]