function cee = cv_auto_covariance(image_in)
%
% cv_auto_covariance - compute covariance of image
% Call: cee = cv_auto_covariance(image_in);
% On input:
%     image_in: input image
% On output:
%     cee: covariance of image: i,j is covariance offset  
% Author:
%     Tom Henderson
%     7 January 2000
% Custom Functions Used:
%     none
% Method:
%     see p. 33 of Trucco and Verri, Intro Techniques for 3D Computer Vision
% Testing:
% > C = randn(5,5)
%
% C =
%
%    0.8580   -0.3999    0.6686   -1.6041    0.5287
%    1.2540    0.6900    1.1908    0.2573    0.2193
%   -1.5937    0.8156   -1.2025   -1.0565   -0.9219
%   -1.4410    0.7119   -0.0198    1.4151   -2.1707
%    0.5711    1.2902   -0.1567   -0.8051   -0.0592
%
% > Cee = cv_auto_covariance(C)
%
% Cee =
%
%    1.0529   -0.1555    0.2251   -0.1737    0.2077
%    0.0150   -0.1565    0.0452   -0.1924    0.0974
%   -0.0669    0.1094   -0.1314    0.0363   -0.1406
%   -0.1449    0.2286   -0.1101    0.0427   -0.0775
%    0.0468    0.0290    0.0063   -0.0272   -0.0007
%

[num_rows, num_cols] = size(image_in);
N = min(num_rows,num_cols);   % make sure image is square
image_mean = mean2(image_in);
c = 1/(N^2);
cee = zeros(N,N);
for step_row_index = 0:N-1          % row offset loop
   max_row = N - step_row_index;    
   for step_col_index = 0:N-1       % col offset loop
      max_col = N - step_col_index;
      temp_sum = 0;
      for row = 1:max_row
         for col = 1:max_col
            temp_sum = temp_sum + (image_in(row,col)-image_mean)* ...
               (image_in(row+step_row_index,col+step_col_index)-image_mean);
         end
      end
      cee(step_row_index+1,step_col_index+1) = c*temp_sum;
   end
end
