test_performance_img

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

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


 params_img =  test_performance( imdls );
  params_img = {AR, PE, RES, SD, RNG, NOI}
  is a cell array of 3D (Np x Np x N_models) arrays, each representing 
  the given perfomance metrics as an image computed using regularly spaced
  targets.
 

 Example
   i_bp = mk_common_gridmdl('backproj');
   i_gp = mk_common_gridmdl('GREITc1');
   test_performance( { i_gp, i_gr })

0001 function [params_img] =  test_performance_img( imdls, fmdl );
0002 % TEST_PERFORMANCE: test of difference reconstruction algorithms
0003 %   in terms of the performance parameters defined in GREIT, image output
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 %
0010 % params_img =  test_performance( imdls );
0011 %  params_img = {AR, PE, RES, SD, RNG, NOI}
0012 %  is a cell array of 3D (Np x Np x N_models) arrays, each representing
0013 %  the given perfomance metrics as an image computed using regularly spaced
0014 %  targets.
0015 %
0016 %
0017 % Example
0018 %   i_bp = mk_common_gridmdl('backproj');
0019 %   i_gp = mk_common_gridmdl('GREITc1');
0020 %   test_performance( { i_gp, i_gr })
0021 
0022 % (c) 2011 Andy Adler and Barlomiej Grychtol. Licenced under GPL v2 or v3
0023 % $Id: test_performance_img.m 5112 2015-06-14 13:00:41Z aadler $
0024 
0025 if ischar(imdls) && strcmp(imdls,'UNIT_TEST'); do_unit_test; return; end
0026 
0027 if nargin==1
0028    fmdl = ng_mk_cyl_models([2,1,0.05],[16,1],[0.05]); 
0029    fmdl.stimulation = mk_stim_patterns(16,1,[0,1],[0,1],{},1);
0030 end
0031 switch fmdl.type
0032     case 'fwd_model'
0033         imgs= mk_image( fmdl, 1);
0034     case 'image'
0035         imgs = fmdl;
0036         fmdl = fmdl.fwd_model;
0037     otherwise
0038         error('Second parameter not understood');
0039 end
0040 % if ~iscell(imdls)
0041 %     imdls = {imdls};
0042 % end
0043 
0044 if (0) % old code
0045 Nsim = 100;
0046 r =  linspace(0,0.9,Nsim);
0047 ctr = [0,0, mean(fmdl.nodes(:,3))];  % Assume x,y centre is zero
0048 maxx = max(abs(fmdl.nodes(:,1) - ctr(1)));
0049 maxy = max(abs(fmdl.nodes(:,2) - ctr(2)));
0050 
0051 th=  r*4321; % want object to jump around in radius
0052 xyzr = [maxx*r.*cos(th); maxy*r.*sin(th); ctr(3)*ones(1,Nsim); 
0053         0.05/mean([maxx,maxy])*ones(1, Nsim)];
0054 else 
0055 %%
0056     N = 128;
0057      % this bit is copied from mdl_slice_mapper to make sure
0058      xmin = min(fmdl.nodes(:,1));    xmax = max(fmdl.nodes(:,1));
0059      xmean= mean([xmin,xmax]); xrange= xmax-xmin;
0060 
0061      ymin = min(fmdl.nodes(:,2));    ymax = max(fmdl.nodes(:,2));
0062      ymean= mean([ymin,ymax]); yrange= ymax-ymin;
0063 
0064      range= max([xrange, yrange]);
0065      xspace = linspace( xmean - range*0.5, xmean + range*0.5, N );
0066      yspace = linspace( ymean + range*0.5, ymean - range*0.5, N );
0067      % end of copied block
0068     bnd_nodes = unique(fmdl.boundary);
0069     min_bb = min(fmdl.nodes(bnd_nodes,:));
0070     max_bb = max(fmdl.nodes(bnd_nodes,:));
0071     h = mean([min_bb(3),max_bb(3)]);
0072     r = 0.025 * range;
0073 %     xspace = linspace(min_bb(1),max_bb(1),N);
0074 %     yspace = linspace(max_bb(2),min_bb(2),N);
0075     [X Y] = meshgrid(xspace,yspace);
0076     img = mk_image(fmdl,1);
0077     img.calc_colours.npoints = N;
0078     M = calc_slices(img,1);
0079     IN = M==1;
0080     n_px = nnz(IN);
0081     OUT = isnan(M);
0082     xyzr = [X(IN)'; Y(IN)'; h*ones(1,n_px);  r*ones(1,n_px); ];
0083 %%
0084 end
0085     
0086     
0087 [vh,vi] = simulate_movement(imgs, xyzr);
0088 IN6 = repmat(IN,[1 1 6]); IN6 = shiftdim(IN6,2);
0089 for i= 1:length(imdls);
0090     if iscell(imdls)
0091         imgr = inv_solve(imdls{i}, vh, vi);
0092         pnoise = calc_noise_figure( imdls{i}, vh, vi );
0093     else
0094         imgr = inv_solve(imdls(i), vh, vi);
0095         pnoise = calc_noise_figure( imdls(i), vh, vi );
0096     end
0097    imgr.calc_colours.npoints = N;
0098    param_GR = eval_GREIT_fig_merit(imgr, xyzr);
0099    params{i} = [param_GR; pnoise];
0100    p_img = nan([6 size(IN)]);
0101    p_img(IN6) = params{i};
0102    params_img{i}=p_img;
0103 end
0104 
0105 
0106 
0107 
0108 return
0109 clf; Nparams = 6;
0110 for j=1:Nparams;
0111    p = [];
0112    for i= 1:length(params);
0113       p = [p, params{i}(j,:)'];
0114    end
0115 
0116 %   subplot(5,1,j);
0117    hig = 0.95/Nparams;
0118    axes('position',[0.1,0.05 + hig*(Nparams-j),.88, hig]);
0119    plot(r, p);
0120    if j==1;     axis([0,0.9,0,2.1]);        ylabel('AR');
0121    elseif j==2; axis([0,0.9,-0.16,0.16]);   ylabel('PE');
0122    elseif j==3; axis([0,0.9,0,0.41]);       ylabel('RES');
0123    elseif j==4; axis([0,0.9,0,0.31]);       ylabel('SD');
0124    elseif j==5; axis([0,0.9,0,0.61]);       ylabel('RNG');
0125    elseif j==6; set(gca,'YScale','log');    ylabel('NF');
0126    end
0127    if j<Nparams; set(gca,'XTickLabel',[]);end
0128 end
0129 
0130 
0131 
0132 
0133 
0134 
0135 
0136 
0137 
0138 function do_unit_test;
0139 % Reconstruct GREIT Images
0140 imdl_gr = mk_common_gridmdl('GREITc1');
0141 
0142 % Reconstruct backprojection Images
0143 imdl_bp = mk_common_gridmdl('backproj');
0144 
0145 % Reconstruct GN Images
0146 imdl_gn = select_imdl( mk_common_model('d2c2', 16), {'Basic GN dif','Choose NF=0.5'});
0147 
0148 test_performance( { imdl_gr, imdl_bp, imdl_gn } );
0149 
0150 
0151 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0152 function [vh,vi,param]=test_performance_old(fmdl,imdl,r,N)
0153 
0154 debug = false;
0155 
0156 %% 1. Figure out the limits of fmdl
0157 boundary = fmdl.boundary;
0158 b_ver = unique(boundary);
0159 boundary = fmdl.nodes(b_ver,:);
0160 min_b = min(boundary)+repmat(0.5*r,1,3);
0161 max_b = max(boundary)-repmat(0.5*r,1,3);
0162 vol = get_elem_volume(fmdl);
0163 
0164 %% 2. create N samples within the volume
0165 rand('state', sum(100*clock));
0166 
0167 n_elems = size(fmdl.elems,1);
0168 
0169 
0170 c2f = eidors_obj('get-cache', {fmdl, N, r}, 'samples');
0171 xyzr_pt = eidors_obj('get-cache', {fmdl, N, r}, 'points');
0172 
0173 count = 0; %count points already created
0174 
0175 if isempty(c2f) || isempty(xyzr_pt);
0176     disp('Generating sample points.');
0177     c2f = sparse(n_elems,N);
0178     while count < N
0179         n_pts = ceil(1.5 * (N - count) ); % more than needed to limit iterations
0180         vec = rand_pt ( repmat(min_b,n_pts,1), repmat(max_b,n_pts,1) );
0181         map = mk_c2f_circ_mapping(fmdl, [vec repmat(r,n_pts,1)]');
0182         frac = vol'*map/(4/3*pi*r^3);
0183         idx = find(frac > 0.9, N - count, 'first');
0184         
0185         xyzr_pt = [xyzr_pt, vec(idx,:)'];
0186         newcount = count + length(idx);
0187         c2f(:, (count+1):newcount ) = map(:,idx);
0188         count = newcount;
0189     end
0190     eidors_obj('set-cache', {fmdl, N, r}, 'samples', c2f);
0191     eidors_obj('set-cache', {fmdl, N, r}, 'points', xyzr_pt );
0192     disp('Done.');
0193 end
0194  
0195 
0196 
0197 % x=110;y=150;r=20;
0198 % x=110;y=80;r=20;
0199 % map = mk_c2f_circ_mapping(fmdl, [x,y,r]');
0200 % img= mk_image(map);img.fwd_model = fmdl;
0201 % phi= linspace(0,2*pi,20);xr=r*cos(phi)+x; yr=r*sin(phi)+y;
0202 % show_fem(img,[0,1,1]); hold on; plot(xr,yr,'r');hold off
0203 % vol = get_elem_volume(fmdl); vol'*map/(pi*r^2)
0204 
0205 %% 3. reconstruct the images
0206 % create a homegenous image
0207    elem_data = ones(size(fmdl.elems,1),1);
0208    img = mk_image(fmdl, elem_data);
0209    disp('Calculating Jacobian. This may take a (long) while.');
0210    J = calc_jacobian( img );
0211    disp('Done.');
0212    vh= fwd_solve(img);
0213    vh = vh.meas;
0214    vi= vh*ones(1,N) + J*c2f;
0215 
0216 %% 4. calculate parameters
0217 
0218 
0219    
0220    img = inv_solve(imdl,vh,vi);
0221    img.calc_colours.npoints=32;
0222    imgs = calc_slices(img);
0223 
0224        param= [calc_noise_figure(imdl, vh, vi); ... % noise parameters
0225            GREIT_sim_params(  imgs, xyzr_pt)];            % image parameters
0226        
0227        r = sqrt(sum(xyzr_pt(1:2,:).^2));
0228        
0229        names = {'Noise','Amplitude','Position Error','Resolution', ...
0230            'Shape Deformation', 'Ringing'};
0231        for j=1:size(param,1)
0232            figure
0233            scatter(r, param(j,:));
0234            title(names{j});
0235        end
0236 
0237        
0238        
0239        
0240 function x = rand_pt(x_min, x_max)
0241 x = rand(size(x_min)) .* (x_max - x_min) + x_min;