Euclidean distance between images

Question

I have two images, say P and S, of size 8192×200, and I want to calculate a custom "Euclidean distance" between them. Currently I use the following steps:

1. Reshape the images into a pair of column and row vectors:

   Ip = Ip(:).';
   Is = Is(:);
   

2. Calculate a metric matrix, G, whose entries are given by the formula

   G(i,j) = 1/(2*pi*r*r) * exp((-d*d)/(2*r*r));
   

   where r is a global parameter that varies from 0 to 20, say, and d is the distance between pixel i and pixel j. E.g., if pixel i is (k,l) and pixel j is (k1,l1), then d = sqrt((k-k1)^2 + (l-l1)^2);. Pixel 1 will be (1,1), Pixel 2 will be (1,2), and so on. Therefore, the size of matrix G will be 1638400×1638400.

3. Compute the final (scalar) Euclidean distance between two images, using:

   ImEuDist = sqrt( (Ip-Is) * G * (Ip-Is).' );  
   

I have already written some code using a mex function, but it is taking too long before giving the results (5-6 Hours) - see this SO question for code and more discussion on this.

Please help me to optimize this; I would ideally like it to run in a matter of seconds. Note that I am not interested in solutions involving the GPU.

Answer 1

If I've understood correctly, you should be able to do the following, and have it run in under 2s:

sample data:

s1 = 8192; s2 = 200;
img_a = rand(s1, s2);
img_b = rand(s1, s2);
r = 2;

and the calculation itself:

img_diff = img_a - img_b;
kernel = bsxfun(@plus, (-s1:s1).^2.', (-s2:s2).^2);
kernel = 1/(2/pi/r^2) * exp(-kernel/ (2*r*2));
g = conv2(img_diff, kernel, 'same');
res = g(:)' * img_diff(:); 
res = sqrt(res);

The above takes about 25s. To get down to 2s, you need to replace the standard conv2 with a faster, fft based convolution. See this and this:

function c = conv2fft(X, Y)
    % ignoring small floating-point differences, this is equivalent
    % to the inbuilt Matlab conv2(X, Y, 'same')
    X1 = [X zeros(size(X,1), size(Y,2)-1);
          zeros(size(Y,1)-1, size(X,2)+size(Y,2)-1)];
    Y1 = zeros(size(X1)); 
    Y1(1:size(Y,1), 1:size(Y,2)) = Y;
    c = ifft2(fft2(X1).*fft2(Y1));
    c = c(size(X,1)+1:size(X,1)+size(X,1), size(X,2)+1:size(X,2)+size(X,2));
end

Incidentally, if you still want it to go faster, you could make use of the fact that exp(-d^2/r^2) gets very close to zero for fairly small d: so you can actually crop your kernel to just a tiny rectangle, rather than the huge thing suggested above. A smaller kernel means conv2fft (or especially conv2) will run faster.