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
  param = [AR, PE, RES, SD, RNG, NOI]

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