0

So when I try to find a template B in a big image A, I can do it by finding the maximum of cross-correlation, like this in spatial domain:

%       Finding maximum of correlation:
        phi = normxcorr2(B,A); 
        [ymax, xmax] = find(phi == max(phi(:)));

%       Find position in original image:
        ypeak = ymax - size(B,1);
        xpeak = xmax - size(B,2);

But when I want to do it in frequency domain, I get wrong results:

%       Calculate correlation in frequency domain:
        Af = fft2(A);
        Bf = fft2(B, size(A,1), size(A,2));
        phi2f = conj(Af)'*Bf;

%       Inverse fft to get back to spatial domain:       
        phi2 = real(ifft(fftshift(phi2f)));

%       Now we have correlation matrix, excatly the same as calculated in
%       the spatial domain.
        [ymax2, xmax2] = find(phi2 == max(phi2(:)));

I don't understand what I'm doing wrong in the frequency domain. I've tried it without fftshift, it gives a different result, albeit a wrong one still. How can I do this correctly?

Vidak
  • 1,083
  • 14
  • 29
  • Looking at the definition of `normcorr2`, I guess we can assume that B is your template and A is your image? I think it would be good if you added that, just for clarity. – kkuilla Oct 09 '14 at 08:12
  • And if you have the images, add those as well together with the expected and real output so that it becomes a repeatable example. – kkuilla Oct 09 '14 at 08:19
  • I think `conj(Af)'` takes the non-conjugate transposed of `Af`, is that really what you want? – McMa Oct 09 '14 at 08:23
  • @McMa it takes the conjugate transposed, just checked using (conj(Af))' – Vidak Oct 09 '14 at 08:57
  • I am afraid it doesn't: the transpose operator `'` already performs a conjugated transposition, by using `conj()` you are only canceling the initial conjugation, leaving you with the non-conjugated transpose. That's the first thing that strikes off to me. – McMa Oct 09 '14 at 09:17

1 Answers1

1

This should do the trick:

t   = imread('cameraman.tif');
a   = imtranslate(t, [15, 25]);

% Determine padding size in x and y dimension
size_t      = size(t);
size_a      = size(a);
outsize     = size_t + size_a - 1;

% Determine 2D cross correlation in Fourier domain
Ft = fft2(t, outsize(1), outsize(2));
Fa = fft2(a, outsize(1), outsize(2));
c = abs( fftshift( ifft2(Fa .* conj(Ft))) );

% Find peak
[max_c, imax]   = max(abs(c(:)));
[ypeak, xpeak]  = ind2sub(size(c), imax(1));

% Correct found peak location for image size
corr_offset = round([(ypeak-(size(c, 1)+1)/2) (xpeak-(size(c, 2)+1)/2)]);

% Write out offsets
y_offset = corr_offset(1)
x_offset = corr_offset(2)
Maurits
  • 2,082
  • 3
  • 28
  • 32
  • 1
    interesting! How would you change your code to calculate the Normalized Cross Correlation version? – giò Nov 30 '16 at 15:53