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
 
 See also: mk_GREIT_model

 (C) 2009 Andy Adler. Licenced under GPL v2 or v3
 $Id: calc_GREIT_RM.m 6681 2024-03-14 19:37:03Z aadler $

0001 function [RM, PJt, M, noiselev] = 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 % See also: mk_GREIT_model
0031 %
0032 % (C) 2009 Andy Adler. Licenced under GPL v2 or v3
0033 % $Id: calc_GREIT_RM.m 6681 2024-03-14 19:37:03Z aadler $
0034 
0035    if ischar(vh) && strcmp(vh,'UNIT_TEST'); do_unit_test; return; end
0036 
0037    if ~isstruct(options)
0038        options.normalize = options;
0039    end
0040    opt = parse_options(options);
0041 
0042    if opt.normalize
0043       Y = vi./(vh*ones(1,size(vi,2))) - 1; 
0044    else
0045       Y = vi - (vh*ones(1,size(vi,2)));
0046    end
0047    try 
0048        f = opt.desired_solution_fn;
0049    catch
0050        f = eidors_default('get','GREIT_desired_img');
0051    end
0052 
0053    D = feval(f, xyc, radius, opt);
0054    
0055    % PJt is expensive and doesn't change when optimising NF
0056    copt.cache_obj = {vi,vh,xyc,radius,opt};
0057    copt.fstr = 'calc_GREIT_RM_PJt';
0058    PJt = eidors_cache(@calc_PJt,{Y,D},copt);
0059    
0060    if isscalar(weight)
0061        [RM, M, noiselev] = calc_RM(PJt, Y, weight, opt);
0062    else
0063        error('use of weight matrix is not yet available');
0064    end
0065 
0066 function [RM, M, noiselev] = calc_RM(PJt, Y, noiselev, opt)
0067 
0068    noiselev = noiselev * mean(abs(Y(:)));
0069    % Desired soln for noise is 0
0070    N_meas = size(Y,1);
0071 
0072    % This implements RM = D*Y'/(J*Sx*J + Sn);
0073    Sn = speye(N_meas) * opt.noise_covar; % Noise covariance
0074 %    PJt= D*Y';
0075    % Conjugate transpose here
0076    M  = (Y*Y' + noiselev^2*Sn);
0077    % Ensure to use transpose not conjugate'
0078    RM =  left_divide(M.',PJt.').';    %PJt/M;
0079 
0080 function [RM, M] = calc_RM_testcode(PJt, Y, noiselev, opt)
0081 
0082    noiselev = noiselev * mean(abs(Y(:)));
0083    % Desired soln for noise is 0
0084    N_meas = size(Y,1);
0085 
0086    % This implements RM = D*Y'/(J*Sx*J + Sn);
0087    Sn = speye(N_meas) * opt.noise_covar; % Noise covariance
0088 %    PJt= D*Y';
0089    % Conjugate transpose here
0090    M  = (Y*Y' + noiselev^2*Sn);
0091    % OLD and TESTCODE
0092    Y = [Y, noiselev*eye(N_meas)];
0093 
0094    RM = PJt/(Y*Y');
0095 
0096 function PJt = calc_PJt(Y,D)
0097    PJt = D*Y';
0098       
0099 
0100 function opt = parse_options(opt)
0101    if ~isfield(opt, 'normalize'), opt.normalize = 1; end
0102    if ~isfield(opt, 'meshsz'),    opt.meshsz = [-1 1 -1 1]; end
0103    if ~isfield(opt, 'imgsz'),     opt.imgsz = [32 32]; end
0104    if ~isfield(opt, 'noise_covar'),
0105                  opt.noise_covar = 1;
0106    end
0107 %  options.data_covariance [optional]
0108 
0109 function do_unit_test
0110    img= mk_image(mk_common_model('a2c2',16));
0111    J = calc_jacobian(img);
0112    opt.noise_covar = 1;
0113    noiselev = 1;
0114    [RM, M] = calc_RM(J', J, noiselev, opt);
0115    img.elem_data = RM*J(:,20);
0116    show_fem(img); 
0117    unit_test_cmp('a2c2 noiselev=1', img.elem_data(19:21), ...
0118   [ -0.108227622116164; 0.737404299070634; -0.009842231726770], 1e-12);
0119 
0120    [RMold, M] = calc_RM_testcode(J', J, noiselev, opt);
0121    unit_test_cmp('test RM calc', RM, RMold, 1e-10);