eval_GREIT_fig_merit

 EVAL_GREIT_FIG_MERIT: calculate GREIT figures of merit for images
 USAGE:
 params = eval_GREIT_fig_merit(imgs, xyzr_pt)
  params(1,:) = Image Amplitude
  params(2,:) = Position Error => + toward centre, - toward edge
  params(3,:) = Resolution
  params(4,:) = Shape Deformation
  params(5,:) = Ringing

  imgs:    a sequence of eidors images of single point targets
  xyzr_pt: [x;y;z;radius] of each images point

0001 function params = eval_GREIT_fig_merit(imgs, xyzr_pt)
0002 % EVAL_GREIT_FIG_MERIT: calculate GREIT figures of merit for images
0003 % USAGE:
0004 % params = eval_GREIT_fig_merit(imgs, xyzr_pt)
0005 %  params(1,:) = Image Amplitude
0006 %  params(2,:) = Position Error => + toward centre, - toward edge
0007 %  params(3,:) = Resolution
0008 %  params(4,:) = Shape Deformation
0009 %  params(5,:) = Ringing
0010 %
0011 %  imgs:    a sequence of eidors images of single point targets
0012 %  xyzr_pt: [x;y;z;radius] of each images point
0013 
0014 % (C) 2009 Andy Adler. Licensed under GPL v2 or v3
0015 % $Id: eval_GREIT_fig_merit.html 2819 2011-09-07 16:43:11Z aadler $
0016 mdl = imgs.fwd_model;
0017 imgs = calc_slices(imgs);
0018 map = ~isnan(squeeze(imgs(:,:,1))); %assume all imgs are the same shape
0019 imgs(isnan(imgs)) = 0;
0020 sz = size(imgs,1);
0021 [x,y,bb_min,bb_max]=prepare_grid(sz,mdl);
0022 
0023 N_imgs = size(imgs,3);
0024 for i= 1:N_imgs
0025    [xmean,ymean,equiv_circ,qmi,img] = calc_cofg(imgs(:,:,i),map,x,y);
0026    params(1,i) = calc_amplitude( img );
0027    params(2,i) = calc_posn_error( qmi, xmean, ymean, xyzr_pt(1:2,i) );
0028    params(3,i) = calc_resolution( qmi, map );
0029    params(4,i) = calc_shape_deform( qmi, equiv_circ );
0030    params(5,i) = calc_ringing( img, qmi );
0031 end
0032 
0033 % TODO: Fix this when we start to care about units
0034 ctr = bb_min + 0.5*(bb_max-bb_min);
0035 r = max(0.5*(bb_max-bb_min));
0036 if N_imgs > 10 % doesn't make sense to normalize otherwise
0037     ctr_pts = sum((xyzr_pt(1:mdl_dim(mdl(1)),:)-repmat(ctr',1,size(xyzr_pt,2))).^2) < (0.05*r)^2;
0038     if any(ctr_pts)
0039         params(1,:) = params(1,:)/mean(params(1,ctr_pts));
0040     else
0041         eidors_msg('eval_GREIT_fig_merit: no centre points found to normalize',1);
0042     end
0043 end
0044 
0045 
0046 
0047 
0048 
0049 function ampl = calc_amplitude(img)
0050    ampl = sum(img(:));
0051 
0052 function pe   = calc_posn_error(qmi, xmean, ymean, xy)
0053    pe = sqrt(sum(xy.^2)) - sqrt( xmean^2 + ymean^2);
0054 
0055 function res  = calc_resolution(qmi, map)
0056    res = sqrt( sum(qmi(:)) / sum(map(:)));
0057 
0058 function sd  = calc_shape_deform(qmi, equiv_circ)
0059    not_circ= qmi & ~equiv_circ;
0060    sd = sum(not_circ(:))/sum(qmi(:));
0061 
0062 function rr = calc_ringing(img, qmi );
0063    ring_part =  img .* ( (img<0) & ~qmi);
0064    rr = -sum( ring_part(:) )/sum( img(:).*qmi(:) );
0065    
0066 function [x,y,bb_min,bb_max]=prepare_grid(sz,mdl)
0067       % bounding box
0068    bnd = unique(mdl.boundary);
0069    bb_min = min(mdl.nodes(bnd,:));
0070    bb_max = max(mdl.nodes(bnd,:));
0071    
0072    [x,y]=meshgrid(linspace(bb_min(1),bb_max(1),sz),linspace(bb_min(2),bb_max(2),sz)); 
0073    
0074 function [xmean,ymean,equiv_circ,qmi,img] = calc_cofg(img,map,x,y);
0075 %  if abs(max(img(:))) < abs(min(img(:))); img= -img; end
0076    sz = size(img,1);
0077    qmi = calc_hm_set( img, 0.25 );
0078    if sum(img(:) & qmi(:))<0 ; keyboard ; end
0079    
0080    pix_sz = range(x(:))*range(y(:))/numel(img);
0081 
0082    %map = x.^2+y.^2<1.1;
0083    qmi = qmi.*map; img = img.*map;
0084 
0085    ss_qmi = sum(qmi(:));
0086    xmean =  sum(sum( (qmi.*x) ))/ss_qmi; % centre of gravity
0087    ymean =  sum(sum( (qmi.*y) ))/ss_qmi;
0088    equiv_circ = (x-xmean).^2 + (y-ymean).^2 < pix_sz*ss_qmi/pi;