OpenCV-跟我一塊兒學數字圖像處理之直方圖均衡化

從這篇博文開始，小生正式從一個絕不相干專業轉投數字圖像處理。廢話很少說了，talk is cheap. show me the code.

• 直方圖均衡化目的

因爲一些圖像灰度的分佈過於集中，這樣會致使圖像的層次不夠分明，直方圖均衡化就是爲了讓圖像的灰度分佈更均勻，圖像的層次感更強。

• 數學原理

　　基於連續灰度分佈的結論推導

直方圖均衡化屬於數字圖像處理中灰度變換（intensity transformation）的內容，灰度變換的目的就是找到一個合適的映射函數s=T(r).將原圖像的灰度值映射到新的圖像中，已達到優化圖像的目的。

假設原圖像的灰度統計直方圖標準化後爲p_r(r).原圖像灰度範圍爲（0~L-1)。那麼直方圖均衡化找到的就是這樣一個映射函數：

設映射後的圖像的灰度分佈爲p_s(s),在由機率論相關理論（隨機變量函數的機率密度與隨機變量機率密度的關係）可知：

對映射函數兩邊進行求導

因此咱們能夠獲得變換後的圖像直方圖分佈爲

咱們能夠看到，變換後的圖像灰度直方圖分佈恆爲1/（L-1)，這就達到了上面的目的，使得圖像的灰度分佈更均勻，層次感更強。

 

注：在灰度變換中，變換函數T(r)須要知足下面的兩點要求，

1. 當0≤r≤L-1時，T(r)是一個嚴格遞增函數。
2. 當0≤r≤L-1時，0≤T(r)≤L-1。

第一點要求的緣由是，對於變化前像素和變換後像素灰度的明暗順序不能改變，之因此要嚴格遞增，是爲了確保變化前和變換後像素能夠一一對應。

第二點要求的緣由是，變換後的圖像不能超過原先的灰度級數。

不難發現，其實直方圖均衡化的過程並不必定知足條件1。因此該變換時不可逆的。

　　公式的離散化

設原圖像灰度等級爲0、一、2……L-1.離散化後的映射公式就是

在利用上面公式進行計算的時候，須要把計算的結果s(r)，近似爲最近的整數。

• 基於OpenCV的直方圖均衡化

OpenCV中有專門的直方圖均衡化函數，equalizeHist，定義的頭文件在imaproc/imaproc.hpp中。

基於OpenCV的直方圖均衡化代碼段：

 1 //load the original image and show
 2 Mat src,dst_1;
 3 src = imread("test.jpg",0);
 4 namedWindow("OriginalGrayImage");
 5 imshow("OriginalGrayImage",src);
 6 
 7 //use the OpenCV measure do histogram equalization
 8 equalizeHist(src,dst_1);
 9 namedWindow("histogram equalization_opencv");
10 imshow("histogram equalization_opencv",dst_1);

仿真結果：

原圖：

使用equalizeHist均衡化後的結果：

• 根據推導過程，本身編寫的直方圖均衡化

代碼段以下：

 1 //rewrite the histogram algorithm
 2 //get some needed information
 3 int nr = src.rows;
 4 int nc = src.cols;
 5 int n = nr*nc;
 6 Mat dst_2(nr,nc,CV_8U);
 7 
 8 //get the histogram of original image
 9 uchar *p_1 = NULL;
10 unsigned int hist[256] = {0};
11 for(int i=0;i<nr;i++)
12 {
13     p_1 = src.ptr<uchar>(i);
14     for(int j=0;j<nc;j++)
15     {
16         hist[p_1[j]] = hist[p_1[j]]+1;
17     }
18 }
19 
20 //calculate the transform function
21 uchar transf_fun[256] = {0};
22 transf_fun[0] = (uchar)(255*hist[0]/n);
23 for(int i=1;i<256;i++)
24 {
25     hist[i] = hist[i-1]+hist[i];
26     transf_fun[i] = (uchar)(255*hist[i]/n);
27 }
28    
29 //pad dst_2 the equalized values
30 uchar *p_2 = NULL;
31 for(int i=0;i<nr;i++)
32 {
33     p_2 = dst_2.ptr<uchar>(i);
34     p_1 = src.ptr<uchar>(i);
35     for(int j=0;j<nc;j++)
36     {
37          p_2[j] = transf_fun[p_1[j]];
38     }
39 }
40 
41 //show the results of our own histogram algorithm
42 namedWindow("histogram equalization_own");
43 imshow("histogram equalization_own",dst_2)       

運行結果：