% NBDcircle.m                                       2004/6/12
%
% Make matfile which stores data on neighborhoods with radius^2 = r2-1
% If the matfile exists, then load it into the current directory
%
% Zlabel: Image of size lenv x lenh with labels 0, 1, 2, ..., gg.
%         Pixels with label 0 are ignored.
% r2:     square of radius
% fname:  Matfile name for saving neighborhood data, "r2" will be automatically added.

function out = NBDinfo(Zlabel, fname, r2);
global NBDi NBDj NBDsize ssi ssj len

NBDname = [fname, '_r2_', num2str(r2),'_', num2str(r2-1), '.mat']; out = 0;

if exist(NBDname) > 1
 load(NBDname); out = 1;
 fprintf('*** Stored file is loaded ***\n')
else
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %  Linear discrimination based on ring neighborhoods
 %  radius = sqrt(r2)
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 r1         = ceil(sqrt(r2)); [lenv, lenh] = size(Zlabel);
 [ssi, ssj] = find(Zlabel > 0.5);     len  = length( ssi );
 NBDi    = zeros((2*r1 + 1)^2, len);  NBDj = zeros((2*r1 + 1)^2, len);
 NBDsize = zeros(len, 1);
 for s = 1:len
  i = ssi(s); j = ssj(s); count1 = 0;
  for j1 = -r1:r1; jj1 = j + j1;        % ring neighborhoods with radius = sqrt(r2)
   if jj1 > 0.5 & jj1 < lenh + 0.5
    for i1 = -r1:r1; ii1 = i + i1;
     if ii1 > 0.5 & ii1 < lenv + 0.5
      if i1^2 + j1^2 == r2
       if Zlabel(ii1, jj1) > 0.5
        count1 = count1 + 1;
        NBDi(count1, s) = ii1; NBDj(count1, s) = jj1;
       end % Zlabel
      end % if
     end % ii1
    end % i1
   end % jj1
  end % j1
  NBDsize(s) = count1;
 end % s

end % if

NBDmax = max(NBDsize); % NBDi = NBDi(1:NBDmax, :); NBDj = NBDj(1:NBDmax, :);
save(NBDname, 'NBDi', 'NBDj', 'NBDsize', 'r2', 'ssi', 'ssj', 'len')

fprintf('NBD info with radius^2 = %d saved as file %s\n', r2, NBDname)
fprintf('*** (min, max) of NBD size = ( %d, %d)', min(NBDsize),max(NBDsize))
fprintf(',  average of NBD size = %f  ***\n', mean(NBDsize))
fprintf('*** No of pixels with empty NBD = %d ***\n',length(find(NBDsize==0)))
