图像直方图均衡(代码实现)

直方图均衡

直方图均衡是一种经典的图像增强方法,是灰度变换的一个重要应用,主要作用于提高图像的对比度,或者说是提高图像的动态范围,一幅对比度很低的图像的直方图范围很窄,而且灰度级变换很大,直方图均衡就是使图像经过一些特殊的变换使图像灰度级的概率分布趋近于均匀分布,最直观的表现就是经过处理后的图像直方图的灰度级范围变广,而且频率接近相等,理想情况是图像灰度级的分布是连续的,最终处理后的图像中的灰度的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;
}

运行结果:

求10000的阶乘

大数问题

这学期开了算法课,老师课上出了一个问题,就是求10000的阶乘。

说实话,第一次听到这个问题我很快就想到了去年实验室纳新时的一道笔试题

那道题是这么说的,写一个函数求两个很大很大很大的整数的和(感觉真的暴力)。

当时我大一,c语言刚自己看完,完全没有思路,甚至脑子里面都没有常用数据类型的范围这种概念,后来问了下出题的学长,他就只告诉我把这两个数存到数组里面,然后处理一下进位就行了。

他说的很轻松…

后来就查了一下,发现网上有很多答案,有用BigInteger的,还有用python直接算的…

具体思路

很多语言的int类型的最大值都是2147483647,当然也有什么long类型,long long类型,unsigned long long类型。

这些类型能存的数一个比一个大,但是都存不下这个10000的阶乘。

但是如果你问10000的阶乘有几位,就未必有这么大了,网上也有大佬推导出了这个阶乘的位数,可惜我没有看懂…

我们可以定义一个数组,把我们计算得到的数一位一位的存到数组中。在计算的过程中处理好进位就好了。

假设已经得到的结果是5040,这是7的阶乘,在数组中实际存储的是:

0 4 0 5

下面它要乘8了,按照我们正常乘法的顺序,先用8和该数的个位相乘,这里也正好是数组中的第一位,结果是0,这时还将0写到数组第一位即可。

0

用8乘第二位4时,结果是32,这时候就需要考虑进位了,我们先将32的个位2取出来写到数组第二个位置上,然后将32的最高位3保存起来。

0 2

然后用8乘第三位,结果是0,再加上上一步的进位3,写到数组第三个位置上。

0 2 3

用8乘最后一位,结果是40,同样的方法,先将个位0写在第四位,剩下进位4,这时的数组就要扩充一位来放这个最后的进位4。

0 2 3 0 4

但是这里要考虑全面,如果最后剩下的进位不是一位数,而是很多很多位数,这时数组扩充的位数就不是1了,这时可以用一个循环将最后的进位拆开,同时扩充数组,直到没有进位。

很明显,这里的计算得到的结果是倒序存储的,所以输出最后结果的时候只需要倒序输出一下数组就行了。

代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
#include <iostream>
#define RESULTLEN 50000

using namespace std;


int main()
{
int result[RESULTLEN];//结果数组,50000位
int temp = 0;//每步计算后的临时结果
int carry = 0;//进位
int result_len = 1;//结果的位数

result[0] = 1;

for(int i = 1; i <= 10000; i++)//每次与结果数组相乘的数
{
carry = 0; //每次与一个数相乘过后进位都清零
for(int j = 0; j < result_len; j++)//循环结果数组的每一位,但是是倒序循环
{
temp = carry + result[j] * i; //临时结果就是结果数组中的第j位与i的乘积加上上一次的进位
result[j] = temp % 10; //当前位的数值成为临时结果的个位数
carry = temp / 10; //计算进位,注意这个进位到后面是一个相当大的数字
}

while(carry) //将进位拆开写到对应的位置上
{
result[result_len++] = carry % 10;
carry = carry / 10;
}
}

//输出
for(int i = result_len - 1; i >= 0; i--)
{
cout << result[i];
}

cout << endl ;
cout << "10000!的长度为:"<< result_len << endl;

}

结果:

新的开始

NEW LOAD

我是鲍骞月,我现在已经是大二下学期了,之前的blog因为电脑的原因,被我删掉了。

今天重构完成…

新的目标

这是我在Android实验室的第二年了,但是Android开发水平还没有到达我想要的水平。

算了,目标太多了…

  • 硬磕自定义View
  • 进行图像处理和计算机图形学的学习
  • 尽快涉及游戏开发(我坚信这是我的道路)
  • 提高本学期成绩
  • 和几个实验室成员开发维护一个中北校园用app(类似于之前实验室前辈写的爱中北),我认为这个app应该是由我们实验室开发并维护的。
  • 在博客上提交几篇高质量的文章

加油吧!