直方图均衡

直方图均衡是一种经典的图像增强方法,是灰度变换的一个重要应用,主要作用于提高图像的对比度,或者说是提高图像的动态范围,一幅对比度很低的图像的直方图范围很窄,而且灰度级变换很大,直方图均衡就是使图像经过一些特殊的变换使图像灰度级的概率分布趋近于均匀分布,最直观的表现就是经过处理后的图像直方图的灰度级范围变广,而且频率接近相等,理想情况是图像灰度级的分布是连续的,最终处理后的图像中的灰度的PDF(概率密度函数)是均匀的。
简而言之,直方图均衡是通过拉伸像素强度分布范围来增强图像对比度的一种方法。

图像直方图

灰度级范围为[0,L-1]的数字图像的直方图是一个离散函数h(rk) = nk,其中rk代表某个灰度级,nk为图像中灰度值为rk的像素个数。

直方图通常有两种形式,其中一种横轴是灰度级rk,纵轴数值就是nk。

另一种是前者的归一化处理,首先计算得到该图像的总像素数,设M和N为图像的行和列的维数,所以可以计算每个灰度值在一幅图像中的概率P(rk)=rk/M*N,其中k=0,1,2..L-1,归一化直方图的所有分量之和应该等于一。

直方图

核心

理解直方图均衡首先就要把一幅图像看作是一个样本空间,然后将灰度级看成在这个样本空间中的一个随机变量,这样我们就能使用概率方法来研究问题了。

我们首先来设一个方程s=T(r),其中r为待处理图像的灰度,取值范围为[0,L-1],s就是对应的处理后的灰度值,也就是说我们通过这个T变换,就可以将原图像中每个灰度值是r的像素改成s。

还需要满足两个条件:

  • T(r)在区间0 ≤ r ≤ L-1上为单调递增函数。
  • 0 ≤ r ≤ L-1时,0 ≤ T(r) ≤ L-1

概率推导

我们设灰度级是区间[0,L-1]内的一个随机变量,我们先考虑连续的情况,在概率论中研究连续型随机变量离不开其概率密度函数(PDF),令Pr(r)和Ps(s)分别表示随机变量r和s的概率密度函数。

这里我们提出这个变换公式:

1
s = T(r) = (L-1) * ∫ Pr(w)dw

其中w是积分的假变量,等号右部分有随机变量r的概率分布函数(CDF),我们来考察一下这个变换公式是否符合上一节提出的两个条件:

  • 概率密度函数总为正,Pr的积分的几何意义就是该函数下方的面积,随着r向右增大,面积肯定增大,所以满足第一个条件
  • 右边积分最大上限是r = L-1,下限是0,所以积分最大为1(曲线下总面积为1),s的最大值是L-1,满足第二个条件

为了分析该变换的效果,我们就要把变换前后灰度值的概率分布函数分别求出来,其实就是Pr(r)和Ps(s),其中变换前的分布函数Pr(r)是已知的,根据一个基本概率定理我们可以得到下面的结果(如果没看明白,请看概率统计课本28页复习下):

1
2
3
4
Ps(s) = Pr(r) * |dr / ds|         
ds / dr = T'(r) = (L-1) * Pr(r)
由上两式得:
Ps(s) = Pr(r) * |1 / (L-1)*Pr(r)| = 1 / L-1

很明显,Ps(s)是均匀分布的。
需要指明的一点的,经过这种变换之后的Ps(s)始终是均匀的,它与Pr(r)的形式无关

具体实现

实际中的图像的灰度值大多都是离散值,对于离散值,我们处理其概率(直方图值)与求和来代替处理概率密度函数和积分。

所以上面分析过的变换函数T的离散版本就是:

离散版本
其中MN是图像中像素总数,L是图像中可能的灰度级的数量,nk是灰度为rk的像素个数。

光看公式可能不好理解,下面简单举一个书上的例子:

我们假设一幅图像大小是64x64像素(MN = 4096),灰度级深度为8(L = 8),图像灰度值取值范围为[0,7],下面是灰度分布表:

直方图

应用上面的离散型公式计算处理后的像素值:

s 计算 结果 最终结果(四舍五入)
s0 7Pr(r0) 1.33 1
s1 7Pr(r0)+7Pr(r1) 3.08 3
s2 7Pr(r0)+7Pr(r1)+7Pr(r2) 4.55 5
s3 5.67 6
s4 6.23 6
s5 6.65 7
s6 6.86 7
s7 7.00 7

可以看到最终结果只有5个不同的灰度级,有很多不同的灰度值被转换成对应一个灰度值,图像灰度值深度从8变为了5。

为了保证最终的结果图像中像素值仍为整数,所以这里对每一步求和结果进行了四舍五入。

看完上面这个具体的计算过程,我们就可以编码实现了,这里我使用C++,并使用了Opencv的一些辅助函数。

#include <opencv2/opencv.hpp>

using namespace cv;

/*
直方图均衡
*/
Mat HistogramBalance(Mat& src)
{
    Mat des(src);//为节省空间,只复制矩阵头
    int gray_freq[256] = { 0 };//统计每个灰度级出现的次数
    float gray_prob[256] = { 0 };//每个灰度级在图像中的概率
    float gray_cdf_prob[256] = { 0 };//每个灰度级经过直方图均衡变换T后的积累分布概率之和    
    int gray_des_lut[256] = { 0 };//直方图均衡后的灰度查找表
    int total_pixels = src.rows * src.cols;
    //首先统计每个灰度级在图像中出现的次数
    for (int i = 0; i < src.rows; i++)
    {
        uchar* p = src.ptr<uchar>(i); //获取到行像素指针
        for (int j = 0; j < src.cols; j++)
        {
            int pixel = p[j];
            gray_freq[pixel]++;
        }
    }

    //计算灰度级在图像中出现的概率
    for (int i = 0; i < 256; i++)
        gray_prob[i] = ((float)gray_freq[i] / total_pixels); //每个灰度级占像素总数的比例


    //计算每个灰度级经直方图均衡变换后的累计分布概率之和
    gray_cdf_prob[0] = gray_prob[0];//灰度值为0为第一个值,概率和等于它本身
    for (int i = 1; i < 256; i++)
        gray_cdf_prob[i] = gray_cdf_prob[i - 1] + gray_prob[i];

    //构建直方图均衡后的灰度级查找表,公式s = T(r) = (L-1)ΣP(rj)   
    //最终灰度值四舍五入
    for (int i = 0; i < 256; i++)
        gray_des_lut[i] = (256 - 1) * gray_cdf_prob[i] + 0.5;

    //构建变换后的图像Mat
    for (int i = 0; i < des.rows; i++)
    {
        uchar *p = des.ptr<uchar>(i);
        for (int j = 0; j < des.cols; j++)
        {
            p[j] = gray_des_lut[p[j]]; //从查找表中读取变换后的灰度值
        }
    }

    return des;
}

int main(int argc, char* argv[])
{
    Mat src = imread("C://Users//鲍骞月//Desktop//testImag//Histogramtestsrc1.jpg", CV_LOAD_IMAGE_GRAYSCALE);
    namedWindow("灰度图原图像", WINDOW_AUTOSIZE);
    imshow("灰度图原图像", src);
    Mat des = HistogramBalance(src);
    imshow("灰度图直方图均衡", des);
    Mat src_color = imread("C://Users//鲍骞月//Desktop//testImag//Histogramtestsrc1.jpg", CV_LOAD_IMAGE_COLOR);
    imshow("彩图原图像", src_color);
    Mat imgBlueChannel, imgGreenChannel, imgRedChannel, des_color; //三通道Mat
    std::vector<Mat> channels;//三通道灰度值
    split(src_color, channels);
    imgBlueChannel = channels.at(0);
    imgGreenChannel = channels.at(1);
    imgRedChannel = channels.at(2);
    //分通道进行直方图均衡处理
    imgBlueChannel = HistogramBalance(imgBlueChannel);
    imgGreenChannel = HistogramBalance(imgGreenChannel);
    imgRedChannel = HistogramBalance(imgRedChannel);
    channels.clear();
    channels.push_back(imgBlueChannel);
    channels.push_back(imgGreenChannel);
    channels.push_back(imgRedChannel);
    //合并处理后的三个通道
    merge(channels, des_color);
    imshow("彩图直方图均衡", des_color);
    waitKey(0);
    return 0;
}

这里我分别实现了灰度图像和彩色图像的直方图均衡,彩色图像的实现是直接将图像三通道分开进行均衡,然后再将三个通道合并显示(不知道这样做有没有错…)。

  • 运行结果:

res1

res2

均衡后的图像灰度级跨越了更宽的灰度级范围,并且近似分布均匀,最直观的感受就是对比度的增强。

使用API

Opencv内部封装好了直方图均衡的api,函数声明是:

CV_EXPORTS_W void equalizeHist( InputArray src, OutputArray dst );
  • 第一个参数是输入图像,要求是8位单通道图像。
  • 第二个参数是输出图像,要求和原图像有一样的尺寸和类型。
int main()
{
    //加载图像
    Mat srcImage, dstImage;
    srcImage = imread("C://Users//鲍骞月//Desktop//testImag//Histogramtestsrc1.jpg", CV_LOAD_IMAGE_GRAYSCALE);
    imshow("原图像", srcImage);
    //直方图均衡
    equalizeHist(srcImage, dstImage);
    //显示结果
    imshow("直方图均衡", dstImage);
    waitKey(0);
    return 0;
}

运行结果: