test_performance

 TEST_PERFORMANCE: test of difference reconstruction algorithms
   in terms of the performance parameters defined in GREIT

 Input:
   imdls = cell array of inverse models to test on
   fmdl  = fmdl of the shape to test on (default is 16 electrode tank).

 Model assumes the tank is elliptical with an x,y centre at 0,0

 [r, params] =  test_performance( imdls );
  r =  radii of test objects
  params(1,:) = Image Amplitude
  params(2,:) = Position Error => + toward centre, - toward edge
  params(3,:) = Resolution
  params(4,:) = Shape Deformation
  params(5,:) = Ringing
  params(6,:) = Image Amplitude in ROI (1/4 ampl)
  params(7,:) = Noise Figure

 Example
   i_bp = mk_common_gridmdl('backproj');
   i_gr = mk_common_gridmdl('GREITc1');
   [r,params] = test_performance( { i_bp, i_gr })

0001 function [r, params] =  test_performance( imdls, fmdl );
0002 % TEST_PERFORMANCE: test of difference reconstruction algorithms
0003 %   in terms of the performance parameters defined in GREIT
0004 %
0005 % Input:
0006 %   imdls = cell array of inverse models to test on
0007 %   fmdl  = fmdl of the shape to test on (default is 16 electrode tank).
0008 %
0009 % Model assumes the tank is elliptical with an x,y centre at 0,0
0010 %
0011 % [r, params] =  test_performance( imdls );
0012 %  r =  radii of test objects
0013 %  params(1,:) = Image Amplitude
0014 %  params(2,:) = Position Error => + toward centre, - toward edge
0015 %  params(3,:) = Resolution
0016 %  params(4,:) = Shape Deformation
0017 %  params(5,:) = Ringing
0018 %  params(6,:) = Image Amplitude in ROI (1/4 ampl)
0019 %  params(7,:) = Noise Figure
0020 %
0021 % Example
0022 %   i_bp = mk_common_gridmdl('backproj');
0023 %   i_gr = mk_common_gridmdl('GREITc1');
0024 %   [r,params] = test_performance( { i_bp, i_gr })
0025 
0026 % (c) 2010 Andy Adler. Licenced under GPL v2 or v3
0027 % $Id: test_performance.m 6887 2024-05-12 09:54:48Z aadler $
0028 
0029 if ischar(imdls) && strcmp(imdls,'UNIT_TEST'); do_unit_test; return; end
0030 
0031 if nargin==1
0032    fmdl = ng_mk_cyl_models([2,1,0.05],[16,1],[0.05]); 
0033    fmdl.stimulation = mk_stim_patterns(16,1,[0,1],[0,1],{},1);
0034 end
0035 imgs= mk_image( fmdl, 1);
0036 
0037 if ~iscell(imdls)
0038     imdls = {imdls};
0039 end
0040 
0041 Nsim = 100;
0042 r =  linspace(0,0.9,Nsim);
0043 ctr = [0,0, mean(fmdl.nodes(:,3))];  % Assume x,y centre is zero
0044 maxx = max(abs(fmdl.nodes(:,1) - ctr(1)));
0045 maxy = max(abs(fmdl.nodes(:,2) - ctr(2)));
0046 
0047 th=  r*4321; % want object to jump around in radius
0048 xyzr = [maxx*r.*cos(th); maxy*r.*sin(th); ctr(3)*ones(1,Nsim); 
0049         0.050*mean([maxx,maxy])*ones(1, Nsim)];
0050 
0051 [vh,vi] = simulate_movement(imgs, xyzr);
0052 for i= 1:length(imdls);
0053    imgr = inv_solve(imdls{i}, vh, vi);
0054    imgr.calc_colours.npoints = 64;
0055    param_GR = eval_GREIT_fig_merit(imgr, xyzr);
0056    pnoise = calc_noise_figure( imdls{i}, vh, vi );
0057    params{i} = [param_GR; pnoise];
0058 end
0059 
0060 clf; Nparams = 6;
0061 for j=1:Nparams;
0062    p = [];
0063    for i= 1:length(params);
0064       p = [p, params{i}(j,:)'];
0065    end
0066    m = mean([maxx,maxy]);
0067 %   subplot(5,1,j);
0068    hig = 0.95/Nparams;
0069    axes('position',[0.1,0.05 + hig*(Nparams-j),.88, hig]);
0070    plot(r, p);
0071    if j==1;     axis([0,0.9,0,2.1]);        ylabel('AR');
0072    elseif j==2; axis([0,0.9,-0.16*m,0.16*m]);   ylabel('PE');
0073    elseif j==3; axis([0,0.9,0,0.41]);       ylabel('RES');
0074    elseif j==4; axis([0,0.9,0,0.31]);       ylabel('SD');
0075    elseif j==5; axis([0,0.9,0,0.61]);       ylabel('RNG');
0076    elseif j==6; set(gca,'YScale','log');    ylabel('NF');
0077    end
0078    if j<Nparams; set(gca,'XTickLabel',[]);end
0079 end
0080 
0081 function do_unit_test;
0082 % Reconstruct GREIT Images
0083 imdl_gr = mk_common_gridmdl('GREITc1');
0084 
0085 % Reconstruct backprojection Images
0086 imdl_bp = mk_common_gridmdl('backproj');
0087 
0088 % Reconstruct GN Images
0089 imdl_gn = select_imdl( mk_common_model('d2c2', 16), {'Basic GN dif','Choose NF=0.5'});
0090 
0091 test_performance( { imdl_gr, imdl_bp, imdl_gn } );
0092 
0093 
0094 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0095 function [vh,vi,param]=test_performance_old(fmdl,imdl,r,N)
0096 
0097 debug = false;
0098 
0099 %% 1. Figure out the limits of fmdl
0100 boundary = fmdl.boundary;
0101 b_ver = unique(boundary);
0102 boundary = fmdl.nodes(b_ver,:);
0103 min_b = min(boundary)+repmat(0.5*r,1,3);
0104 max_b = max(boundary)-repmat(0.5*r,1,3);
0105 vol = get_elem_volume(fmdl);
0106 
0107 %% 2. create N samples within the volume
0108 rand('state', sum(100*clock));
0109 
0110 n_elems = size(fmdl.elems,1);
0111 
0112 
0113 c2f = eidors_obj('get-cache', {fmdl, N, r}, 'samples');
0114 xyzr_pt = eidors_obj('get-cache', {fmdl, N, r}, 'points');
0115 
0116 count = 0; %count points already created
0117 
0118 if isempty(c2f) || isempty(xyzr_pt);
0119     disp('Generating sample points.');
0120     c2f = sparse(n_elems,N);
0121     while count < N
0122         n_pts = ceil(1.5 * (N - count) ); % more than needed to limit iterations
0123         vec = rand_pt ( repmat(min_b,n_pts,1), repmat(max_b,n_pts,1) );
0124         map = mk_c2f_circ_mapping(fmdl, [vec repmat(r,n_pts,1)]');
0125         frac = vol'*map/(4/3*pi*r^3);
0126         idx = find(frac > 0.9, N - count, 'first');
0127         
0128         xyzr_pt = [xyzr_pt, vec(idx,:)'];
0129         newcount = count + length(idx);
0130         c2f(:, (count+1):newcount ) = map(:,idx);
0131         count = newcount;
0132     end
0133     eidors_obj('set-cache', {fmdl, N, r}, 'samples', c2f);
0134     eidors_obj('set-cache', {fmdl, N, r}, 'points', xyzr_pt );
0135     disp('Done.');
0136 end
0137  
0138 
0139 
0140 % x=110;y=150;r=20;
0141 % x=110;y=80;r=20;
0142 % map = mk_c2f_circ_mapping(fmdl, [x,y,r]');
0143 % img= mk_image(map);img.fwd_model = fmdl;
0144 % phi= linspace(0,2*pi,20);xr=r*cos(phi)+x; yr=r*sin(phi)+y;
0145 % show_fem(img,[0,1,1]); hold on; plot(xr,yr,'r');hold off
0146 % vol = get_elem_volume(fmdl); vol'*map/(pi*r^2)
0147 
0148 %% 3. reconstruct the images
0149 % create a homegenous image
0150    elem_data = ones(size(fmdl.elems,1),1);
0151    img = mk_image(fmdl, elem_data);
0152    disp('Calculating Jacobian. This may take a (long) while.');
0153    J = calc_jacobian( img );
0154    disp('Done.');
0155    vh= fwd_solve(img);
0156    vh = vh.meas;
0157    vi= vh*ones(1,N) + J*c2f;
0158 
0159 %% 4. calculate parameters
0160 
0161 
0162    
0163    img = inv_solve(imdl,vh,vi);
0164    img.calc_colours.npoints=32;
0165    imgs = calc_slices(img);
0166 
0167        param= [calc_noise_figure(imdl, vh, vi); ... % noise parameters
0168            GREIT_sim_params(  imgs, xyzr_pt)];            % image parameters
0169        
0170        r = sqrt(sum(xyzr_pt(1:2,:).^2));
0171        
0172        names = {'Noise','Amplitude','Position Error','Resolution', ...
0173            'Shape Deformation', 'Ringing'};
0174        for j=1:size(param,1)
0175            figure
0176            scatter(r, param(j,:));
0177            title(names{j});
0178        end
0179 
0180        
0181        
0182        
0183 function x = rand_pt(x_min, x_max)
0184 x = rand(size(x_min)) .* (x_max - x_min) + x_min;