calc_GREIT_RM

PURPOSE ^

CALCULATE GREIT reconstruction matrix

SYNOPSIS ^

function [RM, PJt, M] = calc_GREIT_RM(vh,vi, xyc, radius, weight, options)

DESCRIPTION ^

 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 4821 2015-03-29 14:36:39Z bgrychtol-ipa $

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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 4821 2015-03-29 14:36:39Z bgrychtol-ipa $
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    if size(weight)==[1,1] % Can't use isscalar for compatibility with M6.5
0051        [RM, PJt, M] = calc_RM(Y,D,weight, opt);
0052    else
0053        error('not coded yet');
0054    end
0055 
0056 function [RM, PJt, M] = calc_RM(Y, D, noiselev, opt)
0057 
0058    noiselev = noiselev * mean(abs(Y(:)));
0059    % Desired soln for noise is 0
0060    N_meas = size(Y,1);
0061 
0062    % This implements RM = D*Y'/(J*Sx*J + Sn);
0063    Sn = speye(N_meas) .* opt.noise_covar; % Noise covariance
0064    PJt= D*Y';
0065    M  = (Y*Y' + noiselev^2*Sn);
0066    RM =  left_divide(M',PJt')';    %PJt/M;
0067    % This implements RM = D*Y'/(Y*Y');
0068    if 0
0069       Y = [Y, noiselev*eye(N_meas)];
0070       D = [D,          zeros(size(D,1),N_meas)];
0071 
0072       RMold = D*Y'/(Y*Y');
0073       if norm(RM-RMold,'fro')/norm(RM,'fro') > 1e-10; warning('not OK'); end
0074    end
0075 
0076 
0077 
0078    function opt = parse_options(opt)
0079        if ~isfield(opt, 'normalize'), opt.normalize = 1; end
0080        if ~isfield(opt, 'meshsz'),    opt.meshsz = [-1 1 -1 1]; end
0081        if ~isfield(opt, 'imgsz'),     opt.imgsz = [32 32]; end
0082        if ~isfield(opt, 'noise_covar'),
0083                      opt.noise_covar = 1;
0084        end
0085 %   options.data_covariance [optional]

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