GREIT_desired_img_FEMmesh

GREIT_DESIRED_IMG_FEMmesh  GREIT onto a FEM mesh
 PSF= GREIT_desired_img_FEMmesh(xyc, radius, opt)
   xyc     - array of point centers [2xNpoints]
   radius  - the radius of the target on the desired image as a fraction
             of the model radius (half the larger dimension in xy)
   opt     - a struct with these mandatory fields:
      .rec_model   a 2D model as generated by MK_GRID_MODEL, with some
                   restrictions, see below. 

 Generates the GREIT "desired img" mapped onto the centres
   of the elements specified in opt.rec_model. It is assumed
   the radius is larger than the element size, so only the
   element centres are relevant

 See also: CALC_GREIT_RM, MK_GREIT_MODEL, MK_PIXEL_SLICE,
           GREIT_DESIRED_IMG_SIGMOID GREIT_DESIRED_IMG_ORIGINAL

0001 function PSF= GREIT_desired_img_FEMmesh(xyc, radius, opt)
0002 %GREIT_DESIRED_IMG_FEMmesh  GREIT onto a FEM mesh
0003 % PSF= GREIT_desired_img_FEMmesh(xyc, radius, opt)
0004 %   xyc     - array of point centers [2xNpoints]
0005 %   radius  - the radius of the target on the desired image as a fraction
0006 %             of the model radius (half the larger dimension in xy)
0007 %   opt     - a struct with these mandatory fields:
0008 %      .rec_model   a 2D model as generated by MK_GRID_MODEL, with some
0009 %                   restrictions, see below.
0010 %
0011 % Generates the GREIT "desired img" mapped onto the centres
0012 %   of the elements specified in opt.rec_model. It is assumed
0013 %   the radius is larger than the element size, so only the
0014 %   element centres are relevant
0015 %
0016 % See also: CALC_GREIT_RM, MK_GREIT_MODEL, MK_PIXEL_SLICE,
0017 %           GREIT_DESIRED_IMG_SIGMOID GREIT_DESIRED_IMG_ORIGINAL
0018 
0019 % (C) 2009-2015 Andy Adler and Bartlomiej Grychtol.
0020 % License: GPL v2 or v3
0021 % $Id: GREIT_desired_img_FEMmesh.m 5868 2018-12-15 12:00:35Z aadler $
0022    if ischar(xyc) && strcmp(xyc,'UNIT_TEST'); do_unit_test; return; end
0023 
0024    copt.fstr = 'GREIT_desired_img_FEMmesh';
0025    copt.cache_obj = {xyc, radius, opt};
0026    PSF = eidors_cache(@desired_soln,{xyc, radius, opt},copt);
0027 end
0028 
0029 function PSF = desired_soln(xyc,radius,opt)
0030    
0031    opt.inside = opt.rec_model.inside;
0032    mxyc = interp_mesh(opt.rec_model); x=mxyc(:,1); y=mxyc(:,2);
0033    xmin = opt.meshsz(1); xmax = opt.meshsz(2);
0034    ymin = opt.meshsz(3); ymax = opt.meshsz(4);
0035    % scale radius to half the greater dimension
0036    radius = radius * 0.5 * max(xmax-xmin, ymax-ymin);
0037    PSF = zeros(size(x,1),size(xyc,2));
0038    x_spc = (xmax-xmin)/(opt.imgsz(1)-1) * 0.5;
0039    y_spc = (ymax-ymin)/(opt.imgsz(2)-1) * 0.5;
0040    for i=1:size(xyc,2);
0041       for dx = linspace(-x_spc, x_spc, 5)
0042          for dy = linspace(-y_spc, y_spc, 5)
0043             PSF(:,i) = PSF(:,i) +  1/25*( ...
0044                (dx+x(:)-xyc(1,i)).^2 + (dy+y(:)-xyc(2,i)).^2 ...
0045                         < radius^2 );
0046          end
0047       end
0048 %     PSF(:,i) = PSF(:,i)/sum(PSF(:,i));
0049    end
0050    PSF = PSF(opt.inside,:);
0051 end
0052 
0053 function do_unit_test
0054    img=mk_image(mk_common_model('d2c2',32)); vh=fwd_solve(img);
0055    img.elem_data(206)=1.1;                    vi=fwd_solve(img);
0056 
0057    fm = ng_mk_cyl_models(0.2,[32,0.1],[0.05]);
0058    fm.stimulation = mk_stim_patterns(32,1,[0,1],[0,1],{},1);
0059    imdl= eidors_obj('inv_model','','fwd_model',fm);
0060    imdl.jacobian_bkgnd.value = 1;
0061    opt.desired_solution_fn = @GREIT_desired_img_original;
0062    imdl  = mk_GREIT_model(imdl,0.20,1e1,opt);
0063    imgr = inv_solve(imdl,vh,vi); subplot(121); show_fem(imgr);
0064    disp(calc_noise_figure(imdl,vh,vi))
0065 
0066    [~,idx] = sort(-imgr.elem_data);
0067    unit_test_cmp('[...]_img_original',idx(1:10), ...
0068          [365 366 333 334 398 397 364 367 332 399]');
0069 
0070    imdl.rec_model = img.fwd_model;
0071    imdl.rec_model.inside = true(num_elems(img),1);
0072    opt.desired_solution_fn = @GREIT_desired_img_FEMmesh;
0073    imdl  = mk_GREIT_model(imdl,0.20,1e1,opt);
0074    imgr = inv_solve(imdl,vh,vi); subplot(122); show_fem(imgr);
0075    disp(calc_noise_figure(imdl,vh,vi))
0076 
0077    [~,idx] = sort(-imgr.elem_data);
0078    unit_test_cmp('[...]_img_FEMmesh',idx(1:10), ...
0079          [235 206 267 236 178 301 205 268 152 300]');
0080 
0081 end