# Calculate Earth Mover's Distance for two grayscale images

I am trying to calculate EMD (a.k.a. Wasserstein Distance) for these two grayscale (299x299) images/heatmaps:

Right now, I am calculating the histogram/distribution of both images. The histograms will be a vector of size 256 in which the nth value indicates the percent of the pixels in the image with the given darkness level. Then, using these to histograms, I am calculating the EMD using the function wasserstein_distance from scipy.stats.

However, I am now comparing only the intensity of the images, but I also need to compare the location of the intensity of the images.

How can I do this?

I am thinking about obtaining a histogram for every row of the images (which results in 299 histograms per image) and then calculating the EMD 299 times and take the average of these EMD's to get a final score. Is this the right way to go?

Other methods to calculate the similarity bewteen two grayscale are also appreciated.

Having looked into it a little more than at my initial answer: it seems indeed that the original usage in computer vision, e.g. Peleg et al. (1989), simply matched between pixel values and totally ignored location. Later work, e.g. Rubner et al. (2000), did the same but on e.g. local texture features rather than the raw pixel values. This then leaves the question of how to incorporate location.

Doing it row-by-row as you've proposed is kind of weird: you're only allowing mass to match row-by-row, so if you e.g. slid an image up by one pixel you might have an extremely large distance (which wouldn't be the case if you slid it to the right by one pixel).

A more natural way to use EMD with locations, I think, is just to do it directly between the image grayscale values, including the locations, so that it measures how much pixel "light" you need to move between the two. This is then a 2-dimensional EMD, which scipy.stats.wasserstein_distance can't compute, but e.g. the POT package can with ot.lp.emd2.

Doing this with POT, though, seems to require creating a matrix of the cost of moving any one pixel from image 1 to any pixel of image 2. Since your images each have $$299 \cdot 299 = 89,401$$ pixels, this would require making an $$89,401 \times 89,401$$ matrix, which will not be reasonable.

One option would be to simply compute the distance on images which have been resized smaller (by simply adding grayscales together). If you downscaled by a factor of 10 to make your images $$30 \times 30$$, you'd have a pretty reasonably sized optimization problem, and in this case the images would still look pretty different. However, this is naturally only going to compare images at a "broad" scale and ignore smaller-scale differences.

There are also, of course, computationally cheaper methods to compare the original images. You can think of the method I've listed here as treating the two images as distributions of "light" over $$\{1, \dots, 299\} \times \{1, \dots, 299\}$$ and then computing the Wasserstein distance between those distributions; one could instead compute the total variation distance by simply $$\operatorname{TV}(P, Q) = \frac12 \sum_{i=1}^{299} \sum_{j=1}^{299} \lvert P_{ij} - Q_{ij} \rvert,$$ or similarly a KL divergence or other $$f$$-divergences. These are trivial to compute in this setting but treat each pixel totally separately. There are also "in-between" distances; for example, you could apply a Gaussian blur to the two images before computing similarities, which would correspond to estimating $$L_2(p, q) = \int (p(x) - q(x))^2 \mathrm{d}x$$ between the two densities with a kernel density estimate. What distance is best is going to depend on your data and what you're using it for.

• If I understand you correctly, I have to do the following: Suppose I have two 2x2 images. Then we have: C1=[0, 1, 1, sqrt(2)], C2=[1, 0, sqrt(2), 1], C3=[1, \sqrt(2), 0, 1], C4=[\sqrt(2), 1, 1, 0] The cost matrix is then: C=[C1, C2, C3, C4]. If I need to do this for the images shown above, I need to provide 299x299 cost matrices?! Sounds like a very cumbersome process. Or is there something I do not understand correctly? – LVDW Apr 26 at 14:34
• Yeah, I think you have to make a cost matrix of shape [299*299, 299*299]. It should only be a few lines of code though, e.g. by making the matrix of indices with something like np.indices and then using ot.utils.dist. (This will be computationally wasteful, but for input size 299 it shouldn't really matter.) – Dougal Apr 26 at 19:02
• You said I need a cost matrix for each image location to each other location. Doesnt this mean I need 299*299=89401 cost matrices? – LVDW Apr 29 at 9:30
• @LVDW I updated the answer; you only need one matrix, but it's really big, so it's actually not really reasonable. – Dougal Apr 29 at 13:16
• I think that would be not ridiculous, but it has a slightly weird effect of making the distance very much not invariant to rotating the images 45 degrees. Whether this matters or not depends on what you're trying to do with it. – Dougal Apr 29 at 15:42