Autocorrelation and Self Convolution

inputautocorrelation self convolution

Files used in this study may be downloaded from: autoself.zip

Simple example

input image: single.tif, method: auto.m

» a = imread('single.tif');
» out = auto(a);

The autocorrelation of the cyl(r) function shown above is given in the following routine:

function y = ideal(r)
%
% IDEAL
% y = ideal(r);
%
% calculates the autocorrelation of the cyl(r) function

 r = abs(r);
y = zeros(size(r));
i = find(r<1.0);
y(i) = 2.0*(acos(r(i))- r(i).*sqrt(1.0-r(i).*r(i)))/pi;

A cross-section of the autocorrelation function is shown below

The marked data points are taken from a horizontal cross-section of the output image. The solid curve is the ideal function, which was generated from ideal(r/8). The factor of 8 converts pixel distance to the appropriate normalized distance. 

» x = 0:8;
» y = out(65,x+65);
» plot(x,y);
» xi = linspace(0,8,100);
» yi = ideal(xi/8);
» plot(x,y,'o',xi,yi);
» xlabel('radius (pixels)');
» ylabel('value');

The agreement between theory and calculation is not as good as we might wish. To test the hypothesis that the disagreement is due to sampling error, we increased the size of the circle by four so that its autocorrelation just fits in a 64 x 64 array. The results, shown below, are much more satisfactory.

» a = imread('circ.tif');
» out = auto(a);
» x = 0:31;
» y = out(33,33+x);
» xi = linspace(0,31,100);
» yi = ideal(xi/32);
» plot(x,y,'o',xi,yi);
» xlabel('radius (pixels)');
» ylabel('value');

Second Example

The input function consists of two circles (cylinder functions), one displaced 8 pixels to the left and the other 8 pixels to the right. Each circle is 8 pixels in diameter. A cross-section of the autocorrelation function is shown below

There are three peaks. The central peak is twice the height of its neighbors. The relative peak heights and locations are shown in the diagram (b) below.

(a)(b)

Third Example

If the input function is real and has inversion symmetry (f(-x,-y) = f(x,y)) then the autocorrelation function is the same as the self convolution function. This example shows an input function without inversion symmetry.
input auto.m self.m

The patterns above are simple enough to be reasoned mentally. However, the Matlab function peak.m can be used to search for peaks in an image. The results of applying this program to the autocorrelation image are:

» a = imread('triple.tif');
» b = auto(a);
» peak(b,0.05,8);

image size: 64 x 64
block size: 8
threshold: 0.05
max # peaks: 20

peak	   value	 row	 col
   1	1.000000	  33	  33
   2	0.333333	  17	  17
   3	0.333333	  33	  17
   4	0.333333	  17	  33
   5	0.333333	  49	  33
   6	0.333333	  33	  49
   7	0.333333	  49	  49

Doing an autocorrelation

The autocorrelation of the input image, input, can be constructed from the command

» out = auto(input);

where auto is the following Matlab file

function f = auto(a)
%
% AUTO
% f = auto(a)
% returns as f the auto-correlation of image a.
f = fft2(a);
f = f.*conj(f);
f = ifft2(f);
f = fftshift(f);
f = real(f);
f = f/max(max(f));

Doing a self convolution

The self convolution of the input image, input, can be constructed from the command

» out = self(input);

where self is the following Matlab file:

function f = self(a)
%
% SELF
% f = self(a)
% returns as f the self-convolution of image a.

f = fft2(a);
f = f.*f;
f = ifft2(f);
f = fftshift(f);
f = real(f);
f = f/max(max(f));

Maintained by John Loomis, last updated 8 March 2000