calc_posn_resolution

 CALC_POSN_RESOLUTION: calculate posn and resolution of a target
 [br,xpos,ypos]= calc_posn_resolution( img, fwd_model)
 
 img - eidors img generated by a point (or small) target
 br  - Blur radius measure of resolution
 pos - [x;y] position in terms of % radius

 These parameters are defined in Adler & Graudo (1996),
  IEEE T. Medical Imaging, 5:170-179, 1996.

 br and pos will have columns for each image in img

 In 3D calc_posn_resolution does a 2D z-slice in the 
 centre of the medium. This may not be appropriate

0001 function [br,pos]= calc_posn_resolution( img, fwd_model)
0002 % CALC_POSN_RESOLUTION: calculate posn and resolution of a target
0003 % [br,xpos,ypos]= calc_posn_resolution( img, fwd_model)
0004 %
0005 % img - eidors img generated by a point (or small) target
0006 % br  - Blur radius measure of resolution
0007 % pos - [x;y] position in terms of % radius
0008 %
0009 % These parameters are defined in Adler & Graudo (1996),
0010 %  IEEE T. Medical Imaging, 5:170-179, 1996.
0011 %
0012 % br and pos will have columns for each image in img
0013 %
0014 % In 3D calc_posn_resolution does a 2D z-slice in the
0015 % centre of the medium. This may not be appropriate
0016 
0017 % (C) 2006 Andy Adler. Licensed under GPL version 2
0018 % $Id: calc_posn_resolution.m 3053 2012-06-06 15:55:32Z aadler $
0019 
0020 if nargin>=2
0021    img.fwd_model= fwd_model;
0022 end
0023 
0024 % create full 'red' image
0025 fimg= img;
0026 fimg.elem_data= ones(size(img.elem_data,1),1);
0027 
0028 np= 256;
0029 % FIXME: if the functions have an error, then npoints is changed
0030 np_save= calc_colours('npoints');
0031 calc_colours('npoints',np);
0032 
0033 if     size(img.fwd_model.elems,2) == 3 % 2D
0034    rimg= calc_slices( img );
0035    fimg= calc_slices( fimg );
0036 elseif size(img.fwd_model.elems,2) == 4 % 3D
0037    rimg= calc_slices( img, 1 );
0038    fimg= calc_slices( fimg, 1 );
0039 else
0040    calc_colours('npoints',np_save);
0041    error('N dimensions is not 2 or 3');
0042 end
0043 
0044 
0045 rimg= squeeze(rimg); % for 3D case
0046 n_imgs= size(rimg,3);
0047 % reset npoints
0048 calc_colours('npoints',np_save);
0049 
0050 % Process full image to get area and x,y limits
0051 fimg= ~isnan(fimg);
0052 img_area= sum(fimg(:));
0053 
0054 xfind= find(any( fimg,1));
0055 img_xmin= min(xfind); img_xmax= max(xfind);
0056 img_xctr= mean([img_xmin, img_xmax]);
0057 img_xrad= [img_xmax - img_xmin]/2; % radius
0058 
0059 yfind= find(any( fimg,1));
0060 img_ymin= min(yfind); img_ymax= max(yfind);
0061 img_yctr= mean([img_ymin, img_ymax]);
0062 img_yrad= [img_ymax - img_ymin]/2; % radius
0063 
0064 br= zeros(1,n_imgs);
0065 pos=zeros(2,n_imgs);
0066 for i=1:n_imgs
0067    img_i =       rimg(:,:,i);
0068    % flip so image is +
0069    img_i =       img_i * sign( max(img_i(:)) + min(img_i(:)));
0070    [img_i,idx] = sort(img_i(:));
0071    img_i(idx)  = cumsum(img_i .* (img_i>0));
0072    ha_set = img_i > max(img_i)/2;
0073    br(i) = sqrt(sum(ha_set)/img_area);
0074 
0075 % create HA image
0076    ha_img= zeros(np,np);
0077    ha_img( ha_set) = 1; 
0078    xpos = sum( ha_img, 1);
0079    xpos = sum( (1:np) .* xpos ) / sum( xpos );
0080    xpos = (xpos - img_xctr) / img_xrad;
0081    ypos = sum( ha_img, 2)';
0082    ypos = sum( (1:np) .* ypos ) / sum( ypos );
0083    ypos = -(ypos - img_yctr) / img_yrad; % images y dirn is backward
0084    pos(:,i) = [xpos;ypos];
0085 end
0086    
0087    
0088