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
 

 (C) 2009 Andy Adler. Licenced under GPL v2 or v3
 $Id: calc_GREIT_RM.html 2819 2011-09-07 16:43:11Z aadler $

0001 function RM= 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 %
0021 %
0022 % (C) 2009 Andy Adler. Licenced under GPL v2 or v3
0023 % $Id: calc_GREIT_RM.html 2819 2011-09-07 16:43:11Z aadler $
0024 
0025    if ~isstruct(options)
0026        options.normalize = options;
0027    end
0028    opt = parse_options(options);
0029 
0030    if opt.normalize
0031       Y = vi./(vh*ones(1,size(vi,2))) - 1; 
0032    else
0033       Y = vi - (vh*ones(1,size(vi,2)));
0034    end
0035 
0036    D = desired_soln( xyc, radius, opt);
0037 
0038    if size(weight)==[1,1] % Can't use isscalar for compatibility with M6.5
0039        [RM] = calc_RM(Y,D,weight);
0040    else
0041        error('not coded yet');
0042    end
0043 
0044 function RM = calc_RM(Y, D, noiselev)
0045 
0046    noiselev = noiselev * mean(abs(Y(:)));
0047    % Desired soln for noise is 0
0048    N_meas = size(Y,1);
0049    Y = [Y, noiselev*eye(N_meas)];
0050    D = [D,          zeros(size(D,1),N_meas)];
0051 
0052    RM = D*Y'/(Y*Y');
0053 
0054 function PSF= desired_soln(xyc, radius, opt)
0055    c_obj = {xyc, radius, opt};
0056    PSF = eidors_obj('get-cache', c_obj, 'desired_solution');
0057    if ~isempty(PSF)
0058        return
0059    end
0060         
0061    xsz = opt.imgsz(1); ysz = opt.imgsz(2);
0062    sz= xsz * ysz;
0063    xmin = opt.meshsz(1); xmax = opt.meshsz(2);
0064    ymin = opt.meshsz(3); ymax = opt.meshsz(4);
0065    % scale radius to half the greater dimension
0066    radius = radius * 0.5 * max(xmax-xmin, ymax-ymin);
0067    [x,y]= ndgrid(linspace(xmin,xmax,xsz), linspace(ymin,ymax,ysz));
0068    x_spc = (xmax-xmin)/(xsz-1) * 0.5;
0069    y_spc = (ymax-ymin)/(ysz-1) * 0.5;
0070    PSF = zeros(sz,size(xyc,2));
0071    for i=1:size(xyc,2);
0072       for dx = linspace(-x_spc, x_spc, 5)
0073          for dy = linspace(-y_spc, y_spc, 5)
0074             PSF(:,i) = PSF(:,i) +  1/25*( ...
0075                (dx+x(:)-xyc(1,i)).^2 + (dy+y(:)-xyc(2,i)).^2 ...
0076                         < radius^2 );
0077          end
0078       end
0079 %     PSF(:,i) = PSF(:,i)/sum(PSF(:,i));
0080    end
0081    eidors_obj('set-cache', c_obj, 'desired_solution', PSF);
0082 
0083    function opt = parse_options(opt)
0084        if ~isfield(opt, 'normalize'), opt.normalize = 1; end
0085        if ~isfield(opt, 'meshsz'),    opt.meshsz = [-1 1 -1 1]; end
0086        if ~isfield(opt, 'imgsz'),     opt.imgsz = [32 32]; end