OpenCV之直方图均衡化

1. 本文代码使用OpenCV版本：2.4.13
2. 本文代码在Win10+Visual Studio 2013 Update 3下测试通过

上两篇博客《OpenCV之直方图拉伸》和《OpenCV之查找表与直方图拉伸》讲述了拉伸图像的直方图以增强对比度。然而，在多数情况下，图像在视觉上的缺陷并非源于使用过窄的强度范围，而是由于某些颜色值出现的频率高于另一些。事实上，我们可以认为一幅高质量的图像应该平均使用所有的像素强度，这便是直方图均衡化(Histogram Equalization)背后的理念，即使得图像的直方图尽可能平坦。可以用以下图像表示：

图像对比度增强的方法可以分为两类，一类是直接对比度增强方法，另一类是间接对比度增强方法。直方图拉伸和直方图均衡化是两种最常见的间接对比度增强方法。

OpenCV实现

OpenCV提供了直方图均衡化的接口：

void equalizeHist(InputArray src, OutputArray dst);

其使用方法如下：

/////////////////////////////////////// OpenCV实现直方图均衡化示例#include "opencv2/core/core.hpp"#include "opencv2/highgui/highgui.hpp"#include "opencv2/imgproc/imgproc.hpp"int main(){      // 读入图像，现在是3通道的RGB图像    cv::Mat image = cv::imread("G:/dataset/lena512.bmp");    if (image.empty())    {        return -1;    }    // 转换为单通道的灰度图像    cv::Mat grayImage;    cv::cvtColor(image, grayImage, cv::COLOR_BGR2GRAY);    // 直方图均衡化    cv::Mat resultImage;    cv::equalizeHist(grayImage, resultImage);    // 显示原始图像及经过直方图均衡化的图像    cv::imshow("Original image", grayImage);    cv::imshow("Result image", resultImage);    cv::waitKey();    return 0;}

原始图像及其直方图如下所示：

经过直方图均衡化的图像及其直方图如下所示：

从图像来看，对比度确实增强了；从直方图来看，确实使用了更多的灰度级，且相对来说更加平坦了。需要注意的是，由于灰度的离散性，直方图均衡化只能得到相对平坦的直方图。

数学原理

直方图均衡化事实上是对灰度的重映射，将原本的灰度值映射为新的灰度值。要得到此映射函数，必须先从连续灰度说起。

连续灰度

考虑连续灰度值，用r代表原始图像的灰度，它的取值范围为[0, L-1]，其中0代表黑白，L-1代表白色。用s代表经过直方图均衡化的图像的灰度，也就是r经过映射或变换后的灰度，则有变换函数：

s=T(r),0≤r≤L−1

考虑这个变换。因为我们最终要使用所有的灰度级，所以变换后的取值范围应该为[0, L-1]。同时，为了保证变换不会产生灰度颠倒的情况，较大的灰度r必须映射为较大的灰度z，即变换T(r)必须是单调递增的。

下面重述这两个条件：

(a) T(r)在区间[0, L-1]上为单调递增函数；如果需要得到T(r)的反函数，则要求其为严格单调递增函数。

(b) 当0≤r≤L−1时，0≤T(r)≤L−1。

连续随机变量

为了进一步得到函数关系，我们将一幅图像的灰度级看成是区间[0, L-1]内的随机变量。描述随机变换的基本描绘子是其概率密度函数(Probability Density Function， PDF)。令pr(r)和ps(s)分别表示随机变量r和s的PDF。由基本概率论的知识可以知道，如果pr(r)和T(r)已知，且在感兴趣的区域上T(r)是连续且可微的，则变换后的变量s的PDF可由下面的简单公式得到：

ps(s)=pr(r)|drds|=pr(r)|drdT(r)|

由于要均衡化的直方图是均匀分布的，所以它的PDF是已知的，即

ps(s)=1L−1

于是有：

pr(r)|drdT(r)|=1L−1

|dT(r)dr|=(L−1)pr(r)

T(r)=(L−1)ddr[∫r0pr(w)dw]

最终得到了T(r)的变换函数，它将任意的PDF转换为均匀的PDF。如下图所示：

。

可以验证，它满足条件(a)和(b)。

离散灰度

得到连续灰度的变换函数后，就可以将其应用到离散的情境下。对于离散值，我们处理其概率(直方图值)与求和来替代PDF与积分。正如前面提到的那样，一幅数学图像中灰度级rk出现的概率近似为

pr(rk)=nkMN,k=0,1,2,...,L−1

其中，MN是图像像素的总数。

而变换T(r)的离散形式为

sk=T(rk)=(L−1)∑j=0kpr(rj)=L−1MN∑j=0knj,k=0,1,2,...,L−1

这样，我们就得到了将原始图像的灰度值映射为新灰度值的映射函数。它使用了累积直方图(∑kj=0nj)。

C++实现

以下的C++代码描述了上述过程：

/////////////////////////////////////// OpenCV直方图均衡化C++代码示例#include "opencv2/core/core.hpp"#include "opencv2/highgui/highgui.hpp"#include "opencv2/imgproc/imgproc.hpp"// 计算直方图cv::Mat calcHist(const cv::Mat& image){    CV_Assert(!image.empty() && image.channels() == 1 && image.depth() == CV_8U);    cv::Mat hist;    int channels[1] = { 0 };    int histSize = 256;    float range[2] = { 0, 256 };    const float* ranges[1] = { range };    cv::calcHist(&image, 1, channels, cv::Mat(), hist, 1, &histSize, ranges);    return hist;}// 计算累积计算图(累积分布函数)cv::Mat calcCDF(const cv::Mat& hist){    CV_Assert(!hist.empty() && hist.cols == 1 && hist.depth() == CV_32F);    cv::Mat cdf(hist.size(), CV_32FC1);    cdf.at<float>(0, 0) = hist.at<float>(0, 0);    float total = cdf.at<float>(0, 0);    for (int i = 1; i < hist.rows; ++i)    {        total += hist.at<float>(i, 0);        cdf.at<float>(i, 0) = total;    }    for (int i = 0; i < cdf.rows; ++i)    {        cdf.at<float>(i, 0) = cdf.at<float>(i, 0) / total;    }    return cdf;}// 根据累积直方图计算查找表cv::Mat calcLUT(const cv::Mat& cdf){    CV_Assert(!cdf.empty() && cdf.cols == 1);    cv::Mat lut(1, cdf.rows, CV_8UC1);    for (int i = 0; i < cdf.rows; ++i)    {        lut.at<uchar>(0, i) = static_cast<uchar>(std::round(cdf.at<float>(i, 0) * 255));    }    return lut;}// 直方图均衡化cv::Mat histogramEqulization(const cv::Mat& image){    CV_Assert(!image.empty() && image.channels() == 1 && image.depth() == CV_8U);    cv::Mat lut = calcLUT(calcCDF(calcHist(image)));    cv::Mat eqImage;    cv::LUT(image, lut, eqImage);    return eqImage;}int main(){    ///////////////////////////////////////////////////////    // 直方图均衡化    // 读入图像，此时是3通道的RGB图像    cv::Mat image = cv::imread("G:/dataset/lena512.bmp");    if (image.empty())    {        return -1;    }    // 转换为单通道的灰度图像    cv::Mat grayImage;    cv::cvtColor(image, grayImage, cv::COLOR_BGR2GRAY);    // 直方图均衡化    cv::Mat imageEqualizedMy = histogramEqulization(grayImage);    // 显示图像    cv::imshow("result image my", imageEqualizedMy);    cv::waitKey();    return 0;}

经对比，此代码的结果与OpenCV实现的函数的结果完全相同，不再展示。

参考链接

1. 《OpenCV 2 视觉编程手册》4.4节
2. 《数字图像处理(第三版》3.3.1节
3. 【数字图像处理】直方图均衡化详解及编程实现

(转载请保留作者信息)