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
  params(6,:) = Image Amplitude in ROI (1/4 ampl)

  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 %  params(6,:) = Image Amplitude in ROI (1/4 ampl)
0011 %
0012 %  imgs:    a sequence of eidors images of single point targets
0013 %  xyzr_pt: [x;y;z;radius] of each images point
0014 
0015 % (C) 2009 Andy Adler. Licensed under GPL v2 or v3
0016 % $Id: eval_GREIT_fig_merit.m 6414 2022-11-28 21:09:12Z aadler $
0017 mdl = imgs.fwd_model;
0018 imgs = calc_slices(imgs);
0019 map = ~isnan(squeeze(imgs(:,:,1))); %assume all imgs are the same shape
0020 imgs(isnan(imgs)) = 0;
0021 sz = size(imgs);
0022 [x,y,bb_min,bb_max]=prepare_grid(sz,mdl);
0023 
0024 N_imgs = size(imgs,3);
0025 for i= 1:N_imgs
0026    [xmean,ymean,equiv_circ,qmi,img] = calc_cofg(imgs(:,:,i),map,x,y);
0027    params(1,i) = calc_amplitude( img );
0028    params(6,i) = calc_amplitude( img, qmi );
0029    params(2,i) = calc_posn_error( qmi, xmean, ymean, xyzr_pt(1:2,i) );
0030    params(3,i) = calc_resolution( qmi, map );
0031    params(4,i) = calc_shape_deform( qmi, equiv_circ );
0032    params(5,i) = calc_ringing( img, qmi );
0033 end
0034 
0035 % TODO: Fix this when we start to care about units
0036 ctr = bb_min + 0.5*(bb_max-bb_min);
0037 r = max(0.5*(bb_max-bb_min));
0038 if N_imgs > 10 % doesn't make sense to normalize otherwise
0039     ctr_pts = sum((xyzr_pt(1:mdl_dim(mdl(1)),:)-repmat(ctr',1,size(xyzr_pt,2))).^2) < (0.05*r)^2;
0040     if any(ctr_pts)
0041         params(1,:) = params(1,:)/mean(params(1,ctr_pts));
0042     else
0043         eidors_msg('eval_GREIT_fig_merit: no centre points found to normalize',1);
0044     end
0045 end
0046 
0047 
0048 
0049 
0050 
0051 function ampl = calc_amplitude(img,qmi)
0052    if nargin == 1 
0053        ampl = sum(img(:));
0054    else 
0055        ampl = sum(img(:).*qmi(:));
0056    end
0057 
0058 function pe   = calc_posn_error(qmi, xmean, ymean, xy)
0059    % This definition allows + and - PE, but can also give zero in unexpected places
0060    pe = sqrt(sum(xy.^2)) - sqrt( xmean^2 + ymean^2);
0061    % This definition gives the absolute PE, but can't be negative
0062 %  pe = sqrt((xy(1,:) - xmean).^2 + ...
0063 %            (xy(2,:) - ymean).^2);
0064 
0065 function res  = calc_resolution(qmi, map)
0066    res = sqrt( sum(qmi(:)) / sum(map(:)));
0067 
0068 function sd  = calc_shape_deform(qmi, equiv_circ)
0069    not_circ= qmi & ~equiv_circ;
0070    sd = sum(not_circ(:))/sum(qmi(:));
0071 
0072 function rr = calc_ringing(img, qmi );
0073    ring_part =  img .* ( (img<0) & ~qmi);
0074    rr = -sum( ring_part(:) )/sum( img(:).*qmi(:) );
0075    
0076 function [x,y,bb_min,bb_max]=prepare_grid(sz,mdl)
0077       % bounding box
0078    bnd = unique(mdl.boundary);
0079    bb_min = min(mdl.nodes(bnd,:));
0080    bb_max = max(mdl.nodes(bnd,:));
0081    
0082    [x,y]=ndgrid(linspace(bb_min(1),bb_max(1),sz(1)),linspace(bb_min(2),bb_max(2),sz(2))); 
0083 
0084    
0085 function [xmean,ymean,equiv_circ,qmi,img] = calc_cofg(img,map,x,y);
0086 %  if abs(max(img(:))) < abs(min(img(:))); img= -img; end
0087    qmi = calc_hm_set( img, 0.25 );
0088    if sum(img(:) & qmi(:))<0 ; 
0089      error('problem in CofG calculation');
0090    end
0091    
0092    pix_sz = (max(x(:)) - min(x(:))) *( max(y(:)) - min(y(:))) /numel(img);
0093 
0094    %map = x.^2+y.^2<1.1;
0095    qmi = qmi.*map; img = img.*map;
0096 
0097 %  qmi = qmi .* img;  %USE THE IMAGE AMPLITUDE FOR THE CALCULATION
0098 
0099    ss_qmi = sum(qmi(:));
0100    xmean =  sum(sum( (qmi.*x) ))/ss_qmi; % centre of gravity
0101    ymean =  sum(sum( (qmi.*y) ))/ss_qmi;
0102    equiv_circ = (x-xmean).^2 + (y-ymean).^2 < pix_sz*ss_qmi/pi;