直方图均衡 & MATLAB实现
结合 Gonzalez 的《数字图像处理》第 3.3.1 节，在这里总结一下直方图均衡的原理、具体实现及代码。
直方图均衡
直方图均衡（Histogram Equalization）是一种利用灰度变换自动调节图像对比度的方法，通过灰度级的概率密度函数求出灰度变换函数，它是一种以累积分布函数变换法为基础的直方图修正法。
基本思想： 把原始图像的灰度分布直方图变换为均匀分布的形式（有展开直方图的趋势），扩大像素灰度值的动态范围，从而增强图像对比度。
用途： 在图像处理领域用于增强图像的对比度
缺点： 直方图均衡对处理的数据不加以选择，它可能会增加背景噪声的对比度并且降低有用信息的对比度。 还有，我们的图像灰度级通常是离散的，所以不能保证离散直方图均衡还可以得到均匀的直方图。
1. 数学表示
首先考虑灰度值连续的情况：
r
r
r 表示待处理图像的灰度，取值区间假设为
[
0
,
L
−
1
]
[0,L-1]
[0,L−1]</li><li>
s
s
s 表示处理后图像的灰度，取值区间仍然需要在
[
0
,
L
−
1
]
[0,L-1]
[0,L−1]</li><li>
s
=
T
(
r
)
s = T(r)
s=T(r) 为灰度映射函数，对于输入图像的每一个具有灰度值
r
r
r 的像素，经过变换输出一个灰度值
s
s
s</li>
定义变换函数
T
(
r
)
T(r)
T(r)</mark> 为：<br>
s
=
T
(
r
)
=
(
L
−
1
)
∫
0
r
p
r
(
w
)
d
w
s = T(r) = (L-1)\int^r_0p_r(w){\rm d}w
s=T(r)=(L−1)∫0rpr(w)dw
灰度值
r
r
r 可以看作是区间
[
0
,
L
−
1
]
[0,L-1]
[0,L−1] 内的随机变量，所以可以再记：
p
r
(
r
)
p_r(r)
pr(r) 是随机变量
r
r
r 的概率密度函数</li><li>
p
s
(
s
)
p_s(s)
ps(s) 是随机变量
s
s
s 的概率密度函数</li>
那么从
p
r
(
r
)
p_r(r)
pr(r) 到
p
s
(
s
)
p_s(s)
ps(s) 有一个变换：<br>
p
s
(
s
)
=
p
r
(
r
)
∣
d
r
d
s
∣
p_s(s) = p_r(r)|\frac{dr}{ds}|
ps(s)=pr(r)∣dsdr∣<br> 从上式可以看出，输出灰度
s
s
s 的概率密度函数是由输入灰度
r
r
r 的概率密度函数和所用的变换函数
T
(
r
)
T(r)
T(r) 共同决定的。
经过简单计算：
d
s
d
r
=
d
T
(
r
)
d
r
=
(
L
−
1
)
d
d
r
[
∫
0
r
p
r
(
w
)
d
w
]
=
(
L
−
1
)
p
r
(
r
)
\frac{ds}{dr} = \frac{dT(r)}{dr} = (L-1)\frac{d}{dr}[\int^r_0p_r(w){\rm d}w] = (L-1)p_r(r)
drds=drdT(r)=(L−1)drd[∫0rpr(w)dw]=(L−1)pr(r)
p
s
(
s
)
=
p
r
(
r
)
∣
d
r
d
s
∣
=
p
r
(
r
)
∣
d
s
d
r
∣
−
1
=
p
r
(
r
)
1
(
L
−
1
)
p
r
(
r
)
=
1
L
−
1
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}
ps(s)=pr(r)∣dsdr∣=pr(r)∣drds∣−1=pr(r)(L−1)pr(r)1=L−11<br> 神奇的是，发现输出灰度值
s
s
s 的概率密度竟然和输入灰度值
r
r
r 没有关系，它是服从**均匀分布**的变量！
再看灰度值离散的情况：
一张灰度级范围为
[
0
,
L
−
1
]
[0, L-1]
[0,L−1] 的数字图像中某一级灰度出现的次数是离散函数
h
(
r
k
)
=
n
k
h(r_k) = n_k
h(rk)=nk，k = 0,1,…,L-1
r
k
r_k
rk 是第
k
k
k 级灰度值</li><li>
n
k
n_k
nk 是图像中灰度为
r
k
r_k
rk 的像素个数</li>
归一化：
p
(
r
k
)
=
n
k
/
M
N
p(r_k) = n_k / MN
p(rk)=nk/MN，此时
p
(
r
k
)
p(r_k)
p(rk) 是灰度级
r
k
r_k
rk 在图像中出现的**概率估计（直方图值）**，归一化的直方图所有分量之和等于 1
M
N
MN
MN 表示图像的像素总数</li><li>
M
M
M 是图像的行数</li><li>
N
N
N 是图像的列数</li>
在离散问题里，直方图其实就是概率
p
r
(
r
k
)
p_r(r_k)
pr(rk) 的图形。
1. PMF 概率质量函数（Probability Mass Function）：PMF 是离散随机变量在各特定点取值上的概率。
【概率质量函数与概率密度函数 (PDF) 的不同之处在于前者是对离散型随机变量定义的，而后者是对连续型随机变量】
一个离散的灰度图像
X
X
X，
n
k
n_k
nk 表示灰度
r
k
=
i
r_k = i
rk=i 出现的次数，
M
N
MN
MN 表示图像中所有像素数，
i
i
i 表示某一个灰度值，则图像在灰度
r
k
=
i
r_k = i
rk=i 上的概率质量函数为：<br>
P
M
F
X
(
i
)
=
n
k
M
N
,
0
≤
i
<
L
PMF_X(i)=\frac{n_k}{MN}, 0 \leq i < L
PMFX(i)=MNnk,0≤i<L<br> 其中，
L
L
L 为图像的灰度数，通常为 256。
2. CDF 累积分布函数（Cumulative Distribution Function）：又叫分布函数，是对概率密度函数（概率质量函数）的积分（求和），能完整描述一个实随机变量 X 的概率分布。
C
D
F
X
(
x
)
=
P
(
X
≤
x
)
CDF_X(x)=P(X \leq x)
CDFX(x)=P(X≤x)
对一个离散的灰度图像
X
X
X，对应
P
M
F
X
PMF_X
PMFX 的累积分布函数为：<br>
C
D
F
X
(
i
)
=
∑
j
=
0
i
P
M
F
X
(
j
)
CDF_X(i) = \sum_{j=0}^{i}PMF_X(j)
CDFX(i)=j=0∑iPMFX(j)
将连续情况的下的
T
(
r
)
T(r)
T(r) 转换到离散情况下，得到<mark>变换函数
T
(
r
k
)
T(r_k)
T(rk)</mark>：<br>
s
k
=
T
(
r
k
)
=
(
L
−
1
)
∑
j
=
0
k
p
r
(
r
j
)
=
(
L
−
1
)
M
N
∑
j
=
0
k
n
j
，
k
=
0
,
1
,
2
,
.
.
.
,
L
−
1
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
sk=T(rk)=(L−1)j=0∑kpr(rj)=MN(L−1)j=0∑knj，k=0,1,2,...,L−1<br> 变换（映射）
T
(
r
k
)
T(r_k)
T(rk) 称为<mark>直方图均衡</mark>。
2. 处理步骤
- 求出给定待处理图像的直方图 PMF 1. 由 PMF 求出 (L - 1) CDF 1. 将 (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(\'均衡后的直方图\');
请发表评论