与我最初的答案相比，对其进行了更多的研究：似乎确实是计算机视觉中的原始用法，例如Peleg 等人。(1989)，在像素值和完全忽略的位置之间简单匹配。后来的工作，例如Rubner 等人。(2000)，做了同样的事情，但在例如局部纹理特征而不是原始像素值上。这就留下了如何合并位置的问题。

按照您的建议逐行执行有点奇怪：您只允许质量逐行匹配，因此如果您将图像向上滑动一个像素，您可能会有非常大的距离（如果您将其向右滑动一个像素，则不会出现这种情况）。

我认为，将 EMD 与位置一起使用的一种更自然的方法是直接在图像灰度值（包括位置）之间进行，以便测量您需要在两者之间移动多少像素“光”。这是一个二维 EMD，它scipy.stats.wasserstein_distance无法计算，但例如POT 包可以使用ot.lp.emd2.

但是，使用 POT 执行此操作似乎需要创建一个将任何一个像素从图像 1 移动到图像 2 的任何像素的成本矩阵。由于您的图像每个都有像素，因此需要制作矩阵，这是不合理的。$299 \cdot 299 = 89,401$$89,401 \times 89,401$

---

更新：可能比我下面描述的更好的方法是使用切片的 Wasserstein 距离，而不是普通的 Wasserstein。这利用了 1 维 Wasserstein 的计算效率非常高的事实，并通过取数据的随机一维投影之间的 Wasserstein 距离的平均值来$d$

这类似于您进行行和列传输的想法：对应于两个特定的投影。但是通过对投影进行平均，您可以获得一个真实的距离，这也比完整的 Wasserstein 具有更好的样本复杂性。

在（未经测试的、低效的）Python 代码中，可能如下所示：

import numpy as np
from scipy.stats import wasserstein_distance

def sliced_wasserstein(X, Y, num_proj):
    dim = X.shape[1]
    ests = []
    for _ in range(num_proj):
        # sample uniformly from the unit sphere
        dir = np.random.rand(dim)
        dir /= np.linalg.norm(dir)

        # project the data
        X_proj = X @ dir
        Y_proj = Y @ dir

        # compute 1d wasserstein
        ests.append(wasserstein_distance(X_proj, Y_proj)
    return np.mean(ests)

（这里的循环，至少直到得到X_projand Y_proj，可以被矢量化，这可能会更快。）

---

另一种选择是简单地计算已缩小尺寸的图像上的距离（通过简单地将灰度相加）。如果您将图像缩小 10 倍以使图像成为，那么您将遇到一个相当合理的优化问题，在这种情况下，图像看起来仍然非常不同。但是，这自然只会比较“大”尺度的图像而忽略较小尺度的差异。$30 \times 30$

当然，还有一些计算成本更低的方法来比较原始图像。您可以将我在这里列出的方法视为将两个图像视为“光”在上的分布，然后计算Wasserstein这些分布之间的距离；可以通过简单的 \operatorname{TV}(P, Q) = \ frac12或 类似的KL 散度或其他$\{1, \dots, 299\} \times \{1, \dots, 299\}$

$$\operatorname{TV}(P, Q) = \frac12 \sum_{i=1}^{299} \sum_{j=1}^{299} \lvert P_{ij} - Q_{ij} \rvert,$$ $f$. 在此设置中计算这些是微不足道的，但将每个像素完全分开处理。还有“中间”距离；例如，您可以在计算相似度之前对两个图像应用高斯模糊，这对应于估计 两个密度之间的核密度估计。最佳距离取决于您的数据以及您使用它的目的。 $$L_2(p, q) = \int (p(x) - q(x))^2 \mathrm{d}x$$

我认为对于您的图像大小要求，@Dougal 建议的切片 wasserstein 可能是最适合的，因为 299^4 * 4 字节意味着32 GBs传输矩阵的内存需求为 ~ ，这是非常巨大的。

为了完成回答使用 EMD 比较两个灰度图像的一般问题，如果估计速度是一个标准，还可以考虑POT通过命令工具箱中可用的正则化 OT 距离ot.sinkhorn(a, b, M1, reg)：正则化版本应该优化为比ot.emd(a, b, M1)命令更快的解决方案。