0001 function params = eval_GREIT_fig_merit(imgs, xyzr_pt)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016 mdl = imgs.fwd_model;
0017 imgs = calc_slices(imgs);
0018 map = ~isnan(squeeze(imgs(:,:,1)));
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
0034 ctr = bb_min + 0.5*(bb_max-bb_min);
0035 r = max(0.5*(bb_max-bb_min));
0036 if N_imgs > 10
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
0054 pe = sqrt(sum(xy.^2)) - sqrt( xmean^2 + ymean^2);
0055
0056
0057
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
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
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
0089 qmi = qmi.*map; img = img.*map;
0090
0091
0092
0093 ss_qmi = sum(qmi(:));
0094 xmean = sum(sum( (qmi.*x) ))/ss_qmi;
0095 ymean = sum(sum( (qmi.*y) ))/ss_qmi;
0096 equiv_circ = (x-xmean).^2 + (y-ymean).^2 < pix_sz*ss_qmi/pi;