直方图均衡化

引言

图像的空域处理是一种重要的图像处理技术，这类方法直接以图像的像素操作为基础，主要分为灰度变换和空域滤波两大类，直方图均衡化（Histogram equalization）就是一种常用的灰度变换方法。

直方图

对于灰度级（intensity levels）范围为$[0,L-1]$的数字图像，其直方图可以表示为一个离散函数$h(r_k)=n_k$， 其中$r_k$是第k级灰度值（intensity value），$n_k$是图像中灰度值为$r_k$的像素个数，也就是说，图像的灰度直方图表征的是该图像的灰度分布。在实际应用中，通常对直方图进行归一化再进行后续处理，假设灰度图像的维数是M×N，MN表示图像的像素总数，则归一化直方图可以表示为
$$p(r_k)=n_k/MN，k=0,1,\dots ,L-1$$也就是说，$p(r_k)$表示灰度级$r_k$在图像中出现概率的估计，归一化直方图所有分量之和等于1。

直方图均衡化

1. 问题描述

通常，暗图像直方图的分量集中在灰度较低的一端，而亮图像直方图分量偏向于灰度较高的一端，如下图

从图中可以得到这样的结论：如果一幅图像的灰度直方图几乎覆盖了整个灰度的取值范围，并且除了个别灰度值的个数较为突出，整个灰度值分布近似于均匀分布，那么这幅图像就具有较大的灰度动态范围和较高的对比度，同时图像的细节更为丰富。已经证明，仅仅依靠输入图像的直方图信息，就可以得到一个变换函数，利用该变换函数可以将输入图像达到上述效果，该过程就是直方图均衡化。

下面从数学的角度理解直方图均衡化。假设待处理图像为灰度图像，$r$表示待处理图像的灰度，取值范围为$[0,L-1]$，则$r=0$表示黑色，$r=L-1$表示白色，直方图均衡化的过程对应于一个变换$T$：$$s=T(r),0\le r\le L-1i\text{\hspace{1em}(1)}$$也就是说，对于输入图像的某个灰度值$r$，可以通过变换$T$得到均衡化后的图像对应位置的灰度值$s$。其中变换$T$满足以下条件：
(a) $T(r)$在$[0,L-1]$上严格单调递增；
(b) 当$0\le r\le L-1$时，$0\le T(r)\le L-1$。

条件(a)中要求$T(r)$严格单调递增是为了保证输出灰度值与输入灰度值一一对应，同时像素灰度值之间的相对大小关系不变，这样可以避免反变换时出现问题；条件(b)保证了输出图像的灰度范围与输入图像相同。实际中处理的图像通常是整数灰度值，必须把所有结果四舍五入为最接近的整数值。因此，当严格单调条件不满足时，使用寻找最接近整数匹配的方法解决反变换不唯一的问题。

2. 均衡化推导

首先考虑连续形式，一幅灰度图像的灰度级可以看作区间 $[0,L-1]$ 内的随机变量，因此可用其概率密度函数(PDF)描述。假设 $p_r(r)$ 和 $p_s(s)$ 分别表示随机变量 $r$ 和 $s$ 的PDF，$p_r(r)$ 和变换 $T$ 已知，且 $T(r)$ 在定义域内连续可微，则变换后 $s$ 的PDF可由下式得到：$$p_s(s)=p_r(r)|\frac{dr}{ds}|\text{\hspace{1em}(2)}$$

此处用到了概率论中有关随机变量函数分布的相关定理：设 $\xi$ 是连续型随机变量，其概率密度函数为 $p(x)$，又函数 $y=f(x)$
严格单调，其反函数 $h(y)$ 有连续导数，则 $\eta =f(\xi )$ 也是一个连续型随机变量，其概率密度函数为
$$\psi (y)=\{\begin{array}{ll}p[h(y)]\cdot |h^{\prime }(y)|& ,\alpha <y<\beta \\ 0& ,其他\end{array}$$ 其中 α = m i n { f ( − ∞ ) , f ( + ∞ ) } β = m a x { f ( − ∞ ) , f ( + ∞ ) } \begin{aligned} \alpha=min\{f(-\infty),f(+\infty)\}\\ \beta =max\{f(-\infty),f(+\infty)\} \end{aligned} α=min{f(−∞),f(+∞)}β=max{f(−∞),f(+∞)}​

由此看到，输出图像灰度 $s$ 的PDF就由输入图像灰度 $r$ 的PDF和变换 $T$ 得到。

图像处理中一个重要的变换函数构造如下 $$s=T(r)=(L-1){\int }_0^r{p}_r(w)\mathrm{d}x\text{\hspace{1em}(3)}$$
其中，$w$是形式积分变量，公式右边是随机变量$r$的累积分布函数(CDF)，可以看到，(3)式完全满足前述的条件(a)和(b)。

为寻找变换后随机变量 $s$ 的概率密度函数 $p_s(s)$，由(2)式得，
$$\begin{array}{rl}\frac{\mathrm{d}s}{\mathrm{d}r}& =\frac{\mathrm{d}T(r)}{\mathrm{d}r}\\ & =(L-1)\frac{\mathrm{d}}{\mathrm{d}r}{\int }_0^r{p}_r(w)\mathrm{d}w=(L-1)p_r(r)\end{array}\text{\hspace{1em}(4)}$$将(4)式代入(2)式得
$$\begin{array}{rl}{p}_s(s)& =p_r(r)|\frac{\mathrm{d}r}{\mathrm{d}s}|\\ & =p_r(r)\frac{1}{(L-1)p_r(r)}=\frac{1}{L-1},0\le s\le L-1\end{array}\text{\hspace{1em}(5)}$$由(5)式可知，$p_s(s)$ 为均匀分布，也就是说，输入图像的PDF经过(3)式中的变换 $T$ 后得到的随机变量 $s$ 服从均匀分布。
结论：图像均衡化变换 $T(r)$ 取决于 $p_r(r)$， 但得到的 $p_s(s)$ 始终是均匀的，与 $p_r(r)$ 的形式无关。
对于离散形式，其推导过程与连续形式相似，用概率直方图和求和运算分别代替概率密度函数和积分运算，可得(3)式的离散形式
$$\begin{array}{rl}{s}_k& =T(r_k)\\ & =(L-1)\sum _{j=0}^k{p}_r(r_j)=\frac{L-1}{MN}\sum _{j=0}^k{n}_j,k=0,1,2,\cdots ,L-1\end{array}\text{\hspace{1em}(6)}$$其中，$MN$为图像像素总数，$n_k$ 是灰度为 $r_k$ 的像素个数，$L$ 是图像可能的灰度级数量(例如对于8比特图像 $L=256$)。通过(6)式，输出图像中像素的灰度值可由输入图像中像素灰度 $r_k$ 映射为 $s_k$ 后得到。(6) 式中的变换（映射）$T(r_k)$ 称为直方图均衡或直方图线性变换。

3. 例子

假设一幅大小为 $64\times 64$ 像素的3比特图像 $(L=8)$ 的灰度分布如表1所示，其中灰度级是范围 $[0,L-1]=[0,7]$ 中的整数。

$r_k$	$n_k$	$p_k(r_k)=n_k/MN$
$r_0=0$	790	0.19
$r_1=1$	1023	0.25
$r_2=2$	850	0.21
$r_3=3$	656	0.16
$r_4=4$	329	0.08
$r_5=5$	245	0.06
$r_6=6$	122	0.03
$r_7=7$	81	0.02

则直方图均衡化的值可由 (6) 式计算得到，例如
$$\begin{array}{rl}{s}_0& =T(r_0)=7\sum _{j=0}^0{p}_r(r_j)\\ & =7p_r(r_0)=7\times 0.19=1.33\end{array}$$同理，$$\begin{array}{rl}{s}_1& =T(r_1)=7\sum _{j=0}^1{p}_r(r_j)\\ & =7p_r(r_0)+7p_r(r_1)=7\times 0.19+7\times 0.25=3.08\end{array}$$以此类推计算求得所有 $s_k$ 值并将其近似为最接近的整数，结果如下：
$$\begin{array}{rlrl}& s_0=1.33\to 1,& & s_4=6.23\to 6\\ & s_1=3.08\to 3,& & s_5=6.65\to 7\\ & s_2=4.55\to 5,& & s_6=6.86\to 7\\ & s_3=5.67\to 6,& & s_7=7.00\to 7\end{array}$$可以看到，均衡化后的图像只有5个灰度级，$r_0=0$ 被映射为 $s_0=1$，在均衡化后的图像中有790个像素取该值（见表1），有1023个像素取 $s_1=3$，有850个像素取 $s_2=5$，有656+329=985个像素取6，因为 $r_3$ 和 $r_4$ 都被映射为6这个值，同理，有245+122+81=448个像素取7。
因为直方图是PDF的近似，并且处理后不产生新的灰度级，所以在实际的直方图均衡化很少能够得到完全平坦的直方图。在离散情况下，通常不能证明离散的直方图均衡化能得到均匀的直方图。尽管如此，均衡化有助于图像直方图的延展，均衡化后图像的灰度级范围更宽，有效地增强了图像的对比度。同时，上述方法完全是“自动”的，仅利用输入图像的直方图信息，无需更多的参数。

Reference

Gonzalez, Rafael C., and Richard E. Woods. “Digital image processing.”