function [h,p,r] = hist2(y1,y2,x1,x2); 
% y1 and y2 are vectors with data (must have same dimension)
% x1 and x2 are centers of histogram bins for the 2 coordinates. 
%    should be equally spaced !!!
% h is the 2-D histogram (y1 distrib in cols, y2 in rows)
% 
% p is the estimate of 2-D prob. dens. function
%
% r is he autocorrelation coefficient computed using the theoretical 
% formula (5-8 in Sebesta) 
%
% ... histogram computation is not very optimized ... 

L1 = length(x1); L2 = length(x2); 
N = length(y1); 

% alloc for hist 
h = zeros(L1,L2); out = 0; % counter of bad samples...
for ii=1:N,   % cycle on x1 data
  h1 = hist(y1(ii),x1); % will produce a hist with just one '1' in the output
  h2 = hist(y2(ii),x2); % will produce a hist with just one '1' in the output
  % find indices ... 
  ih1 = find (h1);       % finds the index of '1', if empty, we were out of 
                       % limits ...
  ih2 = find (h2);       % finds the index of '1', if empty, we were out of 
                       % limits ...
  % assign to hist if both non-empty
  if (~isempty(ih1) & ~isempty(ih1))
    h(ih1,ih2) =  h(ih1,ih2) +1; 
  else 
    out=out+1;
  end
end
%disp(['hist2: ' num2str(out) ' out of ' num2str(N) ' were out of bins']); 
% hm hm ... probably can't happen, it looks like the edge bins are going to 
% collect it ... 

%%% prob dens: will have to normalize
surf = (x1(2) - x1(1)) * (x2(2) - x2(1));  % surface of one tile
p = h / N / surf;  

%%%% autocor coeff 
% make col vector out of x1 and clone it L2 times. 
x1 = x1(:); X1 = repmat(x1,1,L2);
% make row vector out of x2 and clone it L1 times. 
x2 = x2(:)'; X2 = repmat(x2,L1,1); 
% now put it together, don't forget to multipl by tile surface
r = sum(sum (X1 .* X2 .* p)) * surf;

%%% check ... 
check = sum(sum (p)) * surf; 
disp(['hist2: check-should be 1:' num2str(check)]);