see http://stackoverflow.com/q/36977008/2399799

so_answer.m

% 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'); % see below for faster, fft based convolution
res = g(:)' * img_diff(:); 


%%
function c = conv2fft(X, Y)
    % ignoring small floating-point differences, this is equivalent
    % to the inbuilt Matlab conv2(X, Y, 'same').
    % https://en.wikipedia.org/wiki/Convolution_theorem
    % http://stackoverflow.com/a/21757317/2399799
    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