function [asm,con,cor,ent]=cooc(x,d,N_block)
% COOC	Compute the four major co-occurrence measures of an image
%	[asm,con,cor,ent]=cooc(x,[d_v d_h],[N_1 N_2]) computes 
%	ASM, CON, COR, and ENT measures, with displacement d_v pixels
%	vertically and d_h pixels horizontally, from N_1 by N_2 
%	subimages. The input image is assumed to be integer valued, and 
%	most meaningful results may be expected it it's in the range from
%	0 to somewhere between 7 and 31. If the maximum is above 63,
%	the image is rescaled to 0 to 7.

[M,N]=size(x);
if length(N_block==1)
  N_block=[N_block N_block];
end

% Build up a list of gray levels in the image
AllSamples=x;
GrayLevels=[];
while ~isempty(AllSamples)
  GrayLevels=[GrayLevels AllSamples(1)];
  AllSamples=AllSamples(find(AllSamples~=AllSamples(1)));
end
GrayLevels=sort(GrayLevels);
N_gray=length(GrayLevels);
if N_gray>16
  fprintf(2,'\07WARNING: more than 16 gray levels. Inefficient!\n');
end

% Displacement
x_d=vshift(hshift(x,d(2)),d(1));

% First pass: estimate mu_x and mu_y
mu_x=0; mu_y=0;
for i=GrayLevels % 			Gray level of base image
  fprintf(1,'''');
  for j=GrayLevels % 			Gray level of displaced image
    match = (x==i) & (x_d==j);
    p = avgsmth(match,N_block(1),N_block(2))/N_gray;
    mu_x=mu_x + i*p; % 			Estimate mu_x
    mu_y=mu_y + j*p; % 			Estimate mu_y
  end
end

% Second pass: estimate sigma_x and sigma_y
sigma_x=0; sigma_y=0;
for i=GrayLevels % 			Gray level of base image
  fprintf(1,'\b`\b');
  for j=GrayLevels % 			Gray level of displaced image
    match = (x==i) & (x_d==j);
    p = avgsmth(match,N_block(1),N_block(2))/N_gray;
    sigma_x=sigma_x + (i-mu_x).^2.*p; % Estimate sigma_x
    sigma_y=sigma_y + (j-mu_y).^2.*p; % Estimate sigma_y
  end
end 
sigma_x=sqrt(sigma_x);
sigma_y=sqrt(sigma_y);

% Third pass: estimate measures
asm=zeros(M,N); con=zeros(M,N); cor=zeros(M,N); ent=zeros(M,N);
for i=GrayLevels % 			Gray level of base image
  fprintf(1,'.');
  for j=GrayLevels % 			Gray level of displaced image
    match = (x==i) & (x_d==j);
    p = avgsmth(match,N_block(1),N_block(2));
    % Angular second moment
    asm=asm + p.^2; 
    % Contrast
    con=con + (i-j)^2 * p; 
    % Correlation
    cor=cor + (i*j*p - mu_x.*mu_y)./(sigma_x.*sigma_y); 
    % Entropy
    ent=ent - p .* log(p+(p==0)); 
  end
end 



