影象直方圖均衡化詳解
直方圖均衡化是一種經典的影象處理演算法，用以改善影象的亮度和對比度。如下圖，亮度整體偏暗，灰度值大致分佈在較低的範圍內（雖然是作者刻意畫成這樣的（￣▽￣）“）：
對影象進行直方圖均衡化的目的是，使其原本分佈集中的畫素值，均衡的分佈到所有可取值的範圍，這樣，影象就既有明亮也有灰暗，對比度和亮度就得到了改善。
直方圖
影象的直方圖就是畫素值出現的總的頻數，其公式如下：
一般情況下，通常將其除以總的畫素數，透過機率的形式進行表述：
全域性直方圖均衡化
1.原理
直方圖均衡化可以看作是一種對映，對原圖中的畫素點施加對映使其轉換為圖。在直方圖均衡化中，我們保持畫素值的單調性，即：
其中和為中的畫素值，和為中的畫素值。由和之間的畫素值是一一對應且單調的，可以得出：
即區間中的畫素數和區間中的畫素數是相等的。因此，每個對應的畫素值的累積機率函式是相等的：
我們希望對映後的圖直方圖分佈是均衡的，即是理想情況下，圖中畫素值的機率分佈是服從均勻分佈的，即，則：
這樣我們就得到了實現均衡化的對映方式。
2. C++實現
實現均衡化的程式碼如下：
// Statistical histogram
void
statHist
（
const
cv
：：
Mat
&
image
，
int
*
hist
）
{
int
idx
=
0
；
for
（
int
i
=
0
；
i
<
image
。
rows
；
i
++
）
{
for
（
int
j
=
0
；
j
<
image
。
cols
；
j
++
）
{
idx
=
int
（
image
。
at
<
uchar
>
（
i
，
j
））；
hist
［
idx
］
++
；
}
}
}
// Get map table
void
mapTable
（
const
int
*
hist
，
const
int
size
，
uchar
*
table
）
{
float
factor
=
256。f
/
size
；
int
map_value
，
cumu_num
=
0
；
for
（
int
i
=
0
；
i
<
256
；
i
++
）
{
cumu_num
+=
hist
［
i
］；
map_value
=
int
（
cumu_num
*
factor
）；
if
（
map_value
>
255
）
map_value
=
255
；
table
［
i
］
=
uchar
（
map_value
）；
}
}
// Mapping
void
map
（
const
cv
：：
Mat
&
src
，
const
uchar
*
table
，
cv
：：
Mat
&
dst
）
{
src
。
copyTo
（
dst
）；
for
（
int
i
=
0
；
i
<
src
。
rows
；
i
++
）
{
for
（
int
j
=
0
；
j
<
src
。
cols
；
j
++
）
{
dst
。
at
<
uchar
>
（
i
，
j
）
=
table
［（
int
）
src
。
at
<
uchar
>
（
i
，
j
）］；
}
}
}
// Global histogram equlization
void
histEqualization
（
const
cv
：：
Mat
&
src
，
cv
：：
Mat
&
dst
）
{
int
hist
［
256
］
=
{
0
}；
uchar
table
［
256
］
=
{
0
}；
statHist
（
src
，
hist
）；
mapTable
（
hist
，
src
。
rows
*
src
。
cols
，
table
）；
map
（
src
，
table
，
dst
）；
}
3. 處理效果
從下圖中可以看出，全域性直方圖均衡化確實增強了對比度和亮度，但是得到的效果有些過分亮了，這是因為原圖中大多畫素值集中於較小的區間內，稍微大一點兒的畫素值對應的累積機率都接近1了，對映過後的畫素值自然會很亮。除此之外，原本柔和的較暗區域內出現了明顯的邊緣現象，這是因為原本相鄰的相差較小畫素，對映後的差值被拉得過大，自然出現邊緣效應。為了避免較亮的畫素點被過分增強，可以考慮對該畫素點周圍區域性區域進行直方圖均衡化，從而不被其他區域較暗的畫素點影響。
區域性直方圖均衡化
1.簡單分塊直方圖均衡化
一種簡單的區域性直方圖均衡化的方法是，直接將影象分為若干的小塊，每個小塊分別進行直方圖均衡化。程式碼如下：
void
localHistEqualization
（
const
cv
：：
Mat
&
src
，
cv
：：
Mat
&
dst
，
const
int
m
）
{
// Padding
cv
：：
Mat
im_pad
；
const
int
h
=
src
。
rows
，
w
=
src
。
cols
；
const
int
rdh
=
h
%
m
，
rdw
=
w
%
m
；
const
int
pdh
=
rdh
==
0
？
h
：
h
+
m
-
rdh
；
const
int
pdw
=
rdw
==
0
？
w
：
w
+
m
-
rdw
；
const
int
res_h
=
pdh
-
h
，
res_w
=
pdw
-
w
；
const
int
top
=
res_h
/
2
，
left
=
res_w
/
2
；
const
int
bottom
=
res_h
-
top
，
right
=
res_w
-
left
；
cv
：：
copyMakeBorder
（
src
，
im_pad
，
top
，
bottom
，
left
，
right
，
cv
：：
BORDER_DEFAULT
）；
// Equalization
int
i
=
0
，
j
=
0
，
ni
=
m
，
nj
=
m
；
while
（
i
<
pdh
）
{
ni
=
i
+
m
；
while
（
j
<
pdw
）
{
nj
=
j
+
m
；
equalizeHist
（
im_pad
。
rowRange
（
i
，
ni
）。
colRange
（
j
，
nj
），
im_pad
。
rowRange
（
i
，
ni
）。
colRange
（
j
，
nj
））；
j
=
nj
；
}
i
=
ni
；
j
=
0
；
}
im_pad
。
rowRange
（
top
，
h
+
top
）。
colRange
（
left
，
w
+
left
）。
copyTo
（
dst
）；
}
結果那是稍微有點兒抽象：
2.區域性自適應直方圖均衡化
雖然我還挺喜歡上面這種風格的，但是很明顯，直接分為若干小塊分別進行直方圖均衡化並不是一個好辦法。不僅存在一些過度增強，塊與塊之間還存在非常顯著的差異。那如果我們不區分塊與塊，而是以每個畫素為中心，根據其周圍區域性區域對該畫素進行直方圖均衡化如何？這其實就是區域性自適應直方圖均衡化。使用此方法，每一個畫素都需要計算次運算，但是還好每次移動塊，只有塊的一邊出去，一邊進來，所以此方法可以加速。話不多說，直接上程式碼：
void
localAdaptiveHistEqualization
（
const
cv
：：
Mat
&
src
，
cv
：：
Mat
&
dst
，
const
int
m
）
{
// Padding
cv
：：
Mat
src_pad
；
const
int
psz
=
m
/
2
；
const
int
hend
=
src
。
rows
+
psz
；
const
int
wend
=
src
。
cols
+
psz
；
cv
：：
copyMakeBorder
（
src
，
src_pad
，
psz
，
psz
，
psz
，
psz
，
cv
：：
BORDER_DEFAULT
）；
src
。
copyTo
（
dst
）；
// Equalization
int
hist
［
256
］
=
{
0
}；
float
factor
=
256。f
/
（
m
*
m
）；
int
pi
=
0
，
pj
=
0
，
val
=
0
，
cumv
=
0
；
for
（
int
i
=
psz
；
i
<
hend
；
i
++
）
{
for
（
int
j
=
psz
；
j
<
wend
；
j
++
）
{
if
（
j
==
psz
）
{
// Reset
for
（
int
pi
=
0
；
pi
<
256
；
pi
++
）
{
hist
［
pi
］
=
0
；
}
// Statistical histogram
for
（
pi
=
i
-
psz
；
pi
<=
i
+
psz
；
pi
++
）
{
for
（
pj
=
j
-
psz
；
pj
<=
j
+
psz
；
pj
++
）
{
val
=
（
int
）
src_pad
。
at
<
uchar
>
（
pi
，
pj
）；
hist
［
val
］
++
；
}
}
}
else
{
// Pop left
pj
=
j
-
psz
-
1
；
for
（
pi
=
i
-
psz
；
pi
<=
i
+
psz
；
pi
++
）
{
val
=
（
int
）
src_pad
。
at
<
uchar
>
（
pi
，
pj
）；
hist
［
val
］
——
；
}
// Push right
pj
=
j
+
psz
；
for
（
pi
=
i
-
psz
；
pi
<=
i
+
psz
；
pi
++
）
{
val
=
（
int
）
src_pad
。
at
<
uchar
>
（
pi
，
pj
）；
hist
［
val
］
++
；
}
}
// Equlization
cumv
=
0
；
val
=
（
int
）
src_pad
。
at
<
uchar
>
（
i
，
j
）；
for
（
pi
=
0
；
pi
<=
val
；
pi
++
）
{
cumv
+=
hist
［
pi
］；
}
dst
。
at
<
uchar
>
（
i
-
psz
，
j
-
psz
）
=
uchar
（
cumv
*
factor
）；
}
}
}
結果如下：
啊這，，，塊效應是減輕了不少，可是過度增強的情況還是存在著。此外，即使減少了一些統計直方圖時的操作，對影象的每個元素進行操作的計算代價仍然過高。針對這些情況，研究人員提出限制對比度的直方圖均衡化方法。
限制對比度的直方圖均衡化
1. 限制對比度的全域性直方圖均衡化
影象中的畫素分佈並不均勻，可能存在個別畫素值的畫素數目特別多，進而使直方圖均衡化時比該值更大的畫素值被過度增強。因此，可以從這類畫素中“挪”出一些平均分配到其他畫素中。下面對全域性直方圖均衡化方法限制對比度：
void
che
（
const
cv
：：
Mat
&
src
，
cv
：：
Mat
&
dst
，
const
float
clip_limit
）
{
// Statistical histogram
const
int
h
=
src
。
rows
；
const
int
w
=
src
。
cols
；
int
hist
［
256
］
=
{
0
}；
uchar
table
［
256
］
=
{
0
}；
statHist
（
src
，
hist
）；
// Limitation
int
steal
=
0
；
const
int
limit
=
h
*
w
*
clip_limit
/
256
；
for
（
int
k
=
0
；
k
<
256
；
k
++
）
{
if
（
hist
［
k
］
>
limit
）
{
steal
+=
hist
［
k
］
-
limit
；
hist
［
k
］
=
limit
；
}
}
// Hand out the steals averagely
const
int
bonus
=
steal
/
256
；
for
（
int
k
=
0
；
k
<
256
；
k
++
）
{
hist
［
k
］
+=
bonus
；
}
// Get mapping table
mapTable
（
hist
，
h
*
w
，
table
）；
// Map
map
（
src
，
table
，
dst
）；
}
結果如下：
相比原圖對比度被增強了不少，相比全域性直方圖均衡化對比度被限制了不少，總的來說效果還行。
2. 限制對比度的自適應直方圖均衡化
同樣地，我們也想知道，如果在區域性自適應直方圖均衡化的基礎之上限制對比度會更好嗎？動手試試看！前面已經嘗試過，使用在區域性使用直方圖會產生塊效應，為了減輕塊效應，我們對相鄰塊之間進行了插值運算。
void
clahe
（
const
cv
：：
Mat
&
src
，
cv
：：
Mat
&
dst
，
const
float
clip_limit
，
const
int
m
）
{
// Padding
cv
：：
Mat
im_pad
；
const
int
h
=
src
。
rows
，
w
=
src
。
cols
；
const
int
pdh
=
minInteger
（
h
，
m
），
pdw
=
minInteger
（
w
，
m
）；
const
int
res_h
=
pdh
-
h
，
res_w
=
pdw
-
w
；
const
int
left
=
res_w
/
2
，
top
=
res_h
/
2
；
const
int
right
=
res_w
-
left
，
bottom
=
res_h
-
top
；
cv
：：
copyMakeBorder
（
src
，
im_pad
，
top
，
bottom
，
left
，
right
，
cv
：：
BORDER_DEFAULT
）；
// Initialize histogram
const
int
hnum
=
pdh
/
m
；
const
int
wnum
=
pdw
/
m
；
const
int
bnum
=
hnum
*
wnum
；
int
**
hists
=
new
int
*
［
bnum
］；
float
**
table
=
new
float
*
［
bnum
］；
for
（
int
i
=
0
；
i
<
bnum
；
i
++
）
{
hists
［
i
］
=
new
int
［
256
］；
table
［
i
］
=
new
float
［
256
］；
for
（
int
j
=
0
；
j
<
256
；
j
++
）
{
hists
［
i
］［
j
］
=
0
；
table
［
i
］［
j
］
=
0。f
；
}
}
// Statistical histogram， Get mapping table
const
int
benum
=
m
*
m
；
const
int
limit
=
benum
*
clip_limit
/
256
；
int
idx
=
0
，
steal
=
0
，
bonus
=
0
；
for
（
int
i
=
0
；
i
<
hnum
；
i
++
）
{
for
（
int
j
=
0
；
j
<
wnum
；
j
++
）
{
// Statistical
idx
=
i
*
wnum
+
j
；
statHist
（
im_pad
。
rowRange
（
i
*
m
，
（
i
+
1
）
*
m
）。
colRange
（
j
*
m
，
（
j
+
1
）
*
m
），
hists
［
idx
］）；
// Limitation
steal
=
0
；
for
（
int
k
=
0
；
k
<
256
；
k
++
）
{
if
（
hists
［
idx
］［
k
］
>
limit
）
{
steal
+=
hists
［
idx
］［
k
］
-
limit
；
hists
［
idx
］［
k
］
=
limit
；
}
}
// Hand out the steals averagely
bonus
=
steal
/
256
；
for
（
int
k
=
0
；
k
<
256
；
k
++
）
{
hists
［
idx
］［
k
］
+=
bonus
；
}
// Cumulative ——> table
table
［
idx
］［
0
］
=
hists
［
idx
］［
0
］
*
255。f
/
benum
；
for
（
int
k
=
1
；
k
<
256
；
k
++
）
{
table
［
idx
］［
k
］
=
table
［
idx
］［
k
-
1
］
+
hists
［
idx
］［
k
］
*
255。f
/
benum
；
}
}
}
delete
［］
hists
；
// Equalization and Interpolation
const
int
hm
=
m
/
2
；
//im_pad。copyTo（im_pad）；
int
hbi
=
0
，
wbi
=
0
；
int
bidx
［
4
］
=
{
0
}；
float
p
=
0。f
，
q
=
0。f
；
for
（
int
i
=
0
；
i
<
pdh
；
i
++
）
{
for
（
int
j
=
0
；
j
<
pdw
；
j
++
）
{
//four coners
if
（
i
<=
hm
&&
j
<=
hm
）
{
idx
=
0
；
im_pad
。
at
<
uchar
>
（
i
，
j
）
=
（
uchar
）（
table
［
idx
］［
im_pad
。
at
<
uchar
>
（
i
，
j
）］）；
}
else
if
（
i
<=
hm
&&
j
>=
pdw
-
hm
）
{
idx
=
wnum
-
1
；
im_pad
。
at
<
uchar
>
（
i
，
j
）
=
（
uchar
）（
table
［
idx
］［
im_pad
。
at
<
uchar
>
（
i
，
j
）］）；
}
else
if
（
i
>=
pdh
-
hm
&&
j
<=
hm
）
{
idx
=
bnum
-
wnum
；
im_pad
。
at
<
uchar
>
（
i
，
j
）
=
（
uchar
）（
table
［
idx
］［
im_pad
。
at
<
uchar
>
（
i
，
j
）］）；
}
else
if
（
i
>=
pdh
-
hm
&&
j
>=
pdw
-
hm
）
{
idx
=
bnum
-
1
；
im_pad
。
at
<
uchar
>
（
i
，
j
）
=
（
uchar
）（
table
［
idx
］［
im_pad
。
at
<
uchar
>
（
i
，
j
）］）；
}
//four edges except coners —— linear interpolation
else
if
（
i
<=
hm
）
{
// hbi = 0；
wbi
=
（
j
-
hm
）
/
m
；
bidx
［
0
］
=
wbi
；
bidx
［
1
］
=
bidx
［
0
］
+
1
；
p
=
（
float
）（
j
-
（
wbi
*
m
+
hm
））
/
m
；
q
=
1
-
p
；
im_pad
。
at
<
uchar
>
（
i
，
j
）
=
（
uchar
）（
q
*
table
［
bidx
［
0
］］［
im_pad
。
at
<
uchar
>
（
i
，
j
）］
+
p
*
table
［
bidx
［
1
］］［
im_pad
。
at
<
uchar
>
（
i
，
j
）］）；
}
else
if
（
i
>=
（（
hnum
-
1
）
*
m
+
hm
））
{
hbi
=
hnum
-
1
；
wbi
=
（
j
-
hm
）
/
m
；
bidx
［
0
］
=
hbi
*
wnum
+
wbi
；
bidx
［
1
］
=
bidx
［
0
］
+
1
；
float
p
=
（
float
）（
j
-
（
wbi
*
m
+
hm
））
/
m
；
float
q
=
1
-
p
；
im_pad
。
at
<
uchar
>
（
i
，
j
）
=
（
uchar
）（
q
*
table
［
bidx
［
0
］］［
im_pad
。
at
<
uchar
>
（
i
，
j
）］
+
p
*
table
［
bidx
［
1
］］［
im_pad
。
at
<
uchar
>
（
i
，
j
）］）；
}
else
if
（
j
<=
hm
）
{
hbi
=
（
i
-
hm
）
/
m
；
//wbi = 0；
bidx
［
0
］
=
hbi
*
wnum
；
bidx
［
1
］
=
bidx
［
0
］
+
wnum
；
p
=
（
float
）（
i
-
（
hbi
*
m
+
hm
））
/
m
；
q
=
1
-
p
；
im_pad
。
at
<
uchar
>
（
i
，
j
）
=
（
uchar
）（
q
*
table
［
bidx
［
0
］］［
im_pad
。
at
<
uchar
>
（
i
，
j
）］
+
p
*
table
［
bidx
［
1
］］［
im_pad
。
at
<
uchar
>
（
i
，
j
）］）；
}
else
if
（
j
>=
（（
wnum
-
1
）
*
m
+
hm
））
{
hbi
=
（
i
-
hm
）
/
m
；
wbi
=
wnum
-
1
；
bidx
［
0
］
=
hbi
*
wnum
+
wbi
；
bidx
［
1
］
=
bidx
［
0
］
+
wnum
；
p
=
（
float
）（
i
-
（
hbi
*
m
+
hm
））
/
m
；
q
=
1
-
p
；
im_pad
。
at
<
uchar
>
（
i
，
j
）
=
（
uchar
）（
q
*
table
［
bidx
［
0
］］［
im_pad
。
at
<
uchar
>
（
i
，
j
）］
+
p
*
table
［
bidx
［
1
］］［
im_pad
。
at
<
uchar
>
（
i
，
j
）］）；
}
// Double linear interpolation
else
{
hbi
=
（
i
-
hm
）
/
m
；
wbi
=
（
j
-
hm
）
/
m
；
bidx
［
0
］
=
hbi
*
wnum
+
wbi
；
bidx
［
1
］
=
bidx
［
0
］
+
1
；
bidx
［
2
］
=
bidx
［
0
］
+
wnum
；
bidx
［
3
］
=
bidx
［
1
］
+
wnum
；
p
=
（
float
）（
i
-
（
hbi
*
m
+
hm
））
/
m
；
q
=
（
float
）（
j
-
（
wbi
*
m
+
hm
））
/
m
；
im_pad
。
at
<
uchar
>
（
i
，
j
）
=
（
uchar
）（
（
1
-
p
）
*
（
1
-
q
）
*
table
［
bidx
［
0
］］［
im_pad
。
at
<
uchar
>
（
i
，
j
）］
+
（
1
-
p
）
*
q
*
table
［
bidx
［
1
］］［
im_pad
。
at
<
uchar
>
（
i
，
j
）］
+
p
*
（
1
-
q
）
*
table
［
bidx
［
2
］］［
im_pad
。
at
<
uchar
>
（
i
，
j
）］
+
p
*
q
*
table
［
bidx
［
3
］］［
im_pad
。
at
<
uchar
>
（
i
，
j
）］）；
}
im_pad
。
at
<
uchar
>
（
i
，
j
）
=
im_pad
。
at
<
uchar
>
（
i
，
j
）
+
（
im_pad
。
at
<
uchar
>
（
i
，
j
）
<<
8
）
+
（
im_pad
。
at
<
uchar
>
（
i
，
j
）
<<
16
）；
}
}
im_pad
。
rowRange
（
top
，
top
+
h
）。
colRange
（
left
，
left
+
w
）。
copyTo
（
dst
）；
delete
［］
table
；
}
結果如下：
既不過分灰暗，也沒有過分增強，效果可以說非常不錯了！
Opencv介面
其實opencv已經實現了全域性直方圖均衡化和限制對比度的自適應直方圖均衡化：
// 全域性直方圖均衡化
cv
：：
equalizeHist
（
src
，
dst
）；
// 限制對比度的自適應直方圖均衡化
// @param clip_limit： double，限制對比度的因子
// @param size： cv：：Size，影象分成網格的尺寸
cv
：：
Ptr
<
cv
：：
CLAHE
>
clahe
=
cv
：：
createCLAHE
（
clip_limit
，
size
）；
clahe
->
apply
（
src
，
dst
）；
原理和之前介紹的一致，這裡就不再細緻介紹了。
多通道問題
還有一點值得注意，直方圖均衡化通常來說都是對單通道影象進行的，如果對多通道影象的每個通道分別進行均衡化，可能會出現一些色彩被增強得過分鮮豔的情況。比如：
對R、G、B個通道分別進行均衡化：
只對明度（HSV色彩空間下的V通道）進行均衡化，色彩則不回被過分改變：
文章所用影象截自B站《凡人修仙傳》，侵刪。
參考：
《計算機視覺：演算法與應用》
CLAHE的實現和研究 - jsxyhelu - 部落格園
