function [A,R] = co_occurrence(I, d,phi)
% 
% Written by:
% -- 
% John L. Weatherwax                2005-08-04
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
% Epage 
% 
%-----

assert( d>0, 'negative d' );
[di,dj] = phi_to_xy_increment(phi);

% Make sure that the pixels of I have values that are indexed from 1,2,...N_g:
uElts = unique( I(:) );
N_g   = length(uElts);
Icp   = zeros(size(I));
for ii=1:N_g
  inds = find( I(:)==uElts(ii) );
  Icp(inds) = ii;
end
I = Icp;

A = zeros(N_g,N_g);
[M,N] = size(I);

for ii=1:M % row index 
  for jj=1:N % column index 
    
    I_1 = I(ii,jj);
    
    % step in the positive d direction: 
    ii_new = ii+di*d;
    in_I_Range = (ii_new > 0) && (ii_new < M+1);
    jj_new = jj+dj*d;
    in_J_Range = (jj_new > 0) && (jj_new < N+1);
    
    if( in_I_Range && in_J_Range )
      I_2 = I(ii_new,jj_new);
      A(I_1,I_2) = A(I_1,I_2) + 1;
    end
    
    % step in the negative d direction: 
    ii_new = ii-di*d;
    in_I_Range = (ii_new > 0) && (ii_new < M+1);
    jj_new = jj-dj*d;
    in_J_Range = (jj_new > 0) && (jj_new < N+1);
    
    if( in_I_Range && in_J_Range )
      I_2 = I(ii_new,jj_new);
      A(I_1,I_2) = A(I_1,I_2) + 1;
    end
    
  end
end

R = sum(A(:));

A = A/R;