直方图均衡原理 & MATLAB实现

结合 Gonzalez 的《数字图像处理》第 3.3.1 节，在这里总结一下直方图均衡的原理、具体实现及代码。

直方图均衡

直方图均衡（Histogram Equalization）是一种利用灰度变换自动调节图像对比度的方法，通过灰度级的概率密度函数求出灰度变换函数，它是一种以累积分布函数变换法为基础的直方图修正法。

基本思想： 把原始图像的灰度分布直方图变换为均匀分布的形式（有展开直方图的趋势），扩大像素灰度值的动态范围，从而增强图像对比度。

用途： 在图像处理领域用于增强图像的对比度

缺点： 直方图均衡对处理的数据不加以选择，它可能会增加背景噪声的对比度并且降低有用信息的对比度。
另外，直方图处理不允许造成新的灰度级，所以在实际的直方图均衡应用中很少有完美平坦的直方图。
还有，我们的图像灰度级通常是离散的，所以不能保证离散直方图均衡还可以得到均匀的直方图。

1. 数学表示

首先考虑灰度值连续的情况：

• $r$ 表示待处理图像的灰度，取值区间假设为 $[0,L-1]$
• $s$ 表示处理后图像的灰度，取值区间仍然需要在 $[0,L-1]$
• $s=T(r)$ 为灰度映射函数，对于输入图像的每一个具有灰度值 $r$ 的像素，经过变换输出一个灰度值 $s$

定义变换函数 $T(r)$ 为：
$$s=T(r)=(L-1){\int }_0^r{p}_r(w)\mathrm{d}w$$

灰度值 $r$ 可以看作是区间 $[0,L-1]$ 内的随机变量，所以可以再记：

• $p_r(r)$ 是随机变量 $r$ 的概率密度函数
• $p_s(s)$ 是随机变量 $s$ 的概率密度函数

那么从 $p_r(r)$ 到 $p_s(s)$ 有一个变换：
$$p_s(s)=p_r(r)|\frac{dr}{ds}|$$
从上式可以看出，输出灰度 $s$ 的概率密度函数是由输入灰度 $r$ 的概率密度函数和所用的变换函数 $T(r)$ 共同决定的。

经过简单计算：
$$\frac{ds}{dr}=\frac{dT(r)}{dr}=(L-1)\frac{d}{dr}[{\int }_0^r{p}_r(w)\mathrm{d}w]=(L-1)p_r(r)$$
可得：
$$p_s(s)=p_r(r)|\frac{dr}{ds}|=p_r(r)|\frac{ds}{dr}{|}^{-1}=p_r(r)\frac{1}{(L-1)p_r(r)}=\frac{1}{L-1}$$
神奇的是，发现输出灰度值 $s$ 的概率密度竟然和输入灰度值 $r$ 没有关系，它是服从均匀分布的变量！

再看灰度值离散的情况：

一张灰度级范围为 $[0,L-1]$ 的数字图像中某一级灰度出现的次数是离散函数 $h(r_k)=n_k$，k = 0,1,…,L-1

• $r_k$ 是第 $k$ 级灰度值
• $n_k$ 是图像中灰度为 $r_k$ 的像素个数

归一化：$p(r_k)=n_k/MN$，此时 $p(r_k)$ 是灰度级 $r_k$ 在图像中出现的概率估计（直方图值），归一化的直方图所有分量之和等于 1

• $MN$ 表示图像的像素总数
• $M$ 是图像的行数
• $N$ 是图像的列数

在离散问题里，直方图其实就是概率 $p_r(r_k)$ 的图形。

1. PMF 概率质量函数（Probability Mass Function）：PMF 是离散随机变量在各特定点取值上的概率。

【概率质量函数与概率密度函数 (PDF) 的不同之处在于前者是对离散型随机变量定义的，而后者是对连续型随机变量】

一个离散的灰度图像 $X$，$n_k$ 表示灰度 $r_k=i$ 出现的次数，$MN$ 表示图像中所有像素数，$i$ 表示某一个灰度值，则图像在灰度 $r_k=i$ 上的概率质量函数为：
$$PM{F}_X(i)=\frac{n_k}{MN},0\le i<L$$
其中，$L$ 为图像的灰度数，通常为 256。

2. CDF 累积分布函数（Cumulative Distribution Function）：又叫分布函数，是对概率密度函数（概率质量函数）的积分（求和），能完整描述一个实随机变量 X 的概率分布。
$$CD{F}_X(x)=P(X\le x)$$

对一个离散的灰度图像 $X$，对应 $PM{F}_X$ 的累积分布函数为：
$$CD{F}_X(i)=\sum _{j=0}^{i}PM{F}_X(j)$$

将连续情况的下的 $T(r)$ 转换到离散情况下，得到变换函数 $T(r_k)$：
$$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,...,L-1$$
变换（映射）$T(r_k)$ 称为直方图均衡。

2. 处理步骤

1. 求出给定待处理图像的直方图 PMF

2. 由 PMF 求出 (L - 1) CDF

3. 将 (L - 1) CDF 映射为新的灰度值

MATLAB实现

MATLAB 实现直方图均衡化处理的工具是：J = histeq（I,n）

• I 为输入的原图像
• J 为直方图均衡化后的图像
• n 为均衡化后的灰度级数，默认值为 64

利用 MATLAB 工具包进行直方图处理：

% 利用工具 histeq() 实现直方图均衡化

close all;
clear all;
clc;

% 读取原灰度图
R = imread('lena_gray.jpg');
% 直方图均衡化
S = histeq(R);

figure(1)
subplot(121),imshow(uint8(R)),title('原灰度图');
subplot(122),imshow(uint8(S)),title('均衡化后');

figure(2)
subplot(121),imhist(R,64),title('原图像直方图');
subplot(122),imhist(S,64),title('均衡化后的直方图');

自己实现一下直方图均衡代码：

第一种：

% 读取原图
R = imread('lena_gray.jpg');
[row, col] = size(R);

% 显示原图和原图对应的直方图
subplot(2,2,1), imshow(R), title('原灰度图');
subplot(2,2,2), imhist(R), title('原灰度图的直方图');

% 计算原图像素各灰度值的PMF
PMF = zeros(1, 256);
for i = 1:row
    for j = 1:col
        PMF(R(i,j) + 1) = PMF(R(i,j) + 1) + 1; % R(i,j)为像素的灰度值
    end
end

% 计算CDF
CDF = zeros(1, 256);
CDF(1) = PMF(1);
for i = 2:256
    CDF(i) = CDF(i - 1) + PMF(i);
end

% 将CDF映射到新的灰度值
for i = 1:256
    Map(i) = round((CDF(i) - 1) * 255 / (row * col));
end
for i = 1:row
    for j = 1:col
        R(i,j) = Map(R(i,j) + 1);
    end
end

subplot(2,2,3), imshow(R), title('直方图均衡');
subplot(2,2,4), imhist(R), title('均衡后的直方图');

第二种：（感觉更好理解一些）

% 读取原图
R = imread('lena_gray.jpg');
[row, col] = size(R);

% 显示原图和原图对应的直方图
subplot(2,2,1), imshow(R), title('原灰度图');
subplot(2,2,2), imhist(R), title('原灰度图的直方图');

% 计算PMF，即统计各灰度值的像素数量
PMF = zeros(1, 256); %注意MATLAB的数组索引从1开始
for i = 1:row
    for j = 1:col
        PMF(R(i,j) + 1) = PMF(R(i,j) + 1) + 1; % R(i,j)为像素的灰度值
    end
end
PMF = PMF / (row * col);

% 计算CDF
CDF = zeros(1,256);
CDF(1) = PMF(1);
for i = 2:256
    CDF(i) = CDF(i - 1) + PMF(i);
end

% 计算均衡后的像素值
Sk = zeros(1,256);
for i = 1:256
    Sk(i) = CDF(i) * 255;
end

% 映射到新的像素值
Sk = round(Sk);
for i = 1:row
    for j = 1:col
        R(i,j) = Sk(R(i,j) + 1);
    end
end

% 绘制直方图均衡后的图像
subplot(2,2,3), imshow(R), title('直方图均衡');
subplot(2,2,4), imhist(R), title('均衡后的直方图');