图像之间的欧几里得距离

5

enter image description here我有两张图片,假设分别为PS,大小均为8192×200,我想计算它们之间的自定义“欧几里得距离”。目前我使用以下步骤:

  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).' );  
    
我已经使用了mex函数编写了一些代码,但是在给出结果之前需要花费太长时间(5-6小时)-请参见这个SO问题,其中包含了代码和更多讨论。
请帮助我优化这个程序;我希望它能在几秒钟内运行。请注意,我不对涉及GPU的解决方案感兴趣。

我仍然认为你应该提供一个图表..只需在画图中制作一些简单的东西。 - dan-man
感谢您的编辑,但在第一步中,我同时将两个图像重塑为行向量或列向量,而不是成对的行向量和列向量。请问您这样做的想法是什么? - nagarwal
如果你还需要一张图片,请告诉我。我可以放一张简单的。 - nagarwal
同时,我的意思是IpIs都是行向量或列向量。这不包括Ip是行向量而Is是列向量以及反之的情况。 - nagarwal
不,你原本写的一个是行,一个是列。但我认为这并不重要,因为解决问题与此无关。 - dan-man
显示剩余2条评论
1个回答

3
如果我理解正确,您应该能够做到以下操作,并且在不到2秒的时间内运行它:
示例数据:
s1 = 8192; s2 = 200;
img_a = rand(s1, s2);
img_b = rand(s1, s2);
r = 2;

以及计算本身:

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);

上述过程大约需要25秒。要将其缩短到2秒,您需要用一种更快的基于fft的卷积替换标准的conv2。请参见此处此处

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

顺便提一下,如果你还想让它跑得更快,你可以利用这个事实:exp(-d^2/r^2) 对于相当小的 d 得到一个非常接近于零的值:因此,你实际上可以将你的卷积核裁剪成一个小矩形,而不是上面建议的巨大的东西。较小的卷积核意味着 conv2fft(或尤其是conv2)会运行得更快。


非常感谢您提供的出色解决方案。如果有任何疑问,我将在此线程中再次与您联系。再次感谢。 - nagarwal
如果您可以确认此解决方案有效(或者可能有效),请将其标记为正确答案。另外,将此问题链接到您之前的问题可能会更好。 - dan-man
能否请您解释一下 kernel = bsxfun(@plus, (-s1:s1).^2.', (-s2:s2).^2); 这个步骤的想法?此外,在这一步之后,数字非常大,MATLAB 对表达式 exp(-kernel/ (2*r*2)) 的结果为0。因此,我认为这不符合我的目的。 - nagarwal
@nagarwal - 第一行生成了一个正确的大型二维数组,给出了中心点与每个其他点之间的平方距离。卷积将核移动到每个可能的位置,并为该位置乘法和求和,这正是您所需要的。关于“Matlab产生[几乎] 0”,我认为您会发现这正是您想要的,我在答案的结尾做了记录。 - dan-man
我对你的代码进行了彻底的分析,发现你的想法真的很好。实际上,内核将每个像素与所有其他像素之间的距离存储为多个8152 x 200数组的形式,然后卷积完成乘法和加法的其余工作。我只觉得在第一行中,内核应该被定义为 kernel = bsxfun(@plus, (-s1+1:s1-1).^2.', (-s2+1:s2-1).^2); 因为我不认为这些角落的行和列会产生任何距离。虽然我仍然相信,由于我们在 conv2期间使用了标志 same,因此我们的结果不会受到影响。 - nagarwal
@dan-man 你好,希望你一切都好。最近,我想从使用高斯核函数转换为双曲正切(Sigmoid)核函数,其表达式在此处定义[http://crsouza.com/2010/03/kernel-functions-for-machine-learning-applications/]。根据这个表达式,我想重新定义上述问题中的矩阵G。你能否请指点一下如何实现。先谢谢了。 - nagarwal

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接