图像处理创建直方图

2012-05-25 11:05:02 xiaowei_cqu 阅读数 82050
  • 爬虫简介

    快速入门掌握Python爬虫技术,包括: 1.Python爬虫原理,入门 2.基本用法 3.爬虫应用

    2905人学习 汤小洋
    免费试看

灰度直方图是数字图像中最简单且有用的工具,这一篇主要总结OpenCV中直方图CvHistogram的结构和应用。

灰度直方图的定义

灰度直方图是灰度级的函数,描述图像中该灰度级的像素个数(或该灰度级像素出现的频率):其横坐标是灰度级,纵坐标表示图像中该灰度级出现的个数(频率)。
一维直方图的结构表示为

高维直方图可以理解为图像在每个维度上灰度级分布的直方图。常见的是二维直方图。如红-蓝直方图的两个分量分别表示红光图像的灰度值和蓝光图像灰度值的函数。其图像坐标(Dr,Db)处对应在红光图像中具有灰度级Dr同时在蓝光图像中具有灰度级Db的像素个数。这是基于多光谱——每个像素有多个变量——的数字图像,二维中对应每个像素统计个变量。

OpenCV中的直方图CvHistogram

注意我们在上面理解直方图的意义时更多把他想象成一幅“图”,继而理解图中横坐标,纵坐标的意义。而在OpenCV中,应该更多把直方图看做“数据结构”来理解。
typedef struct CvHistogram
{
    int     type;
    CvArr*  bins;  //存放每个灰度级数目的数组指针
    float   thresh[CV_MAX_DIM][2];  //均匀直方图
    float** thresh2; //非均匀直方图
    CvMatND mat;  //直方图数组的内部数据结构
}
CvHistogram;
这个结构看起来简单(比IplImage*元素少多了。。。)其实并不太好理解。
第一个成员type用来指定第二个成员bins的类型。OpenCv中常见到CvArr*的接口,可以用以指定诸如CvMat、CvMatND、IplImage的类型,其实CvArr*的是一个指向void的指针。在函数内部有时需要得到确切的指向类型,这就需要type来指定。
thresh用来指定统计直方图分布的上下界。比如[0 255]表示用来统计图像中像素分别在灰度级[0 255]区间的分布情况,CV_MAX_DIM对应直方图的维数,假如设定二维红-蓝直方图的thresh为[0 255;100 200],就是分别统计红色图像灰度级在[0 255]以及蓝色图像在灰度级[100 200]的分布情况。
thresh用以指定均匀直方图的分布,我们按每个像素理解自然是“均匀分布”,其实也可以统计像素在几个区间的分布。如果统计像素在2个区间的分布,则对应[0 255]的上下界,均匀分布统计的区间即[0 127] [127 255]分布的概率,这也是为什么thresh第二个维数默认为2——会自动均分上下界;而thresh2指定非均匀的分布,这就需要指定每个区间的上下界,如果要统计直方图在区间(0,10,100,255)的分布,那需要指定thresh2的一个维度为[0 10 100 255],所以用float**形式表示。
mat简单说就是存储了直方图的信息,即我们统计的直方图分布概率。

创建直方图 cvCreateHist()

OpenCV中用cvCreateHist()创建一个直方图:
CvHistogram* cvCreateHist( 
	int dims, //直方图维数 
	int* sizes,//直翻图维数尺寸
	int type, //直方图的表示格式
        float** ranges=NULL, //图中方块范围的数组
	int uniform=1 //归一化标识
	);
size数组的长度为dims,每个数表示分配给对应维数的bin的个数。如dims=3,则size中用[s1,s2,s3]分别指定每维bin的个数。
type有两种:CV_HIST_ARRAY 意味着直方图数据表示为多维密集数组 CvMatND; CV_HIST_TREE 意味着直方图数据表示为多维稀疏数组 CvSparseMat。
ranges就是那个复杂的不好理解的thresh的范围,他的内容取决于uniform的值。uniform为0是均匀的,非0时不均匀。
计算图像直方图的函数为CalcHist():
void cvCalcHist( 
	IplImage** image, //输入图像(也可用CvMat**)
	CvHistogram* hist, //直方图指针
                 int accumulate=0, //累计标识。如果设置,则直方图在开始时不被清零。
	const CvArr* mask=NULL //操作 mask, 确定输入图像的哪个象素被计数
	);
要注意的是这个函数用来计算一张(或多张)单通道图像的直方图,如果要计算多通道,则用这个函数分别计算图像每个单通道。

实践:一维直方图

下面实践一下用OpenCV生成图像的一维直方图
int main( )
{
	IplImage * src= cvLoadImage("baboon.jpg");
	IplImage* gray_plane = cvCreateImage(cvGetSize(src),8,1);
	cvCvtColor(src,gray_plane,CV_BGR2GRAY);

	int hist_size = 256;    //直方图尺寸
	int hist_height = 256;
	float range[] = {0,255};  //灰度级的范围
	float* ranges[]={range};
	//创建一维直方图,统计图像在[0 255]像素的均匀分布
	CvHistogram* gray_hist = cvCreateHist(1,&hist_size,CV_HIST_ARRAY,ranges,1);
	//计算灰度图像的一维直方图
	cvCalcHist(&gray_plane,gray_hist,0,0);
	//归一化直方图
	cvNormalizeHist(gray_hist,1.0);

	int scale = 2;
	//创建一张一维直方图的“图”,横坐标为灰度级,纵坐标为像素个数(*scale)
	IplImage* hist_image = cvCreateImage(cvSize(hist_size*scale,hist_height),8,3);
	cvZero(hist_image);
	//统计直方图中的最大直方块
	float max_value = 0;
	cvGetMinMaxHistValue(gray_hist, 0,&max_value,0,0);
	
	//分别将每个直方块的值绘制到图中
	for(int i=0;i<hist_size;i++)
	{
		float bin_val = cvQueryHistValue_1D(gray_hist,i); //像素i的概率
		int intensity = cvRound(bin_val*hist_height/max_value);  //要绘制的高度
		cvRectangle(hist_image,
			cvPoint(i*scale,hist_height-1),
			cvPoint((i+1)*scale - 1, hist_height - intensity),
			CV_RGB(255,255,255));  
	}
	cvNamedWindow( "GraySource", 1 );
	cvShowImage("GraySource",gray_plane);
	cvNamedWindow( "H-S Histogram", 1 );
	cvShowImage( "H-S Histogram", hist_image );

	cvWaitKey(0);
}
试验结果:

对应的,我们可以用一样的思路统计每个通道的直方图,并绘制图像每个通道像素的分布:


实践:二维直方图

我们也可以结合OpenCV的例子生成二维直方图:
IplImage* r_plane  = cvCreateImage( cvGetSize(src), 8, 1 );
	IplImage* g_plane  = cvCreateImage( cvGetSize(src), 8, 1 );
	IplImage* b_plane  = cvCreateImage( cvGetSize(src), 8, 1 );
	IplImage* planes[] = { r_plane, g_plane };
	//将HSV图像分离到不同的通道中
	cvCvtPixToPlane( src, b_plane, g_plane, r_plane, 0 );

	// 生成二维直方图数据结构
	int r_bins =256, b_bins = 256; 
	CvHistogram* hist;
	{
		int    hist_size[] = { r_bins, b_bins };
		float  r_ranges[]  = { 0, 255 };          // hue is [0,180]
		float  b_ranges[]  = { 0, 255 }; 
		float* ranges[]    = { r_ranges,b_ranges };
		hist = cvCreateHist( 2, hist_size, CV_HIST_ARRAY, ranges, 1); 
	}
	//计算一张或多张单通道图像image(s) 的直方图
	cvCalcHist( planes, hist, 0, 0 );
刚才的图我们是对应每个横坐标绘制纵坐标的直方块,二维的图需要绘制每个点:
for( int h = 0; h < r_bins; h++ ) {
  for( int s = 0; s < b_bins; s++ ) {
    float bin_val = cvQueryHistValue_2D( hist, h, s ); //查询直方块的值
        int intensity = cvRound( bin_val * 255 / max_value );
        cvRectangle( hist_img, 
          cvPoint( h*scale, s*scale ),
          cvPoint( (h+1)*scale - 1, (s+1)*scale - 1),
          CV_RGB(intensity,intensity,intensity), 
          CV_FILLED);
    }
}


最终生成二维直方图:

直方图的应用以后再讨论。

转载请注明出处:http://blog.csdn.net/xiaowei_cqu/article/details/7600666

实验代码下载:http://download.csdn.net/detail/xiaowei_cqu/4328881


Mat格式的参考这里:

http://blog.csdn.net/xiaowei_cqu/article/details/8833799



(转载请注明作者和出处:http://blog.csdn.net/xiaowei_cqu 未经允许请勿用于商业用途)







2016-11-24 00:58:35 qq_26399665 阅读数 2250
  • 爬虫简介

    快速入门掌握Python爬虫技术,包括: 1.Python爬虫原理,入门 2.基本用法 3.爬虫应用

    2905人学习 汤小洋
    免费试看

  Matlab 直方图均衡化 

  直方图均衡化又称直方图修平,是一种很重要的非线性点运算。使用该方法可以加强图像的局部对比度,尤其是当图像的有用数据的对比度相当接近的时候。通过这种方法,亮度可以更好的在直方图上分布。

  直方图均衡化的基本思想是把原始图像的直方图变换为均匀分布的形式。这样增加了灰度值的动态范围,从而达到增强图像整体对比度的效果,直方图均衡化是图像增强的一种基本方法,可提高图像的对比度,即:将较窄的图像灰度范围以一定规则拉伸至较大(整个灰度级范围内)的范围。目的是在得到在整个灰度级范围内具有均匀分布的图像。

注:对比度指的是一幅图像中明暗区域最亮的白和最暗的黑之间不同亮度层级的测量,即指一幅图像灰度反差的大小。差异范围越大代表对比越大,差异范围越小代表对比越小,好的对比率120:1就可容易地显示生动、丰富的色彩,当对比率高达300:1时,便可支持各阶的颜色。

      实行步骤:

对图像(灰度图)进行直方图均衡化主要有一下几个步骤:

1、计算各个灰度值(0-255)出现的次数,即统计直方图.

2、用统计直方图计算各个灰度值的累积分布图,就是从1到这个灰度的所有点的个数。

2、根据累积分布图s(k)归一化计算出原来各灰度值的均衡化之后的新的值。新的灰度值=255*(s(k)/nm)。

    实现代码:

Image=imread('sc.png');
% 以对话框的形式选择打开一幅图像
[m,n,nDims]=size(Image);
subplot(3,3,1);imshow(uint8(Image));title('原始图像');%显示原始图像
Image=rgb2gray(Image);
subplot(3,3,4);imshow(Image);title('灰度图像');
subplot(3,3,5);imhist(Image);title('灰度图像直方图');%显示灰度图像及直方图
r=zeros(256,1);
for i=1:m
   for j=1:n
       r(Image(i,j)+1,1)= r(Image(i,j)+1,1)+1;
    end
end  %计算灰度直方图中的数值:每个灰度级对应的像素数目。
subplot(3,3,6);stem(r);
title('计算所得灰度图像直方图');
s=zeros(256,1);
s(1,1)=r(1,1);
for i=2:256
    s(i,1)=s(i-1,1)+r(i,1);
end         %累积分布函数:对应的也是像素的数目。
subplot(3,3,9);stem(s);title('累积分布直方图');

for i=1:256
    s(i,1)=floor(255*s(i,1)/(m*n));
end   % s(i)/(height*width) 为频率,*256 为归一到0—255之间,floor为取整(整数部分)函数。【round为四舍五入函数,ceil为取整数部分加1】
subplot(3,3,8);stem(s);title('累积分布归一化直方图');

I_HE=Image;
for i=1:m
   for j=1:n
       I_HE(i,j)= s(Image(i,j)+1);
    end
end%得到均衡化后的图像。s(1)~s(256)里的数值即为灰度值,1~256标号对应的是原始灰度图像的0~255的灰度值。

k=zeros(256,1);
for i=1:m
   for j=1:n
       k(I_HE(i,j)+1)= k(I_HE(i,j)+1)+1;
    end
end  %计算直方图中的数值:每个灰度级对应的像素数目。
subplot(3,3,7);stem(k);
title('计算所得均衡化后图像直方图');

subplot(3,3,2);imshow(I_HE);                          
title('均衡化后的图像');
subplot(3,3,3);imhist(I_HE);
title('均衡化后图像直方图');


     结果:



       注解:

10 matlab中小数取整的函数大约有四个:floor、ceil、round、fix
   floor:朝负无穷方向靠近最近的整数;
   ceil:朝正无穷方向靠近最近的整数;
   round:取最近的整数(相当于四舍五入)
   fix:取离0最近的整数


11 rgb2gray() rgb图像灰度化,也就是变成单通道图像。size()可以获取图像的长、宽和波段数。


13 imhist()用来显示图像数据灰度级统计直方图。 imhist(I,n)  计算和显示图像I的直方图,n为指定的灰度级数目,默认为256。如果I是二值图像,那么n仅有两个值。


14 subplot(m,n,p)、subplot(mnp)是将多个图画到一个平面上的工具。其中,m表示是图排成m行,n表示图排成n列,也就是整个figure中有n个图是排成一行的,一共m行,p是指你现在要把曲线画到figure中哪个图上,最后一个如果是1表示是从左到右第一个位置。


15 plot 是绘制二维图形的最基本函数,它是针对向量或矩阵的列来绘制曲线的。也就是说,使用plot 函数之前,必须首先定义好曲线上每一点的x 及y 坐标,常用格式为:
   (1)plot(x) 当x 为一向量时,以x 元素的值为纵坐标,x 的序号为横坐标值绘制曲线。当x 为一实矩阵时,则以其序号为横坐标,按列绘制每列元素值相对于其序号的曲
        线,
   (2)plot(x,y) 以x 元素为横坐标值,y 元素为纵坐标值绘制曲线。
   (3)plot(x,y1,x,y2,…) 以公共的x 元素为横坐标值,以y1,y2,… 元素为纵坐标值绘制多条曲线。


16 title()在图形正上方加个标题。


17 stem()该句法表示绘制数据序列Y。
   在绘制时,Y的数据是根据等距离的并且自动生成的x轴数据扩展而来的。如果Y是向量,等距离生成一个值为相应值的stem;如果Y是一个矩阵,该句法表示每一行的x值是相同的,每一列的x值间隔是等距的。


18 数组(array)就是我们最熟悉的array,在Matlab可以建立任意尺寸和维数,只要你的内
存足够不够的时候会提示。
   向量(vector)特指1*n或n*1的数组,前者称为行向量,后者称为列向量。
   矩阵(matrix)一般特指二维数组,其它与数组相同。


19 figure 创建一个窗口。
   figure
   figure('PropertyName'.'propertyvalue')
   figure(h)
   h=figure(...)


20 s(1,:)是取矩阵s中的第一行的所有元素,:是所有列的默认;:放前面就是所有行的默认。


参考文章:http://www.cnblogs.com/diligentcalf/p/4063493.html

2019-02-10 23:08:55 charlee44 阅读数 486
  • 爬虫简介

    快速入门掌握Python爬虫技术,包括: 1.Python爬虫原理,入门 2.基本用法 3.爬虫应用

    2905人学习 汤小洋
    免费试看

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。运行结果如下所示,可以发现经过直方图均衡化之后,图像的对比度增强了很多。
图1

2. 原理

直方图均衡化的基本思想是把原始图的直方图尽可能的均匀分布,其数学原理与数学中的概率论相关。注意,我这里很多论述牺牲了数学的严密性来加强可理解性,毕竟作者只是个应用者和使用者。

1) 概率密度函数

具体到一张图像上来说,可以把图像的灰度(像素值)ri看作是随机变量,则可以知道图像灰度的概率为:
f1
对应的,对于一个连续型的随机变量x,如果存在函数f(x)也满足上面两个条件:
f2
则这个函数就是概率密度函数。
离散随机变量的概率有具体的公式让你理解,那么连续随机变量的概率密度函数具体的公式是怎么样的呢?这个概念其实需要下面要介绍的概率分布函数来理解。

2) 概率分布函数

概率分布函数就是是概率密度函数的变上限积分:
f3
通俗来讲,概率分布函数就是所有小于当前随机变量的概率累加。所以,概率分布函数也被叫做累积概率函数。
知道概率分布函数,引用下网上相关论述[1]就能更好的理解概率密度函数了:
图2

3) 原理应用

直方图均衡化变换就是一种灰度级非线性变换,设r和s分别表示变换前和变换后的灰度,且r和s都进行了归一化的处理。则直方图均衡化变换的公式为:
f4
即归一化后,直方图均衡化的结果s就是r的概率分布函数。

4) 原理推导

根据概率论随机变量的函数的分布的相关知识,有s的概率密度函数为
f5
以下[2]具体论述了其应用过程:
图3
图4
继续推导,有:
f6
其中s为r的概率分布函数,则:
f7
变换后变量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)与直方图规定化

2017-06-30 11:57:00 weixin_34291004 阅读数 782
  • 爬虫简介

    快速入门掌握Python爬虫技术,包括: 1.Python爬虫原理,入门 2.基本用法 3.爬虫应用

    2905人学习 汤小洋
    免费试看

本文主要介绍了灰度直方图相关的处理,包括以下几个方面的内容:

  • 利用OpenCV计算图像的灰度直方图,并绘制直方图曲线
  • 直方图均衡化的原理及实现
  • 直方图规定化(匹配)的原理及实现

图像的灰度直方图

一幅图像由不同灰度值的像素组成,图像中灰度的分布情况是该图像的一个重要特征。图像的灰度直方图就描述了图像中灰度分布情况,能够很直观的展示出图像中各个灰度级所占的多少。
图像的灰度直方图是灰度级的函数,描述的是图像中具有该灰度级的像素的个数:其中,横坐标是灰度级,纵坐标是该灰度级出现的频率。
439761-20170630113705633-1850306826.png

不过通常会将纵坐标归一化到\([0,1]\)区间内,也就是将灰度级出现的频率(像素个数)除以图像中像素的总数。灰度直方图的计算公式如下:
\[ p(r_k) = \frac{n_k}{MN} \]
其中,\(r_k\)是像素的灰度级,\(n_k\)是具有灰度\(r_k\)的像素的个数,\(MN\)是图像中总的像素个数。

OpenCV灰度直方图的计算

直方图的计算是很简单的,无非是遍历图像的像素,统计每个灰度级的个数。在OpenCV中封装了直方图的计算函数calcHist,为了更为通用该函数的参数有些复杂,其声明如下:

void calcHist( const Mat* images, int nimages,
                          const int* channels, InputArray mask,
                          OutputArray hist, int dims, const int* histSize,
                          const float** ranges, bool uniform = true, bool accumulate = false );

该函数能够同时计算多个图像,多个通道,不同灰度范围的灰度直方图.
其参数如下:

  • images,输入图像的数组,这些图像要有相同大大小,相同的深度(CV_8U CV_16U CV_32F).
  • nimages ,输入图像的个数
  • channels,要计算直方图的通道个数。
  • mask,可选的掩码,不使用时可设为空。要和输入图像具有相同的大小,在进行直方图计算的时候,只会统计该掩码不为0的对应像素
  • hist,输出的直方图
  • dims,直方图的维度
  • histSize,直方图每个维度的大小
  • ranges,直方图每个维度要统计的灰度级的范围
  • uniform,是否进行归一化,默认为true
  • accumulate,累积标志,默认值为false。

为了计算的灵活性和通用性,OpenCV的灰度直方图提供了较多的参数,但对于只是简单的计算一幅灰度图的直方图的话,又显得较为累赘。这里对calcHist进行一次封装,能够方便的得到一幅灰度图直方图。

class Histogram1D
{
private:
    int histSize[1]; // 项的数量
    float hranges[2]; // 统计像素的最大值和最小值
    const float* ranges[1];
    int channels[1]; // 仅计算一个通道

public:
    Histogram1D()
    {
        // 准备1D直方图的参数
        histSize[0] = 256;
        hranges[0] = 0.0f;
        hranges[1] = 255.0f;
        ranges[0] = hranges;
        channels[0] = 0;
    }

    MatND getHistogram(const Mat &image)
    {
        MatND hist;
        // 计算直方图
        calcHist(&image ,// 要计算图像的
            1,                // 只计算一幅图像的直方图
            channels,        // 通道数量
            Mat(),            // 不使用掩码
            hist,            // 存放直方图
            1,                // 1D直方图
            histSize,        // 统计的灰度的个数
            ranges);        // 灰度值的范围
        return hist;
    }

    Mat getHistogramImage(const Mat &image)
    {
        MatND hist = getHistogram(image);

        // 最大值,最小值
        double maxVal = 0.0f;
        double minVal = 0.0f;

        minMaxLoc(hist, &minVal, &maxVal);

        //显示直方图的图像
        Mat histImg(histSize[0], histSize[0], CV_8U, Scalar(255));

        // 设置最高点为nbins的90%
        int hpt = static_cast<int>(0.9 * histSize[0]);
        //每个条目绘制一条垂直线
        for (int h = 0; h < histSize[0]; h++)
        {
            float binVal = hist.at<float>(h);
            int intensity = static_cast<int>(binVal * hpt / maxVal);
            // 两点之间绘制一条直线
            line(histImg, Point(h, histSize[0]), Point(h, histSize[0] - intensity), Scalar::all(0));
        }
        return histImg;
    }
};

Histogram1D提供了两个方法:getHistogram返回统计直方图的数组,默认计算的灰度范围是[0,255];getHistogramImage将图像的直方图以线条的形式画出来,并返回包含直方图的图像。测试代码如下:

    Histogram1D hist;
    Mat histImg;
    histImg = hist.getHistogramImage(image);

    imshow("Image", image);
    imshow("Histogram", histImg);

其结果如下:
439761-20170630113736899-1774052938.png

直方图均衡化 Histogram Equalization

假如图像的灰度分布不均匀,其灰度分布集中在较窄的范围内,使图像的细节不够清晰,对比度较低。通常采用直方图均衡化直方图规定化两种变换,使图像的灰度范围拉开或使灰度均匀分布,从而增大反差,使图像细节清晰,以达到增强的目的。
直方图均衡化,对图像进行非线性拉伸,重新分配图像的灰度值,使一定范围内图像的灰度值大致相等。这样,原来直方图中间的峰值部分对比度得到增强,而两侧的谷底部分对比度降低,输出图像的直方图是一个较为平坦的直方图。

均衡化算法

直方图的均衡化实际也是一种灰度的变换过程,将当前的灰度分布通过一个变换函数,变换为范围更宽、灰度分布更均匀的图像。也就是将原图像的直方图修改为在整个灰度区间内大致均匀分布,因此扩大了图像的动态范围,增强图像的对比度。通常均衡化选择的变换函数是灰度的累积概率,直方图均衡化算法的步骤:

  • 计算原图像的灰度直方图 \(P(S_k) = \frac{n_k}{n}\),其中\(n\)为像素总数,\(n_k\)为灰度级\(S_k\)的像素个数
  • 计算原始图像的累积直方图 \(CDF(S_k) = \sum\limits^k_{i=0}\frac{n_i}{n}=\sum\limits^k_{i=0}P_s(S_i)\)
  • \(D_j = L\cdot CDF(S_i)\),其中 \(D_j\)是目的图像的像素,\(CDF(S_i)\)是源图像灰度为i的累积分布,L是图像中最大灰度级(灰度图为255)

其代码实现如下:

  • 在上面中封装了求灰度直方图的类,这里直接应用该方法得到图像的灰度直方图;
  • 将灰度直方图进行归一化,计算灰度的累积概率;
  • 创建灰度变化的查找表
  • 应用查找表,将原图像变换为灰度均衡的图像

具体代码如下:

void equalization_self(const Mat &src, Mat &dst)
{
    Histogram1D hist1D;
    MatND hist = hist1D.getHistogram(src);

    hist /= (src.rows * src.cols); // 对得到的灰度直方图进行归一化
    float cdf[256] = { 0 }; // 灰度的累积概率
    Mat lut(1, 256, CV_8U); // 灰度变换的查找表
    for (int i = 0; i < 256; i++)
    {
        // 计算灰度级的累积概率
        if (i == 0)
            cdf[i] = hist.at<float>(i);
        else
            cdf[i] = cdf[i - 1] + hist.at<float>(i);

        lut.at<uchar>(i) = static_cast<uchar>(255 * cdf[i]); // 创建灰度的查找表
    }

    LUT(src, lut, dst); // 应用查找表,进行灰度变化,得到均衡化后的图像

}

上面代码只是加深下对均衡化算法流程的理解,实际在OpenCV中也提供了灰度均衡化的函数equalizeHist,该函数的使用很简单,只有两个参数:输入图像,输出图像。下图为,上述代码计算得到的均衡化结果和调用equalizeHist的结果对比
439761-20170630113810774-1163114065.png

最左边为原图像,中间为OpenCV封装函数的结果,右边为上面代码得到的结果。

直方图规定化

从上面可以看出,直方图的均衡化自动的确定了变换函数,可以很方便的得到变换后的图像,但是在有些应用中这种自动的增强并不是最好的方法。有时候,需要图像具有某一特定的直方图形状(也就是灰度分布),而不是均匀分布的直方图,这时候可以使用直方图规定化
直方图规定化,也叫做直方图匹配,用于将图像变换为某一特定的灰度分布,也就是其目的的灰度直方图是已知的。这其实和均衡化很类似,均衡化后的灰度直方图也是已知的,是一个均匀分布的直方图;而规定化后的直方图可以随意的指定,也就是在执行规定化操作时,首先要知道变换后的灰度直方图,这样才能确定变换函数。规定化操作能够有目的的增强某个灰度区间,相比于,均衡化操作,规定化多了一个输入,但是其变换后的结果也更灵活。

在理解了上述的均衡化过程后,直方图的规定化也较为简单。可以利用均衡化后的直方图作为一个中间过程,然后求取规定化的变换函数。具体步骤如下:

  • 将原始图像的灰度直方图进行均衡化,得到一个变换函数\(s = T(r)\),其中s是均衡化后的像素,r是原始像素
  • 对规定的直方图进行均衡化,得到一个变换函数\(v = G(z)\),其中v是均衡化后的像素,z是规定化的像素
  • 上面都是对同一图像的均衡化,其结果应该是相等的,\(s = v,且 z = G^{-1}(v) = G^{-1}(T(r))\)

通过,均衡化作为中间结果,将得到原始像素\(r\)\(z\)规定化后像素之间的映射关系。

详解规定化过程

对图像进行直方图规定化操作,原始图像的直方图和以及规定化后的直方图是已知的。假设\(P_r(r)\)表示原始图像的灰度概率密度,\(P_z(z)\)表示规定化图像的灰度概率密度(r和z分别是原始图像的灰度级,规定化后图像的灰度级)。

  • 对原始图像进行均衡化操作,则有\(s_k = T(r_k) = L \cdot \sum\limits_{i=0}^{i=k}P_r(r_k)\)
  • 对规定化的直方图进行均衡化操作,则\(v_k = G(z_m) = L \cdot \sum\limits_{j=0}^{j=m}P_z(z_m)\)
  • 由于是对同一图像的均衡化操作,所以有\(s_k = v_m\)
  • 规定化操作的目的就是找到原始图像的像素\(s_k\)到规定化后图像像素的\(z_k\)之间的一个映射。有了上一步的等式后,可以得到\(s_k = G(z_k)\),因此要想找到\(s_k\)想对应的\(z_k\)只需要在\(z\)进行迭代,找到使式子\(G(z_m)-s_k\)的绝对值最小即可。
  • 上述描述只是理论的推导过程,在实际的计算过程中,不需要做两次的均衡化操作,具体的推导过程如下:\[ \begin{array}{c} s_k = v_k \\ L \cdot \sum\limits_{i=0}^{i=k}P_r(r_k) = L \cdot \sum\limits_{j=0}^{j=m}P_z(z_m) \\ \sum\limits_{i=0}^{i=k}P_r(r_k) = \sum\limits_{j=0}^{j=m}P_z(z_m) \end{array} \]
    上面公式表示,假如\(s_k\) 规定化后的对应灰度是\(z_m\)的话,需要满足的条件是\(s_k\)的累积概率和\(z_m\)的累积概率是最接近的
    下面是一个具体计算的例子:
    439761-20170630114016602-177041903.png

首先得到原直方图的各个灰度级的累积概率\(V_s\)以及规定化后直方图的各个灰度级的累积概率\(V_z\),那么确定\(s_k\)\(z_m\)之间映射关系的条件就是:\[\mid V_s - V_z \mid\]的值最小。
\(k = 2\)为例,其原始直方图的累积概率是:0.65,在规定化后的直方图的累积概率中和0.65最接近(相等)的是灰度值为5的累积概率密度,则可以得到原始图像中的灰度级2,在规定化后的图像中的灰度级是5

直方图规定化的实现

直方图规定化的实现可以分为一下三步:

  • 计算原图像的累积直方图
  • 计算规定直方图的累积直方图
  • 计算两累积直方图的差值的绝对值
  • 根据累积直方图差值建立灰度级的映射

具体代码实现如下:

void hist_specify(const Mat &src, const Mat &dst,Mat &result)
{
    Histogram1D hist1D;
    MatND src_hist = hist1D.getHistogram(src);
    MatND dst_hist = hist1D.getHistogram(dst);

    float src_cdf[256] = { 0 };
    float dst_cdf[256] = { 0 };

    // 源图像和目标图像的大小不一样,要将得到的直方图进行归一化处理
    src_hist /= (src.rows * src.cols);
    dst_hist /= (dst.rows * dst.cols);

    // 计算原始直方图和规定直方图的累积概率
    for (int i = 0; i < 256; i++)
    {
        if (i == 0)
        {
            src_cdf[i] = src_hist.at<float>(i);
            dst_cdf[i] = dst_hist.at<float>(i);
        }
        else
        {
            src_cdf[i] = src_cdf[i - 1] + src_hist.at<float>(i);
            dst_cdf[i] = dst_cdf[i - 1] + dst_hist.at<float>(i);
        }
    }

    // 累积概率的差值
    float diff_cdf[256][256];
    for (int i = 0; i < 256; i++)
        for (int j = 0; j < 256; j++)
            diff_cdf[i][j] = fabs(src_cdf[i] - dst_cdf[j]);

    // 构建灰度级映射表
    Mat lut(1, 256, CV_8U);
    for (int i = 0; i < 256; i++)
    {
        // 查找源灰度级为i的映射灰度
        // 和i的累积概率差值最小的规定化灰度
        float min = diff_cdf[i][0];
        int index = 0;
        for (int j = 1; j < 256; j++)
        {
            if (min > diff_cdf[i][j])
            {
                min = diff_cdf[i][j];
                index = j;
            }
        }
        lut.at<uchar>(i) = static_cast<uchar>(index);
    }

    // 应用查找表,做直方图规定化
    LUT(src, lut, result);
}

上面函数的第二个参数的直方图就是规定化的直方图。代码比较简单,这里就不一一解释了。其结果如下:
439761-20170630114404664-1309652532.png

左边是原图像,右边是规定化的图像,也就是上面函数的第一个和第二个输入参数。原图像规定化的结果如下:
439761-20170630115536164-248882789.png

原图像规定化后的直方图和规定化的图像的直方图的形状比较类似, 并且原图像规定化后整幅图像的特征和规定化的图像也比较类似,例如:原图像床上的被子,明显带有规定化图像中水的波纹特征。

直方图规定化过程中,在做灰度映射的时候,有两种常用的方法:

  • 单映射 Single Mapping Law,SML,这种方法也是上面使用的方法,根据累积直方图的差值,从原图像中找到其在规定化图像中的映射。
  • 组映射 Group Mapping Law,GML 这种方法较上述方法复杂不少,但是处理效果较好。

对于GML的映射方法,一直没有很好的理解,但是根据其算法描述实现了该方法,代码这里先不放出,其处理结果如下:
439761-20170630115831571-1515651251.png

其结果较SML来说更为亮一些,床上的波浪特征也更为明显,但是其直方图形状,和规定化的直方图对比,第一个峰不是很明显。

总结

  • 图像的灰度直方图能够很直观的展示图像中灰度级的整体分布情况,对图像的后续处理有很好的指导作用。
  • 直方图的均衡化的是将一幅图像的直方图变平,使各个灰度级的趋于均匀分布,这样能够很好的增强图像对比度。直方图均衡化是一种自动化的变换,仅需要输入图像,就能够确定图像的变换函数。但是直方图的均衡化操作也有一定的确定,在均衡化的过程中对图像中的数据不加选择,这样有可能会增强图像的背景;变换后图像的灰度级减少,有可能造成某些细节的消失;会压缩图像直方图中的高峰,造成处理后图像对比度的不自然等。
  • 直方图规定化,也称为直方图匹配,经过规定化处理将原图像的直方图变换为特定形状的直方图(上面中的示例,就是将图像的直方图变换为另一幅图像的直方图)。它可以按照预先设定的某个形状来调整图像的直方图,运用均衡化原理的基础上,通过建立原始图像和期望图像之间的关系,选择地控制直方图,使原始图像的直方图变成规定的形状它可以按照预先设定的某个形状来调整图像的直方图。直方图规定化是在运用均衡化原理的基础上,通过建立原始图像和期望图像之间的关系,选择地控制直方图,使原始图像的直方图变成规定的形状,从而弥补直方图均衡化的一些缺点.
2020-04-07 12:53:39 qq_41650371 阅读数 386
  • 爬虫简介

    快速入门掌握Python爬虫技术,包括: 1.Python爬虫原理,入门 2.基本用法 3.爬虫应用

    2905人学习 汤小洋
    免费试看

实验题目:

编程实现灰度和彩色图像的直方图均衡化处理。要求给出原始图像的直方图、均衡化图像及其直方图和直方图均衡化时所用的灰度级变换曲线图。注意彩色图像需要编程实现图像由RGB色彩空间到HSI的变换,然后在亮度通道上进行直方图均衡化。

源码:

对灰度图像的处理,给出原始图像的直方图、均衡化图像及其直方图和直方图均衡化时所用的灰度级变换曲线图

// 灰度图像
clear 
%读入彩色图像将其灰度化
pic=imread('grey.png');	
imshow(pic)		%显示出来
title('输入的灰色png图像')

%绘制直方图
[m,n]=size(pic);	%测量图像尺寸参数
GP=zeros(1,256);	%预创建存放灰度出现概率的向量
N=zeros(1,256);
for k=0:255
    N(k+1) = length(find(pic==k));    
    GP(k+1)=length(find(pic==k))/(m*n); 	% 计算每级灰度出现的概率,将其存入 GP 中相应位置
end
figure,bar(0:255,GP,'b')	%绘制直方图title('原图像直方图')
xlabel('灰度值') 
ylabel('出现概率')

%直方图均衡化
CF=zeros(1,256);
for i=1:256
    for j=1:i
        CF(i)=GP(j)+CF(i);	%计算 Sk
    end
end
S2=round((CF*256)+0.5);	%将 Sk 归到相近级的灰度
for i=1:256
    GPeq(i)=sum(GP(find(S2==i)));	%计算现有每个灰度级出现的概率
end
figure,bar(0:255,GPeq,'b')	%显示均衡化后的直方图title('均衡化后的直方图')
xlabel('灰度值') 
ylabel('出现概率')

%图像均衡化
PA=pic;
for i=0:255
PA(find(pic==i))=S2(i+1);	%将各个像素归一化后的灰度值赋给这个像素
end
figure,imshow(PA)	

%绘制灰度级变换曲线
figure
plot(0:255,S2,'r')
xlabel('均值化前')
ylabel('均值化后')
grid on
legend('灰度级变换曲线');

对彩色图像的处理,给出原始图像的直方图、均衡化图像及其直方图和直方图均衡化时所用的灰度级变换曲线图。注意彩色图像需要编程实现图像由RGB色彩空间到HSI的变换,然后在亮度通道上进行直方图均衡化。

// 彩色图像
function pic2()
clear;
pic=imread('color.jpg');	%读入 JPG 彩色图像文件
imshow(pic)		%显示出来
title('原图像');

[M,N,D] = size(pic);
%提取单通道分量
pic=double(pic);
r=pic(:,:,1);
g=pic(:,:,2);
b=pic(:,:,3);
%实现转换
angle=acos(0.5*((r-g)+(r-b))./(sqrt((r-g).^2+(r-b).*(g-b))));
if (b>=g)
    H = 2*pi-angle;
else
    H = angle;
end
H=H/(2*pi);
S=1-3.*(min(min(r,g),b))./(r+g+b);
H(S==0)=0;
I=(r+g+b)/3;
I=uint8(I);%范围不超[0,255]
%统计像素亮度
INumPixel = zeros(1,256);%用长度为256的一维数组统计各灰度值的数目
for i = 1:M  
    for j = 1: N  
        INumPixel(I(i,j) + 1) = INumPixel(I(i,j) + 1) + 1;  
    end  
end  
%亮度直方图
GP = zeros(1,256);  
for i = 1:256  
    GP(i) = INumPixel(i) / (M * N * 1.0);  
end  
figure,bar(0:255,GP,'b')	%绘制直方图title('原图像直方图')
xlabel('亮度值') 
ylabel('出现概率')

%直方图均衡化
S2 = zeros(1,256);  
for i = 1:256  
    if i == 1  
        S2(i) = GP(i);  
    else  
        S2(i) = S2(i - 1) + GP(i);  
    end  
end  
                                %S2乘上最大灰度255并且向上取整  
S2 = uint8(255 .* S2 + 0.5);  
                                %将原图像各个位置灰度值映射到新值  
for i = 1:M  
    for j = 1: N  
        Inew(i,j) = S2(I(i,j)+1);  
    end 
end
for i=1:256
    GPeq(i)=sum(GP(find(S2==i)));	%计算现有每个灰度级出现的概率
end
figure
plot(0:255,GPeq,'b')	%显示均衡化后的直方图title('均衡化后的直方图')
xlabel('亮度值') 
ylabel('出现概率')

%绘制亮度级变换曲线
figure
plot(0:1:255,S2,'r')
xlabel('均值化前')
ylabel('均值化后')
grid on
legend('亮度级变换曲线');

I=Inew;
H=H*2*pi;
I=double(I);
S=double(S);
H=double(H);
[m,n]=size(H);
%转换
for i = 1:m
    for j = 1:n
        if (0<=H(i,j))&(H(i,j)<2*pi/3)
            B(i,j)=I(i,j).*(1-S(i,j));
            R(i,j)=I(i,j).*(1+S(i,j).*cos(H(i,j))./cos(pi/3-H(i,j)));
            G(i,j)=3*I(i,j)-(R(i,j)+B(i,j));
        end
        if (2*pi/3<=H(i,j))&(H(i,j)<4*pi/3)
            H(i,j)=H(i,j)-2*pi/3;
            R(i,j)=I(i,j).*(1-S(i,j));
            G(i,j)=I(i,j).*(1+S(i,j).*cos(H(i,j)-2*pi/3)./cos(pi-H(i,j)));
            B(i,j)=3*I(i,j)-(R(i,j)+G(i,j));
        end
        if (4*pi/3<=H(i,j))& (H(i,j)<2*pi)
            H(i,j)=H(i,j)-4*pi/3;
            G(i,j)=I(i,j).*(1-S(i,j));
            B(i,j)=I(i,j).*(1+S(i,j).*cos(H(i,j)-4*pi/3)./cos(5*pi/3-H(i,j)));
            R(i,j)=3*I(i,j)-(G(i,j)+B(i,j));
        end
    end
end

output3=cat(3,R,G,B);
output3=uint8(output3);
figure
imshow(output3)

实验结果

灰度图像
灰度原图 在这里插入图片描述
在这里插入图片描述 在这里插入图片描述
在这里插入图片描述

彩色图像
在这里插入图片描述晋亨话之秋

在这里插入图片描述 在这里插入图片描述
在这里插入图片描述

通过原始图像和均衡化后的图像对比,原始图像直方图与均衡化后直方图对比,均衡化后的直方图灰度值更加平均,在整幅图像中不再集中。