图像背后的数学--直方图

从最开始接触图像处理，到现在，也有很多年了，现在回过来看以前学的很多东西，会开始慢慢尝试理解现象背后的本质，以前只是单纯地从技术的角度去学习图像处理的很多算法，随着知识的积累，会感到很多算法的背后其实都有着看似简单却又深刻的数学与物理原理。今天简单介绍图像中常用的一个统计–直方图。

直方图统计是图像处理中非常基本的一种统计，简单来说，就是给定一张图像 I∈[0,255]，我们要统计这张图像里，像素值为 0 的像素有多少个，像素值为 1 的有多少个，以此类推，直到像素值为 255，我们可以假设图像总的像素个数为 N，像素值为 i 的像素个数为 ni， 那么我们有：

N=∑i=0255ni

我们以前学概率的时候，都知道一个经典的抛硬币的思想实验，抛硬币是二分类 {0, 1} 的情形，还有一个经典的抛骰子的实验，这是多分类的情形 {1, …, 6}，我们可以把这个思想实验再进一步的拓展，考虑更多分类的情形，比如 {0, 1, 2, …, 255}，这样的话，就是 256 种情况。每一个像素值，都相当于一个类别，一张图像总的像素数就是一个样本总体，我们想知道每一个像素值所占的比重，或者说概率，就用这个像素值的像素个数除以总的像素数，即:

pi=niN

所以说，直方图，反映的其实是像素值的一种概率密度分布。

我们以前学过，给定概率密度函数 p(x)，想求某一个区间的概率，可以用如下的积分表达式：

p(x∈(a,b))=∫bap(x)dx

对于离散概率分布，我们可以直接用求和的表达式：

p(x∈(a,b))=∑abp(x)

所以说，要求直方图某一段像素值区间的累积概率，可以表示为：

p(i∈(a,b))=∑i=abpi

而这个就是我们常说的累积直方图，一般我们都是从 0 开始， 比如说，我们要求像素值 0 到 k 之间的累积直方图，可以表示为：

ck=∑i=0kpi

现在我们有了直方图，有了累积直方图，接下来，我们来看看直方图均衡，直方图均衡是图像动态范围拉伸的一种方法，我们知道，要将一个区间 [a, b] 映射到另外一个区间 [c, d]， 最简单的方法是线性映射:

y=x−ab−a×(d−c)+c,x∈[a,b],y∈[c,d]

目前最常用的图像编码都是 8 bits，所以图像的动态范围都是 [0, 255]，一个图像的对比度要好，动态范围要尽可能地大，同时直方图的分布不能太窄，如果都挤在一起，图像的对比度肯定不好。直方图均衡是图像处理中，非常经典的一个拉伸动态范围调整图像对比度的算法。就是利用累积直方图的性质，我们知道累积直方图的区间为 [0, 1]，并且是单调递增的，我们假设图像原来的灰度值为 i，其对应的累积直方图为 ci， 0≤ci≤1, 那么我们可以求出其映射到新的灰度值为:

y=ci−01−0×(255−0)+0=ci∗255

其实，就是将累积直方图的区间 [0, 1] 映射到新的区间 [0, 255]，利用累积直方图单调递增的特性，保证了映射前后灰度值的大小关系没有发生改变，就是说映射前 i<j，那么映射后的灰度值同样满足 yi<yj, 我们来看一小段程序：

import numpy as np
import matplotlib.pyplot as plt
from skimage import data

img = data.moon()

rows, cols = img.shape
N = rows * cols - 1

Hist_ori = np.zeros((256, 1))
Hist_cmu = np.zeros((256, 1))

for ii in range(256):
    mask_1 = (img == ii)
    Hist_ori[ii] = np.sum(mask_1) * 1.0 / N
    if ii == 0:
        Hist_cmu[ii] = Hist_ori[ii]
    else:
        Hist_cmu[ii] = Hist_cmu[ii-1] + Hist_ori[ii]

Hist_eq = np.zeros((256, 1))
img_out = img.copy()

for ii in range(rows):
    for jj in range(cols):
        img_out[ii, jj] = math.floor(Hist_cmu[img[ii, jj]] * 255.0)
        Hist_eq[img_out[ii, jj]] = Hist_eq[img_out[ii, jj]] + 1

plt.figure()
plt.axis('off')
plt.subplot(1,2,1)
plt.imshow(img, plt.cm.gray)
plt.subplot(1,2,2)
plt.imshow(img_out, plt.cm.gray)
plt.subplot(2,2,3)
plt.plot(Hist_ori)
plt.subplot(2,2,4)
plt.plot(Hist_eq*1.0/N)

plt.show()

看看效果图：

可以看到直方图被明显的拉伸开来了。