圖像處理之直方圖均衡化拉伸

目錄

• 1. OpenCV實現
• 2. 原理
  • 1) 機率密度函數
  • 2) 機率分佈函數
  • 3) 原理應用
  • 4) 原理推導
• 3. 具體實現
• 4. 參考文獻

1. OpenCV實現

在OpenCV中，實現直方圖均衡化比較簡單，調用equalizeHist函數便可。具體代碼以下：

#include <iostream>
#include <opencv2\opencv.hpp>

using namespace std;
using namespace cv;

int main()
{
    Mat srcImage;
    srcImage = imread("D:\\Data\\imgDemo\\lena512color.bmp", IMREAD_GRAYSCALE);
    imshow("原圖像", srcImage);
    
    Mat dstImage;
    equalizeHist(srcImage, dstImage);
    imshow("均衡後", dstImage);

    waitKey();

    return 0;
}

注意equalizeHist函數處理的是8位單波段的mat。運行結果以下所示，能夠發現通過直方圖均衡化以後，圖像的對比度加強了不少。

2. 原理

直方圖均衡化的基本思想是把原始圖的直方圖儘量的均勻分佈，其數學原理與數學中的機率論相關。注意，我這裏不少論述犧牲了數學的嚴密性來增強可理解性，畢竟做者只是個應用者和使用者。

1) 機率密度函數

具體到一張圖像上來講，能夠把圖像的灰度（像素值）r_i看做是隨機變量，則能夠知道圖像灰度的機率爲：

對應的，對於一個連續型的隨機變量x，若是存在函數f(x)也知足上面兩個條件：

則這個函數就是機率密度函數。
離散隨機變量的機率有具體的公式讓你理解，那麼連續隨機變量的機率密度函數具體的公式是怎麼樣的呢？這個概念其實須要下面要介紹的機率分佈函數來理解。

2) 機率分佈函數

機率分佈函數就是是機率密度函數的變上限積分：

通俗來說，機率分佈函數就是全部小於當前隨機變量的機率累加。因此，機率分佈函數也被叫作累積機率函數。
知道機率分佈函數，引用下網上相關論述^[1]就能更好的理解機率密度函數了：

3) 原理應用

直方圖均衡化變換就是一種灰度級非線性變換，設r和s分別表示變換前和變換後的灰度，且r和s都進行了歸一化的處理。則直方圖均衡化變換的公式爲：

即歸一化後，直方圖均衡化的結果s就是r的機率分佈函數。

4) 原理推導

根據機率論隨機變量的函數的分佈的相關知識，有s的機率密度函數爲

如下^[2]具體論述了其應用過程：

繼續推導，有：

其中s爲r的機率分佈函數，則：

變換後變量s的機率密度爲常數，說明其機率密度爲均勻分佈的。

3. 具體實現

根據第二節的論述，就知道直方圖均衡化的具體操做了，能夠分紅如下幾步：

1. 讀取源圖像，統計源圖像的直方圖。
2. 歸一化直方圖，統計源圖像每一個像素的機率密度值和機率分佈值。
3. 將每一個像素的機率分佈值恢復到 0 到 255 的區間，做爲目標圖像的像素。
4. 寫出目標圖像。

其具體代碼實現以下，我這裏是採用 GDAL 來讀取影像的，由於我想直接操做讀
取的內存 buf，這樣更底層一些。若是你不會使用 GDAL 也沒有關係，你只須要
知道 GDAL 讀取的是按照 RGBRGBRGB…排序的內存 buf。

#include <iostream>
#include <algorithm>
#include <gdal_priv.h>

using namespace std;

//直方圖均衡化
void GetHistAvgLut(GUIntBig* anHistogram, int HistNum, vector<uint8_t > &lut)
{
    //統計像素總的個數
    size_t sum = 0;
    for (int ci = 0; ci < HistNum; ci++)
    {
        sum = sum + anHistogram[ci];
    }

    //
    vector<double> funProbability(HistNum, 0.0); //機率密度函數
    vector<double> funProbabilityDistribution(HistNum, 0.0);    //機率分佈函數

    //計算機率分佈函數
    double dsum = (double)sum;
    double accumulation = 0;
    for (int ci = 0; ci < HistNum; ci++)
    {
        funProbability[ci] = anHistogram[ci] / dsum;
        accumulation = accumulation + funProbability[ci];
        funProbabilityDistribution[ci] = accumulation;
    }

    //歸一化的值擴展爲0~255的像素值，存到顏色映射表
    lut.resize(HistNum, 0);
    for (int ci = 0; ci < HistNum; ci++)
    {
        double value = std::min<double>(std::max<double>(255 * funProbabilityDistribution[ci], 0), 255);
        lut[ci] = (unsigned char)value;
    }
}

//計算16位的顏色映射表
bool CalImgLut(GDALDataset* img, vector<vector<uint8_t>>& lut)
{
    int bandNum = img->GetRasterCount();    //波段數
    lut.resize(bandNum);

    //
    for (int ib = 0; ib < bandNum; ib++)
    {
        //計算該通道的直方圖
        int HistNum = 256;
        GUIntBig* anHistogram = new GUIntBig[HistNum];
        int bApproxOK = FALSE;
        img->GetRasterBand(ib + 1)->GetHistogram(-0.5, 255.5, HistNum, anHistogram, TRUE, bApproxOK, NULL, NULL);

        //直方圖均衡化    
        GetHistAvgLut(anHistogram, HistNum, lut[ib]);

        //
        delete[] anHistogram;
        anHistogram = nullptr;
    }

    return true;
}


int main()
{
    //
    GDALAllRegister();          //GDAL全部操做都須要先註冊格式
    CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");  //支持中文路徑

    //讀取
    const char* imgPath = "D:\\Data\\imgDemo\\lena512color.bmp";
    GDALDataset* img = (GDALDataset *)GDALOpen(imgPath, GA_ReadOnly);
    if (!img)
    {
        cout << "Can't Open Image!" << endl;
        return 1;
    }

    //
    int imgWidth = img->GetRasterXSize();   //圖像寬度
    int imgHeight = img->GetRasterYSize();  //圖像高度
    int bandNum = img->GetRasterCount();    //波段數
    int depth = GDALGetDataTypeSize(img->GetRasterBand(1)->GetRasterDataType()) / 8;    //圖像深度

    //建立顏色映射表
    vector<vector<uint8_t>> lut;
    CalImgLut(img, lut);

    //建立
    GDALDriver *pDriver = GetGDALDriverManager()->GetDriverByName("BMP"); //圖像驅動
    char** ppszOptions = NULL;
    const char* dstPath = "D:\\Data\\imgDemo\\dst.bmp";
    int bufWidth = imgWidth;
    int bufHeight = imgHeight;
    GDALDataset* dst = pDriver->Create(dstPath, bufWidth, bufHeight, bandNum, GDT_Byte, ppszOptions);
    if (!dst)
    {
        printf("Can't Write Image!");
        return false;
    }

    //讀取buf
    size_t imgBufNum = (size_t)bufWidth * bufHeight * bandNum * depth;
    GByte *imgBuf = new GByte[imgBufNum];
    img->RasterIO(GF_Read, 0, 0, bufWidth, bufHeight, imgBuf, bufWidth, bufHeight,
        GDT_Byte, bandNum, nullptr, bandNum*depth, bufWidth*bandNum*depth, depth);

    //迭代經過顏色映射表替換值
    for (int yi = 0; yi < bufHeight; yi++)
    {
        for (int xi = 0; xi < bufWidth; xi++)
        {
            for (int bi = 0; bi < bandNum; bi++)
            {
                size_t m = (size_t)bufWidth * bandNum * yi + bandNum * xi + bi;
                imgBuf[m] = lut[bi][imgBuf[m]];
            }
        }
    }
    
    //寫入
    dst->RasterIO(GF_Write, 0, 0, bufWidth, bufHeight, imgBuf, bufWidth, bufHeight,
        GDT_Byte, bandNum, nullptr, bandNum*depth, bufWidth*bandNum*depth, depth);

    //釋放
    delete[] imgBuf;
    imgBuf = nullptr;
    GDALClose(dst);
    dst = nullptr;
    GDALClose(img);
    img = nullptr;

    return 0;
}

能夠看到我這裏統計了0到255的直方圖以後，歸一化計算每一個像素的分佈機率，再還原成0到255的值並預先生成了一個顏色映射表，最後直接經過這個顏色映射表進行灰度變換。這是圖像處理的一種加速辦法。最終獲得的結果對比：

其直方圖對比：

4. 參考文獻

[1] 應該如何理解機率分佈函數和機率密度函數
[2] 直方圖均衡化的數學原理
[3] 理解機率密度函數
[4] 直方圖均衡化的數學原理
[5] 直方圖均衡化（Histogram equalization）與直方圖規定化