使用NumPy进行灰度图像的直方图均衡化

Question

如何轻松对存储在NumPy数组中的多个灰度图像进行直方图均衡化？

我有96x96像素的NumPy数据，格式为4D：

(1800, 1, 96,96)

- pbu

---

这个有用吗：https://web.archive.org/web/20151219221513/http://www.janeriksolem.net/2009/06/histogram-equalization-with-python-and.html？ - Martin Thoma

谢谢你，moose。当然，这对PIL有帮助。是否有任何纯numpy的解决方案？因为我要处理1万张图片，所以速度可能更快。PIL和Skimage通常需要很长时间。 - pbu

@pbu 在这个例子中，只有在读取图像数据时才需要PIL。由于您已经拥有了数据，因此您不需要它，只需要NumPy即可。请查看我在答案中编写的代码。同时非常感谢moose提供的解决方案链接！ - Trilarion

Answer 1

Moose的comment指向这个blog entry非常好。为了完整起见，我在此提供一个示例，使用更好的变量名，并对1000个96x96图像进行循环执行，这些图像在一个4D数组中，与问题中一样。它很快（在我的电脑上只需1-2秒钟），并且只需要NumPy。

import numpy as np

def image_histogram_equalization(image, number_bins=256):
    # from http://www.janeriksolem.net/histogram-equalization-with-python-and.html

    # get image histogram
    image_histogram, bins = np.histogram(image.flatten(), number_bins, density=True)
    cdf = image_histogram.cumsum() # cumulative distribution function
    cdf = (number_bins-1) * cdf / cdf[-1] # normalize

    # use linear interpolation of cdf to find new pixel values
    image_equalized = np.interp(image.flatten(), bins[:-1], cdf)

    return image_equalized.reshape(image.shape), cdf

if __name__ == '__main__':

    # generate some test data with shape 1000, 1, 96, 96
    data = np.random.rand(1000, 1, 96, 96)

    # loop over them
    data_equalized = np.zeros(data.shape)
    for i in range(data.shape[0]):
        image = data[i, 0, :, :]
        data_equalized[i, 0, :, :] = image_histogram_equalization(image)[0]

- Trilarion

---

对于那些想知道的人，图像的归一化直方图是指将图像的直方图除以图像中像素的总数，可以被视为每个灰度级的概率密度函数，这正是 density=True 所做的。因此，在这里 image_histogram 实际上就是归一化直方图。 - mh sattarian

这太棒了，感谢您提供可复制的代码。正是我所需要的。另外加上一，因为这可以适用于N维数组，而不仅仅是图像。 - Nick Crews

cdf = 255 * cdf / cdf[-1] 中的 255 应该改为 (number_bins - 1)，这样更加一致。 - Dominik Stańczak

@DominikStańczak0 谢谢您的观察。您是正确的，我已经编辑了代码。 - Trilarion

Answer 2

使用skimage模块提供的累积分布函数非常简单快捷。基本上，你需要数学证明它。

from skimage import exposure
import numpy as np
def histogram_equalize(img):
    img = rgb2gray(img)
    img_cdf, bin_centers = exposure.cumulative_distribution(img)
    return np.interp(img, bin_centers, img_cdf)

- Dodo

---

这对我很有用。它返回一个0-1浮点矩阵，您可以将其应用为掩码。 - Rutrus

Answer 3

截至今天，janeriksolem的网址已经失效。

然而，我找到了this gist，它链接到相同的页面，并声称可以在不计算直方图的情况下执行直方图均衡化。

代码如下：

img_eq = np.sort(img.ravel()).searchsorted(img)

- Diego Ferri

---

这种方法比直方图方法慢大约3倍。 - RuRo

Answer 4

这是一个针对单通道图像的快速替代实现。请参考skimage.exposure.histogram。使用timeit，Trilarion答案中的'image_histogram_equalization'的平均执行时间为0.3696秒，而这个函数的平均执行时间为0.0534秒。但是，这个实现也依赖于skimage。

import numpy as np
from skimage import exposure

def hist_eq(image):
    hist, bins = exposure.histogram(image, nbins=256, normalize=False)
    # append any remaining 0 values to the histogram
    hist = np.hstack((hist, np.zeros((255 - bins[-1])))) 
    cdf = 255*(hist/hist.sum()).cumsum()
    equalized = cdf[image].astype(np.uint8)

    return equalized

- Isaac Berrios