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

詳細介紹了圖像進行中直方圖均衡的原理，并借助OpenCV和GDAL，通過C/C++來直接操作記憶體buf來底層實作。

目錄

• ​​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) 機率密度函數

具體到一張圖像上來說，可以把圖像的灰階（像素值）ri看作是随機變量，則可以知道圖像灰階的機率為：

對應的，對于一個連續型的随機變量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的值并預先生成了一個顔色映射表，最後直接通過這個顔色映射表進行灰階變換。這是圖像處理的一種加速辦法。最終得到的結果對比：

其直方圖對比：