calc_sim_params

 params = GREIT_sim_params(imgs)
  params(1,:) = Image Amplitude
  params(2,:) = Position Error => + toward centre, - toward edge
  params(3,:) = Resolution
  params(4,:) = Shape Deformation
  params(5,:) = Ringing

0001 function params = calc_sim_params(imgs, xyzr_pt)
0002 % params = GREIT_sim_params(imgs)
0003 %  params(1,:) = Image Amplitude
0004 %  params(2,:) = Position Error => + toward centre, - toward edge
0005 %  params(3,:) = Resolution
0006 %  params(4,:) = Shape Deformation
0007 %  params(5,:) = Ringing
0008 
0009 % (C) 2008 Andy Adler. Licensed under GPL v2 or v3
0010 % $Id: calc_sim_params.html 2819 2011-09-07 16:43:11Z aadler $
0011 
0012 N_imgs = size(imgs,3);
0013 for i= 1:N_imgs
0014    [xmean,ymean,equiv_circ,map,qmi,img] = calc_cofg(imgs(:,:,i));
0015    params(1,i) = calc_amplitude( img );
0016    params(2,i) = calc_posn_error( qmi, xmean, ymean, xyzr_pt(1:2,i) );
0017    params(3,i) = calc_resolution( qmi, map );
0018    params(4,i) = calc_shape_deform( qmi, equiv_circ );
0019    params(5,i) = calc_ringing( img, qmi );
0020 end
0021 
0022 % TODO: Fix this when we start to care about units
0023 if N_imgs > 20
0024     params(1,:) = params(1,:)/mean(params(1,1:10));
0025 end
0026 function ampl = calc_amplitude(img)
0027    ampl = sum(img(~isnan(img)));
0028 
0029 function pe   = calc_posn_error(qmi, xmean, ymean, xy)
0030    pe = sqrt(sum(xy.^2)) - sqrt( xmean^2 + ymean^2);
0031 
0032 function res  = calc_resolution(qmi, map)
0033    res = sqrt( sum(qmi(:)) / sum(map(:)));
0034 
0035 function sd  = calc_shape_deform(qmi, equiv_circ)
0036    not_circ= qmi & ~equiv_circ;
0037    sd = sum(not_circ(:))/sum(qmi(:));
0038 
0039 function rr = calc_ringing(img, qmi );
0040    img(isnan(img)) = 0;
0041    ring_part =  img .* ( (img<0) & ~qmi);
0042    rr = -sum( ring_part(:) )/sum( img(:).*qmi(:) );
0043 
0044 function [xmean,ymean,equiv_circ,map,qmi,img] = calc_cofg(img);
0045 %  if abs(max(img(:))) < abs(min(img(:))); img= -img; end
0046    qmi = calc_hm_set( img, 0.25 );
0047    [x,y]=meshgrid(linspace(-1,1,32),linspace(-1,1,32)); map = x.^2+y.^2<1.1;
0048    if sum(img(map) & qmi(map))<0 ; keyboard ; end
0049    qmi = qmi.*map; img = img.*map;
0050 
0051    ss_qmi = sum(qmi(:));
0052    xmean =  sum(sum( (qmi.*x) ))/ss_qmi; % centre of gravity
0053    ymean =  sum(sum( (qmi.*y) ))/ss_qmi;
0054    equiv_circ = (x-xmean).^2 + (y-ymean).^2 < ss_qmi/pi/(32/2)^2;