【问题标题】:Euclidean distance between images图像之间的欧几里得距离
【发布时间】:2023-03-12 22:16:01
【问题描述】:

我有两个图像,例如PS,大小为 8192×200,我想计算它们之间的自定义“欧几里得距离”。目前我使用以下步骤:

  1. 将图像重塑为一对列和行向量:

    Ip = Ip(:).';
    Is = Is(:);
    
  2. 计算一个度量矩阵G,其条目由公式给出

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

    其中r 是一个从0 到20 变化的全局参数,例如,d 是像素i 和像素j 之间的距离。例如,如果像素i(k,l) 并且像素j(k1,l1),那么d = sqrt((k-k1)^2 + (l-l1)^2);。像素 1 将是 (1,1),像素 2 将是 (1,2),依此类推。因此,矩阵G 的大小将是1638400×1638400

  3. 计算两个图像之间的最终(标量)欧几里得距离,使用:

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

我已经使用 mex 函数编写了一些代码,但是在给出结果之前需要很长时间(5-6 小时) - 请参阅 this SO question 获取代码和更多关于此的讨论。

请帮我优化一下;理想情况下,我希望它在几秒钟内运行。请注意,我对涉及 GPU 的解决方案不感兴趣。

【问题讨论】:

  • 我还是觉得你应该给个图..画个简单点就行了
  • 感谢您的编辑,但在第一步中,我将两个图像同时重塑为行向量或列向量,而不是一对行向量和列向量。我可以知道你背后的想法吗?
  • 如果您还需要图片,请告诉我。我可以放一个简单的。
  • 图像会很好,是的。鉴于您正在尝试做的事情,我认为您的第一步不是正确的方法。我的编辑确实稍微改变了它(您在转换为矢量之前翻转了两个图像)但这并不是特别有趣/与解决整个问题无关。我不确定您在评论中所说的“同时”是什么意思。
  • 同时,我的意思是IpIs 都是行向量或列向量。排除Ip为行向量,Is为列向量的情况,反之亦然。

标签: image algorithm matlab


【解决方案1】:

如果我理解正确的话,你应该可以做到以下几点,并让它在 2s 内运行:

样本数据:

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。见thisthis

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)运行速度会更快。

【讨论】:

  • 非常感谢您提供如此出色的解决方案。如果出现任何疑问,我将在此线程中再次与您联系。再次感谢。
  • 如果您可以检查此解决方案是否有效(或者可能几乎有效),那么请将其标记为正确。此外,编辑您之前的问题以链接到这个问题可能会很好。
  • 您能否解释一下步骤kernel = bsxfun(@plus, (-s1:s1).^2.', (-s2:s2).^2); 背后的想法。此外,在这一步之后,数字非常大,以至于 MATLAB 为表达式 exp(-kernel/ (2*r*2)) 生成了 0。因此,我认为这不符合我的目的。
  • @nagarwal - 第一行产生(一个正确的)大二维数组,给出中心点和每隔一个点之间的平方距离。卷积将内核移动到每个可能的位置,并为该位置相乘和求和,这正是您所需要的。关于“Matlab 正在产生 [几乎] 0”,我想您会发现这正是您想要的,我在答案的末尾记下了这一点。
  • 我对您的代码进行了彻底的分析,发现您的想法非常好。内核实际上是以多个 8152 x 200 阵列的形式存储每个像素与所有其他像素的距离,然后卷积完成其余的乘法和加法工作。我只觉得在第一行中,内核应该定义为kernel = bsxfun(@plus, (-s1+1:s1-1).^2.', (-s2+1:s2-1).^2); ,因为我认为这些角落的行和列没有任何距离。虽然我仍然相信,由于我们在conv2 期间使用了标志same,因此我们的结果不会受到影响。