eval_GREIT_fig_merit

PURPOSE ^

EVAL_GREIT_FIG_MERIT: calculate GREIT figures of merit for images

SYNOPSIS ^

function params = eval_GREIT_fig_merit(imgs, xyzr_pt)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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.m 4808 2015-03-29 11:54:47Z bgrychtol-ipa $
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);
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    % This definition allows + and - PE, but can also give zero in unexpected places
0054    pe = sqrt(sum(xy.^2)) - sqrt( xmean^2 + ymean^2);
0055    % This definition gives the absolute PE, but can't be negative
0056 %  pe = sqrt((xy(1,:) - xmean).^2 + ...
0057 %            (xy(2,:) - ymean).^2);
0058 
0059 function res  = calc_resolution(qmi, map)
0060    res = sqrt( sum(qmi(:)) / sum(map(:)));
0061 
0062 function sd  = calc_shape_deform(qmi, equiv_circ)
0063    not_circ= qmi & ~equiv_circ;
0064    sd = sum(not_circ(:))/sum(qmi(:));
0065 
0066 function rr = calc_ringing(img, qmi );
0067    ring_part =  img .* ( (img<0) & ~qmi);
0068    rr = -sum( ring_part(:) )/sum( img(:).*qmi(:) );
0069    
0070 function [x,y,bb_min,bb_max]=prepare_grid(sz,mdl)
0071       % bounding box
0072    bnd = unique(mdl.boundary);
0073    bb_min = min(mdl.nodes(bnd,:));
0074    bb_max = max(mdl.nodes(bnd,:));
0075    
0076    [x,y]=ndgrid(linspace(bb_min(1),bb_max(1),sz(1)),linspace(bb_min(2),bb_max(2),sz(2))); 
0077 
0078    
0079 function [xmean,ymean,equiv_circ,qmi,img] = calc_cofg(img,map,x,y);
0080 %  if abs(max(img(:))) < abs(min(img(:))); img= -img; end
0081    qmi = calc_hm_set( img, 0.25 );
0082    if sum(img(:) & qmi(:))<0 ; 
0083      error('problem in CofG calculation');
0084    end
0085    
0086    pix_sz = (max(x(:)) - min(x(:))) *( max(y(:)) - min(y(:))) /numel(img);
0087 
0088    %map = x.^2+y.^2<1.1;
0089    qmi = qmi.*map; img = img.*map;
0090 
0091 %  qmi = qmi .* img;  %USE THE IMAGE AMPLITUDE FOR THE CALCULATION
0092 
0093    ss_qmi = sum(qmi(:));
0094    xmean =  sum(sum( (qmi.*x) ))/ss_qmi; % centre of gravity
0095    ymean =  sum(sum( (qmi.*y) ))/ss_qmi;
0096    equiv_circ = (x-xmean).^2 + (y-ymean).^2 < pix_sz*ss_qmi/pi;

Generated on Wed 21-Jun-2017 09:29:07 by m2html © 2005