圖像處理基礎(8)：圖像的灰度直方圖、直方圖均衡化、直方圖規定化（匹配）

本文主要介紹了灰度直方圖相關的處理，包括如下幾個方面的內容：

• 利用OpenCV計算圖像的灰度直方圖，並繪製直方圖曲線
• 直方圖均衡化的原理及實現
• 直方圖規定化（匹配）的原理及實現

圖像的灰度直方圖

一幅圖像由不一樣灰度值的像素組成，圖像中灰度的分佈狀況是該圖像的一個重要特徵。圖像的灰度直方圖就描述了圖像中灰度分佈狀況，可以很直觀的展現出圖像中各個灰度級所佔的多少。
圖像的灰度直方圖是灰度級的函數，描述的是圖像中具備該灰度級的像素的個數：其中，橫座標是灰度級，縱座標是該灰度級出現的頻率。

不過一般會將縱座標歸一化到\([0,1]\)區間內，也就是將灰度級出現的頻率（像素個數）除以圖像中像素的總數。灰度直方圖的計算公式以下：
\[ p(r_k) = \frac{n_k}{MN} \]
其中，\(r_k\)是像素的灰度級，\(n_k\)是具備灰度\(r_k\)的像素的個數，\(MN\)是圖像中總的像素個數。

OpenCV灰度直方圖的計算

直方圖的計算是很簡單的，無非是遍歷圖像的像素，統計每一個灰度級的個數。在OpenCV中封裝了直方圖的計算函數calcHist，爲了更爲通用該函數的參數有些複雜，其聲明以下：

void calcHist( const Mat* images, int nimages,
                          const int* channels, InputArray mask,
                          OutputArray hist, int dims, const int* histSize,
                          const float** ranges, bool uniform = true, bool accumulate = false );

該函數可以同時計算多個圖像，多個通道，不一樣灰度範圍的灰度直方圖.
其參數以下：

• images，輸入圖像的數組，這些圖像要有相同大大小，相同的深度（CV_8U CV_16U CV_32F）.
• nimages ，輸入圖像的個數
• channels，要計算直方圖的通道個數。
• mask，可選的掩碼，不使用時可設爲空。要和輸入圖像具備相同的大小，在進行直方圖計算的時候，只會統計該掩碼不爲0的對應像素
• hist，輸出的直方圖
• dims，直方圖的維度
• histSize，直方圖每一個維度的大小
• ranges，直方圖每一個維度要統計的灰度級的範圍
• uniform，是否進行歸一化，默認爲true
• accumulate，累積標誌，默認值爲false。

爲了計算的靈活性和通用性，OpenCV的灰度直方圖提供了較多的參數，但對於只是簡單的計算一幅灰度圖的直方圖的話，又顯得較爲累贅。這裏對calcHist進行一次封裝，可以方便的獲得一幅灰度圖直方圖。

class Histogram1D
{
private:
    int histSize[1]; // 項的數量
    float hranges[2]; // 統計像素的最大值和最小值
    const float* ranges[1];
    int channels[1]; // 僅計算一個通道

public:
    Histogram1D()
    {
        // 準備1D直方圖的參數
        histSize[0] = 256;
        hranges[0] = 0.0f;
        hranges[1] = 255.0f;
        ranges[0] = hranges;
        channels[0] = 0;
    }

    MatND getHistogram(const Mat &image)
    {
        MatND hist;
        // 計算直方圖
        calcHist(&image ,// 要計算圖像的
            1,                // 只計算一幅圖像的直方圖
            channels,        // 通道數量
            Mat(),            // 不使用掩碼
            hist,            // 存放直方圖
            1,                // 1D直方圖
            histSize,        // 統計的灰度的個數
            ranges);        // 灰度值的範圍
        return hist;
    }

    Mat getHistogramImage(const Mat &image)
    {
        MatND hist = getHistogram(image);

        // 最大值，最小值
        double maxVal = 0.0f;
        double minVal = 0.0f;

        minMaxLoc(hist, &minVal, &maxVal);

        //顯示直方圖的圖像
        Mat histImg(histSize[0], histSize[0], CV_8U, Scalar(255));

        // 設置最高點爲nbins的90%
        int hpt = static_cast<int>(0.9 * histSize[0]);
        //每一個條目繪製一條垂直線
        for (int h = 0; h < histSize[0]; h++)
        {
            float binVal = hist.at<float>(h);
            int intensity = static_cast<int>(binVal * hpt / maxVal);
            // 兩點之間繪製一條直線
            line(histImg, Point(h, histSize[0]), Point(h, histSize[0] - intensity), Scalar::all(0));
        }
        return histImg;
    }
};

Histogram1D提供了兩個方法：getHistogram返回統計直方圖的數組，默認計算的灰度範圍是[0,255]；getHistogramImage將圖像的直方圖以線條的形式畫出來，並返回包含直方圖的圖像。測試代碼以下：

Histogram1D hist;
    Mat histImg;
    histImg = hist.getHistogramImage(image);

    imshow("Image", image);
    imshow("Histogram", histImg);

其結果以下：

直方圖均衡化 Histogram Equalization

假如圖像的灰度分佈不均勻，其灰度分佈集中在較窄的範圍內，使圖像的細節不夠清晰，對比度較低。一般採用直方圖均衡化及直方圖規定化兩種變換，使圖像的灰度範圍拉開或使灰度均勻分佈，從而增大反差，使圖像細節清晰，以達到加強的目的。
直方圖均衡化，對圖像進行非線性拉伸，從新分配圖像的灰度值，使必定範圍內圖像的灰度值大體相等。這樣，原來直方圖中間的峯值部分對比度獲得加強，而兩側的谷底部分對比度下降，輸出圖像的直方圖是一個較爲平坦的直方圖。

均衡化算法

直方圖的均衡化實際也是一種灰度的變換過程，將當前的灰度分佈經過一個變換函數，變換爲範圍更寬、灰度分佈更均勻的圖像。也就是將原圖像的直方圖修改成在整個灰度區間內大體均勻分佈，所以擴大了圖像的動態範圍，加強圖像的對比度。一般均衡化選擇的變換函數是灰度的累積機率，直方圖均衡化算法的步驟：

• 計算原圖像的灰度直方圖 \(P(S_k) = \frac{n_k}{n}\)，其中\(n\)爲像素總數，\(n_k\)爲灰度級\(S_k\)的像素個數
• 計算原始圖像的累積直方圖 \(CDF(S_k) = \sum\limits^k_{i=0}\frac{n_i}{n}=\sum\limits^k_{i=0}P_s(S_i)\)
• \(D_j = L\cdot CDF(S_i)\)，其中 \(D_j\)是目的圖像的像素，\(CDF(S_i)\)是源圖像灰度爲i的累積分佈，L是圖像中最大灰度級（灰度圖爲255）

其代碼實現以下：

• 在上面中封裝了求灰度直方圖的類，這裏直接應用該方法獲得圖像的灰度直方圖；
• 將灰度直方圖進行歸一化，計算灰度的累積機率；
• 建立灰度變化的查找表
• 應用查找表，將原圖像變換爲灰度均衡的圖像

具體代碼以下：

void equalization_self(const Mat &src, Mat &dst)
{
    Histogram1D hist1D;
    MatND hist = hist1D.getHistogram(src);

    hist /= (src.rows * src.cols); // 對獲得的灰度直方圖進行歸一化
    float cdf[256] = { 0 }; // 灰度的累積機率
    Mat lut(1, 256, CV_8U); // 灰度變換的查找表
    for (int i = 0; i < 256; i++)
    {
        // 計算灰度級的累積機率
        if (i == 0)
            cdf[i] = hist.at<float>(i);
        else
            cdf[i] = cdf[i - 1] + hist.at<float>(i);

        lut.at<uchar>(i) = static_cast<uchar>(255 * cdf[i]); // 建立灰度的查找表
    }

    LUT(src, lut, dst); // 應用查找表，進行灰度變化，獲得均衡化後的圖像

}

上面代碼只是加深下對均衡化算法流程的理解，實際在OpenCV中也提供了灰度均衡化的函數equalizeHist,該函數的使用很簡單，只有兩個參數：輸入圖像，輸出圖像。下圖爲，上述代碼計算獲得的均衡化結果和調用equalizeHist的結果對比

最左邊爲原圖像，中間爲OpenCV封裝函數的結果，右邊爲上面代碼獲得的結果。

直方圖規定化

從上面能夠看出，直方圖的均衡化自動的肯定了變換函數，能夠很方便的獲得變換後的圖像，可是在有些應用中這種自動的加強並非最好的方法。有時候，須要圖像具備某一特定的直方圖形狀（也就是灰度分佈），而不是均勻分佈的直方圖，這時候可使用直方圖規定化。
直方圖規定化，也叫作直方圖匹配，用於將圖像變換爲某一特定的灰度分佈，也就是其目的的灰度直方圖是已知的。這其實和均衡化很相似，均衡化後的灰度直方圖也是已知的，是一個均勻分佈的直方圖；而規定化後的直方圖能夠隨意的指定，也就是在執行規定化操做時，首先要知道變換後的灰度直方圖，這樣才能肯定變換函數。規定化操做可以有目的的加強某個灰度區間，相比於，均衡化操做，規定化多了一個輸入，可是其變換後的結果也更靈活。

在理解了上述的均衡化過程後，直方圖的規定化也較爲簡單。能夠利用均衡化後的直方圖做爲一箇中間過程，而後求取規定化的變換函數。具體步驟以下：

• 將原始圖像的灰度直方圖進行均衡化，獲得一個變換函數\(s = T(r)\)，其中s是均衡化後的像素，r是原始像素
• 對規定的直方圖進行均衡化，獲得一個變換函數\(v = G(z)\)，其中v是均衡化後的像素，z是規定化的像素
• 上面都是對同一圖像的均衡化，其結果應該是相等的,\(s = v,且 z = G^{-1}(v) = G^{-1}(T(r))\)

經過，均衡化做爲中間結果，將獲得原始像素\(r\)和\(z\)規定化後像素之間的映射關係。

詳解規定化過程

對圖像進行直方圖規定化操做，原始圖像的直方圖和以及規定化後的直方圖是已知的。假設\(P_r(r)\)表示原始圖像的灰度機率密度，\(P_z(z)\)表示規定化圖像的灰度機率密度（r和z分別是原始圖像的灰度級，規定化後圖像的灰度級）。

• 對原始圖像進行均衡化操做，則有\(s_k = T(r_k) = L \cdot \sum\limits_{i=0}^{i=k}P_r(r_k)\)
• 對規定化的直方圖進行均衡化操做，則\(v_k = G(z_m) = L \cdot \sum\limits_{j=0}^{j=m}P_z(z_m)\)
• 因爲是對同一圖像的均衡化操做，因此有\(s_k = v_m\)。
• 規定化操做的目的就是找到原始圖像的像素\(s_k\)到規定化後圖像像素的\(z_k\)之間的一個映射。有了上一步的等式後，能夠獲得\(s_k = G(z_k)\)，所以要想找到\(s_k\)想對應的\(z_k\)只須要在\(z\)進行迭代，找到使式子\(G(z_m)-s_k\)的絕對值最小便可。
• 上述描述只是理論的推導過程，在實際的計算過程當中，不須要作兩次的均衡化操做，具體的推導過程以下：\[ \begin{array}{c} s_k = v_k \\ L \cdot \sum\limits_{i=0}^{i=k}P_r(r_k) = L \cdot \sum\limits_{j=0}^{j=m}P_z(z_m) \\ \sum\limits_{i=0}^{i=k}P_r(r_k) = \sum\limits_{j=0}^{j=m}P_z(z_m) \end{array} \]
  上面公式表示，假如\(s_k\) 規定化後的對應灰度是\(z_m\)的話，須要知足的條件是\(s_k\)的累積機率和\(z_m\)的累積機率是最接近的。
  下面是一個具體計算的例子：
  

首先獲得原直方圖的各個灰度級的累積機率\(V_s\)以及規定化後直方圖的各個灰度級的累積機率\(V_z\)，那麼肯定\(s_k\)到\(z_m\)之間映射關係的條件就是：\[\mid V_s - V_z \mid\]的值最小。
以\(k = 2\)爲例，其原始直方圖的累積機率是:0.65,在規定化後的直方圖的累積機率中和0.65最接近（相等）的是灰度值爲5的累積機率密度，則能夠獲得原始圖像中的灰度級2，在規定化後的圖像中的灰度級是5。

直方圖規定化的實現

直方圖規定化的實現能夠分爲一下三步：

• 計算原圖像的累積直方圖
• 計算規定直方圖的累積直方圖
• 計算兩累積直方圖的差值的絕對值
• 根據累積直方圖差值創建灰度級的映射

具體代碼實現以下：

void hist_specify(const Mat &src, const Mat &dst,Mat &result)
{
    Histogram1D hist1D;
    MatND src_hist = hist1D.getHistogram(src);
    MatND dst_hist = hist1D.getHistogram(dst);

    float src_cdf[256] = { 0 };
    float dst_cdf[256] = { 0 };

    // 源圖像和目標圖像的大小不同，要將獲得的直方圖進行歸一化處理
    src_hist /= (src.rows * src.cols);
    dst_hist /= (dst.rows * dst.cols);

    // 計算原始直方圖和規定直方圖的累積機率
    for (int i = 0; i < 256; i++)
    {
        if (i == 0)
        {
            src_cdf[i] = src_hist.at<float>(i);
            dst_cdf[i] = dst_hist.at<float>(i);
        }
        else
        {
            src_cdf[i] = src_cdf[i - 1] + src_hist.at<float>(i);
            dst_cdf[i] = dst_cdf[i - 1] + dst_hist.at<float>(i);
        }
    }

    // 累積機率的差值
    float diff_cdf[256][256];
    for (int i = 0; i < 256; i++)
        for (int j = 0; j < 256; j++)
            diff_cdf[i][j] = fabs(src_cdf[i] - dst_cdf[j]);

    // 構建灰度級映射表
    Mat lut(1, 256, CV_8U);
    for (int i = 0; i < 256; i++)
    {
        // 查找源灰度級爲ｉ的映射灰度
        //　和ｉ的累積機率差值最小的規定化灰度
        float min = diff_cdf[i][0];
        int index = 0;
        for (int j = 1; j < 256; j++)
        {
            if (min > diff_cdf[i][j])
            {
                min = diff_cdf[i][j];
                index = j;
            }
        }
        lut.at<uchar>(i) = static_cast<uchar>(index);
    }

    // 應用查找表，作直方圖規定化
    LUT(src, lut, result);
}

上面函數的第二個參數的直方圖就是規定化的直方圖。代碼比較簡單，這裏就不一一解釋了。其結果以下：

左邊是原圖像，右邊是規定化的圖像，也就是上面函數的第一個和第二個輸入參數。原圖像規定化的結果以下：

原圖像規定化後的直方圖和規定化的圖像的直方圖的形狀比較相似， 而且原圖像規定化後整幅圖像的特徵和規定化的圖像也比較相似，例如：原圖像牀上的被子，明顯帶有規定化圖像中水的波紋特徵。

直方圖規定化過程當中，在作灰度映射的時候，有兩種經常使用的方法：

• 單映射 Single Mapping Law,SML，這種方法也是上面使用的方法，根據累積直方圖的差值，從原圖像中找到其在規定化圖像中的映射。
• 組映射 Group Mapping Law,GML 這種方法較上述方法複雜很多，可是處理效果較好。

對於GML的映射方法，一直沒有很好的理解，可是根據其算法描述實現了該方法，代碼這裏先不放出，其處理結果以下：

其結果較SML來講更爲亮一些，牀上的波浪特徵也更爲明顯，可是其直方圖形狀，和規定化的直方圖對比，第一個峯不是很明顯。

總結

• 圖像的灰度直方圖可以很直觀的展現圖像中灰度級的總體分佈狀況，對圖像的後續處理有很好的指導做用。
• 直方圖的均衡化的是將一幅圖像的直方圖變平，使各個灰度級的趨於均勻分佈，這樣可以很好的加強圖像對比度。直方圖均衡化是一種自動化的變換，僅須要輸入圖像，就可以肯定圖像的變換函數。可是直方圖的均衡化操做也有必定的肯定，在均衡化的過程當中對圖像中的數據不加選擇，這樣有可能會加強圖像的背景；變換後圖像的灰度級減小，有可能形成某些細節的消失；會壓縮圖像直方圖中的高峯，形成處理後圖像對比度的不天然等。
• 直方圖規定化，也稱爲直方圖匹配，通過規定化處理將原圖像的直方圖變換爲特定形狀的直方圖（上面中的示例，就是將圖像的直方圖變換爲另外一幅圖像的直方圖）。它能夠按照預先設定的某個形狀來調整圖像的直方圖,運用均衡化原理的基礎上，經過創建原始圖像和指望圖像之間的關係，選擇地控制直方圖，使原始圖像的直方圖變成規定的形狀它能夠按照預先設定的某個形狀來調整圖像的直方圖。直方圖規定化是在運用均衡化原理的基礎上，經過創建原始圖像和指望圖像之間的關係，選擇地控制直方圖，使原始圖像的直方圖變成規定的形狀，從而彌補直方圖均衡化的一些缺點.