核密度函数加权直方图的计算

在本篇博客中，我将与众位看官分享我自己写的关于用核密度函数加权的直方图的计算

欢饮批评指正！！！

核密度函数加权直方图在基于Mean shift的跟踪算法中经常用到。请客官查看我的博客：

样本均值漂移的原理 Mean Shift 原理

常用的核密度函数有

EpanechnikovKernal：

 

 

 

在Mean Shift的迭代过程中用的是上述两个函数的截面函数：

下面是这两个函数(4)和(5)的程序实现：

 

1.  

2. /*核函数：Epanechnikov Kernal

3. center: 核函数中心坐标

4. xy : 在核区域内的某一个点

5. hx,hy : 核区域在x方向和y方向的半径（或叫 带宽）

6. */

7. static float EpanechnikovKernalFunc(Point2f& center,Point2f& xy,int hx,int hy)

8. {

9. //计算点xy到中心center的归一化距离

10. float distx = (xy.x - center.x)/hx;

11. float disty = (xy.y - center.y)/hy;

12. float dist_square = distx*distx + disty*disty;//距离平方

13.  

14. float result = 0.f; //函数要返回的结果

15.  

16. //核区域就是以hx和hy为边长的矩形的内接圆，在该内接圆中的点

17. //的函数值都不为0，若距离超过核区域，则返回0，

18. //距离center越近，函数值越大

19. if(dist_square>=1)

20. {

21. result = 0.f;

22. }

23. else

24. {

25. //float Cd = CV_PI; //单位圆面积 pi*R^2 ,R==1

26. float Cd_inv = 0.318309886f; // == 1.0/Cd;单位圆的面积的倒数

27. int dimension = 2; //坐标的维数

28. //Epanechnikov核函数的截面断层函数：profile

29. result = 0.5f*Cd_inv*( dimension + 2 )*( 1 - dist_square );

30. }

31.  

32. return result;

33. }

34.  

35. /*核函数：Gaussian Kernal

36. center: 核函数中心坐标

37. xy : 在核区域内的某一个点

38. hx,hy : 核区域在x方向和y方向的半径（或叫 带宽）

39. */

40. static float GaussianKernalFunc(Point2f& center,Point2f& xy,int hx,int hy)

41. {

42. //计算点xy到中心center的归一化距离

43. float distx = (xy.x - center.x)/hx;

44. float disty = (xy.y - center.y)/hy;

45. float dist_square = distx*distx + disty*disty;//距离平方

46.  

47. float result = 0.f; //函数要返回的结果

48.  

49. //核区域就是以hx和hy为边长的矩形的内接圆，在该内接圆中的点

50. //的函数值都不为0，若距离超过核区域，则返回0，

51. //距离center越近，函数值越大

52. if(dist_square>=1)

53. {

54. result = 0.f;

55. }

56. else

57. {

58. float Cd = 6.2831853072f; // ==2*CV_PI

59. float Cd_inv = 1.0f/Cd;

60. int dimension = 2; //坐标的维数

61. result = Cd_inv*expf( -0.5f*dist_square ); //高斯核函数的截面断层函数：profile

62. }

63.  

64. return result;

65. }

66.  

 

下面是调用上述核函数实现的一维加权直方图统计函数 

该函数中的最后一个参数是函数指针类型的参数，用于指定到底用哪个核函数

KernalFunc的定义如下：

//定义一个核函数指针，该函数指针可以指向任意一个符合参数标准的核函数
typedef float (*KernalFunc)(Point2f& center,Point2f& xy,int hx,int hy);

 

1. //计算图像一维加权直方图的函数

2. static void myCalcWeightedHist1D( const Mat& image, vector

3.  

4. & channels,OutputArray _hist,

5. int dims, const int* histSize,const float** ranges,KernalFunc kernal)

6. {

7.  

8. //KernalFunc kernal;//声明一个KernalFunc类型的函数指针变量

9.  

10. int calc_channels = channels.size();//image图像中指定要被统计计算的通道的个数

11. int img_channel = image.channels(); //image图像的总的通道个数

12. if(img_channel < channels[calc_channels - 1] + 1 )

13. {

14. printf("channels中的最大通道索引数不能超过images的通道数\n");

15. getchar();

16. }

17. if(image.depth() != CV_8U)

18. {

19. printf("该函数仅支持CV_8U类的图像\n");

20. getchar();

21. }

22. if(dims != calc_channels)

23. {

24. printf("被索引的通道数目与直方图的维数不相等\n");

25. getchar();

26. }

27. int img_width = image.cols; //图像宽度

28. int img_height = image.rows;//图像高度

29.  

30. //图像区域的中心，即核函数中心

31. Point2f center(0.5f*img_width,0.5f*img_height);

32. int hx = img_width/2 , hy = img_height/2;//核函数带宽

33. float NormalConst =0.f; //用于归一化加权系数

34. //参数_hist是输出参数，保存着计算好的直方图，一维直方图是histSize[0]行1列的矩阵

35. _hist.create(dims, histSize, CV_32F); //为输出直方图开辟空间

36. Mat hist = _hist.getMat();

37. hist.flags = (hist.flags & ~CV_MAT_TYPE_MASK)|CV_32S;//直方图矩阵的类型信息

38. hist = Scalar(0.f); //清空原来的数据，当前直方图不累加到原来的直方图上去

39.  

40. int ch = channels[dims-1];//要统计的通道编号 ==: 0 or 1 or 2

41. float low_range = ranges[0][0] ,high_range = ranges[0][1];//要计算的直方图区间的上下限

42. float a = histSize[0]/(high_range - low_range);/* 单位区间内直方图bin的数目*/

43. float b = -a*low_range;

44.  

45. //外循环按行从上往下扫描

46. for(int y=0; y < img_height ; y++ )

47. {

48. //指向图像第y行的指针

49. const uchar* ptr_y = image.ptr

50.  

51. (y);

52.  

53. //内循环从左往右扫描

54. for(int x=0; x < img_width ; x++)

55. {

56. //取出输入图像的第y行第x列第ch通道的像素值

57. uchar val = ptr_y[x*img_channel + ch];

58.  

59. //灰度值val在输入参数规定的灰度区间内的bin序号的索引值

60. int idx = cvFloor(val*a + b);//该索引值是数学上线性运算后得到的

61. float* hist_ptr = hist.ptr

62.  

63. (idx); //指向输出直方图的第idx个bin

64. float weighted_value = kernal(center,Point2f(x,y),hx,hy); //（x,y）位置处的加权值

65. (*hist_ptr) += weighted_value; //加权累计直方图的每一个bin

66. NormalConst += weighted_value; //加权系数归一化因子

67. }

68. }

69.  

70. float NomalConst_inv = 1.0f/(NormalConst + DBL_EPSILON);

71. //归一化加权直方图

72. for(int ii=0; ii < histSize[0]; ii++)

73. {

74. float* hist_ptr = hist.ptr

75.  

76. (ii); //指向输出直方图的第ii个bin

77. *hist_ptr *= NomalConst_inv;

78. }

79. return;

80. }

81.  
82.  
83.  
84.  

 

 

下面是测试上述各个小函数的主函数：

说明：在主函数中，计算了三种不同的一维直方图。

一种是非加权直方图，它的计算函数可以参见我的博客：自己编写的直方图函数

第二种是Epanikov 核函数加权的直方图

第三种是Gaussian核函数加权的直方图

主函数： 

 

1. //基于核函数的加权直方图计算

2.  

3. #include "stdafx.h"

4. #include

5.  
6.  

7. #include

8.  
9.  

10. #include

11.  
12.  
13.  

14. using namespace std;

15. using namespace cv;

16.  

17. Mat image; //母图像

18. Mat roi_hist; //兴趣区域图像块的直方图

19. Mat roi_whist1,roi_whist2;//兴趣区域图像块的加权直方图

20. Mat roi_img; //兴趣区域的图像块

21. Rect roi_rect(162,159,100,122); //ROI图像块在原图的真正区域

22. Mat ms_hist; //由MeanShift给出的目标区域的直方图

23. Mat ProbImg; //反射概率分布图像

24. bool selectObject = false; //selectObject的初始值为false,

25. int calculateHist = 0; //calculateHist的初值是0

26. bool showHist = false;

27. Point origin;

28. Rect selection;

29.  

30. //指定要统计的通道编号

31. vector

32.  

33. selected_channel(1,0);

34. //设定直方图中bin的数目

35. int histSize = 128;

36. //设定取值范围

37. float range[] = {0,256};//灰度特征的范围都是[0,256)

38. const float* histRange = {range};

39. bool uniform = true; //均匀分布，

40. bool accumu = false;//无累加

41. int hist_w = 400; int hist_h = 400;//直方图图像的宽度和高度

42. int bin_w = cvRound( (double) hist_w/histSize );//直方图中一个矩形条纹的宽度

43. // 创建直方图画布

44. Mat histImage( hist_w, hist_h, CV_8UC3, Scalar( 0,0,0) );//创建画布图像

45. //用鼠标选择目标图像区域

46. static void onMouse( int event, int x, int y, int, void* );

47. //该函数用于计算给定矩形图像块的加权直方图

48. static void CalculateWeightedHistogram(const Mat& img,vector

49.  

50. & channels, Mat& _hist,int weighted_method); //该函数用于计算给定矩形图像块的非加权直方图 static void CalculateHistogram(const Mat& img,vector

51.  

52. & channels, Mat& _hist); //自己写的计算图像一维直方图的函数 void myCalcHist1D( const Mat& image, vector

53.  

54. & channels,OutputArray _hist, int dims, const int* histSize,const float** ranges); //定义一个核函数指针，该函数指针可以指向任意一个标准的核函数 typedef float (*KernalFunc)(Point2f& center,Point2f& xy,int hx,int hy); //计算图像一维加权直方图的函数 static void myCalcWeightedHist1D( const Mat& image, vector

55.  

56. & channels,OutputArray _hist, int dims, const int* histSize,const float** ranges,KernalFunc kernal); /*核函数：Epanechnikov Kernal*/ static float EpanechnikovKernalFunc(Point2f& center,Point2f& xy,int hx,int hy); /*核函数：Gaussian Kernal */ static float GaussianKernalFunc(Point2f& center,Point2f& xy,int hx,int hy); //绘制直方图 void DrawHistogramImage( Mat& _hist , Scalar& color); int _tmain(int argc, _TCHAR* argv[]) { string filename = "meinv1.jpg"; string winname = "MainWindow"; namedWindow(winname,1); namedWindow("Hist Image",1); namedWindow("ROI Image",1); //为 Histogram Demo窗口添加鼠标响应函数，用于选择图像块 setMouseCallback( "MainWindow", onMouse, 0 ); Mat src, dst;//声明原始图像和目标图像 /// 装载图像 src = imread( filename, 1 ); if( !src.data ) //判断图像是否成功读取 { return -1; } //src.copyTo(image); cvtColor(src,image,CV_BGR2GRAY); for(;;) { ///循环显示读入的图像 imshow(winname,src); if(calculateHist) { roi_img = image(selection);//提取鼠标选中的图像块 imshow("ROI Image",roi_img);//显示roi 图像块 calculateHist = 0; //计算选中的图像块的直方图 CalculateHistogram(roi_img,selected_channel,roi_hist); //计算选中图像块的加权直方图 int weighted_method = 0; //指定加权核函数的种类 CalculateWeightedHistogram(roi_img,selected_channel,roi_whist1,weighted_method); weighted_method = 1; //使用高斯核函数 CalculateWeightedHistogram(roi_img,selected_channel,roi_whist2,weighted_method); //清空画布 histImage = Scalar::all(0); //绘制直方图 DrawHistogramImage(roi_hist,Scalar(0,0,255)); //绘制加权直方图 DrawHistogramImage(roi_whist1,Scalar(0,255,0)); //绘制加权直方图 DrawHistogramImage(roi_whist2,Scalar(255,0,0)); cout<<"红色曲线： 未加权直方图"<

57.  

58. 0 && selection.height > 0 ) calculateHist = 1; //当在第一帧用鼠标选定了合适的目标跟踪窗口后，calculateHist的值置为 1 //selection = roi_rect; //放开这一句注释，则鼠标选取被无效化，使用roi_rect的默认位置作为roi区域 cout<<"选中的矩形区域为： "<

59.  

60. <

61.  

62. & channels, Mat& _hist,int weighted_method) { //调用自己写的函数计算选中图像块的加权直方图 KernalFunc kernal;//声明一个KernalFunc类型的函数指针变量 //根据加权参数指定具体的核函数 if(weighted_method == 0) kernal = EpanechnikovKernalFunc; //调用Epanechnikov核函数 else if(weighted_method == 1) kernal = GaussianKernalFunc; //调用高斯Gaussian核函数 myCalcWeightedHist1D(img,channels,_hist,1,&histSize,&histRange,kernal); /// 将直方图归一化到范围 [ 0, 1] normalize(_hist, _hist, 0, 1, NORM_MINMAX, -1, Mat() ); } //该函数用于计算给定矩形图像块的非加权直方图 static void CalculateHistogram(const Mat& img,vector

63.  

64. & channels, Mat& _hist) { //调用OpenCV函数计算选中图像块的灰度直方图 //calcHist(&img,1,0,Mat(),_hist,1,&histSize,&histRange,uniform,accumu); //调用自己写的函数计算选中图像块的灰度直方图 myCalcHist1D(img,channels,_hist,1,&histSize,&histRange); /// 将直方图归一化到范围 [ 0, 1] normalize(_hist, _hist, 0, 1, NORM_MINMAX, -1, Mat() ); } //绘制直方图 void DrawHistogramImage( Mat& _hist , Scalar& color) { //用于绘图的直方图矩阵数组 Mat draw_hist; /// 将直方图归一化到范围 [ 0, histImage.rows ] normalize( _hist , draw_hist , 0, histImage.rows, NORM_MINMAX, -1, Mat() ); /// 在直方图画布上画出直方图 /*histImage = Scalar::all(0);*/ for(int i=1;i<=histSize;i++) { //矩形图表示 /* rectangle( histImage,Point((i-1)*bin_w,hist_h),Point(i*bin_w,hist_h-cvRound(draw_hist.at

65.  

66. (i-1))),Scalar(0,0,255),1,8,0); */ //折线图表示 line( histImage, Point( bin_w*(i-1), hist_h - cvRound(draw_hist.at

67.  

68. (i-1)) ) , Point( bin_w*(i), hist_h - cvRound(draw_hist.at

69.  

70. (i)) ), color, 1, 8, 0 ); } /// 显示直方图 imshow("Hist Image", histImage ); } /** 计算给定图像的一维直方图，图像类型必须是CV_8UCx的,x = 1,or 2, or 3 or ...., channels规定了通道索引顺序，对于1D直方图,channels数组只能有一个值，且必须 < x */ void myCalcHist1D( const Mat& image, vector

71.  

72. & channels,OutputArray _hist, int dims, const int* histSize,const float** ranges) { int calc_channels = channels.size();//image图像中指定要被统计计算的通道的个数 int img_channel = image.channels(); //image图像的总的通道个数 if(img_channel < channels[calc_channels - 1] + 1 ) { printf("channels中的最大通道索引数不能超过images的通道数\n"); getchar(); } if(image.depth() != CV_8U) { printf("该函数仅支持CV_8U类的图像\n"); getchar(); } if(dims != calc_channels) { printf("被索引的通道数目与直方图的维数不相等\n"); getchar(); } int img_width = image.cols; //图像宽度 int img_height = image.rows;//图像高度 //参数_hist是输出参数，保存着计算好的直方图，一维直方图是histSize[0]行1列的矩阵 _hist.create(dims, histSize, CV_32F); //为输出直方图开辟空间 Mat hist = _hist.getMat(); hist.flags = (hist.flags & ~CV_MAT_TYPE_MASK)|CV_32S;//直方图矩阵的类型信息 hist = Scalar(0.f); //清空原来的数据，当前直方图不累加到原来的直方图上去 int ch = channels[dims-1];//要统计的通道编号 ==: 0 or 1 or 2 float low_range = ranges[0][0] ,high_range = ranges[0][1];//要计算的直方图区间的上下限 float a = histSize[0]/(high_range - low_range);/* 单位区间内直方图bin的数目*/ float b = -a*low_range; //外循环按行从上往下扫描 for(int y=0; y < img_height ; y++ ) { //指向图像第y行的指针 const uchar* ptr_y = image.ptr

73.  

74. (y); //内循环从左往右扫描 for(int x=0; x < img_width ; x++) { //取出输入图像的第y行第x列第ch通道的像素值 uchar val = ptr_y[x*img_channel + ch]; //灰度值val在输入参数规定的灰度区间内的bin序号的索引值 int idx = cvFloor(val*a + b);//该索引值是数学上线性运算后得到的 float* hist_ptr = hist.ptr

75.  

76. (idx); //指向输出直方图的第idx个bin (*hist_ptr)++; //累计idx对应的直方图bin上的数量 } } return; } /*核函数：Epanechnikov Kernal center: 核函数中心坐标 xy : 在核区域内的某一个点 hx,hy : 核区域在x方向和y方向的半径（或叫 带宽） */ static float EpanechnikovKernalFunc(Point2f& center,Point2f& xy,int hx,int hy) { //计算点xy到中心center的归一化距离 float distx = (xy.x - center.x)/hx; float disty = (xy.y - center.y)/hy; float dist_square = distx*distx + disty*disty;//距离平方 float result = 0.f; //函数要返回的结果 //核区域就是以hx和hy为边长的矩形的内接圆，在该内接圆中的点 //的函数值都不为0，若距离超过核区域，则返回0， //距离center越近，函数值越大 if(dist_square>=1) { result = 0.f; } else { //float Cd = CV_PI; //单位圆面积 pi*R^2 ,R==1 float Cd_inv = 0.318309886f; // == 1.0/Cd;单位圆的面积的倒数 int dimension = 2; //坐标的维数 //Epanechnikov核函数的截面断层函数：profile result = 0.5f*Cd_inv*( dimension + 2 )*( 1 - dist_square ); } return result; } /*核函数：Gaussian Kernal center: 核函数中心坐标 xy : 在核区域内的某一个点 hx,hy : 核区域在x方向和y方向的半径（或叫 带宽） */ static float GaussianKernalFunc(Point2f& center,Point2f& xy,int hx,int hy) { //计算点xy到中心center的归一化距离 float distx = (xy.x - center.x)/hx; float disty = (xy.y - center.y)/hy; float dist_square = distx*distx + disty*disty;//距离平方 float result = 0.f; //函数要返回的结果 //核区域就是以hx和hy为边长的矩形的内接圆，在该内接圆中的点 //的函数值都不为0，若距离超过核区域，则返回0， //距离center越近，函数值越大 if(dist_square>=1) { result = 0.f; } else { float Cd = 6.2831853072f; // ==2*CV_PI float Cd_inv = 1.0f/Cd; int dimension = 2; //坐标的维数 result = Cd_inv*expf( -0.5f*dist_square ); //高斯核函数的截面断层函数：profile } return result; } //计算图像一维加权直方图的函数 static void myCalcWeightedHist1D( const Mat& image, vector

77.  

78. & channels,OutputArray _hist, int dims, const int* histSize,const float** ranges,KernalFunc kernal) { //KernalFunc kernal;//声明一个KernalFunc类型的函数指针变量 int calc_channels = channels.size();//image图像中指定要被统计计算的通道的个数 int img_channel = image.channels(); //image图像的总的通道个数 if(img_channel < channels[calc_channels - 1] + 1 ) { printf("channels中的最大通道索引数不能超过images的通道数\n"); getchar(); } if(image.depth() != CV_8U) { printf("该函数仅支持CV_8U类的图像\n"); getchar(); } if(dims != calc_channels) { printf("被索引的通道数目与直方图的维数不相等\n"); getchar(); } int img_width = image.cols; //图像宽度 int img_height = image.rows;//图像高度 //图像区域的中心，即核函数中心 Point2f center(0.5f*img_width,0.5f*img_height); int hx = img_width/2 , hy = img_height/2;//核函数带宽 float NormalConst =0.f; //用于归一化加权系数 //参数_hist是输出参数，保存着计算好的直方图，一维直方图是histSize[0]行1列的矩阵 _hist.create(dims, histSize, CV_32F); //为输出直方图开辟空间 Mat hist = _hist.getMat(); hist.flags = (hist.flags & ~CV_MAT_TYPE_MASK)|CV_32S;//直方图矩阵的类型信息 hist = Scalar(0.f); //清空原来的数据，当前直方图不累加到原来的直方图上去 int ch = channels[dims-1];//要统计的通道编号 ==: 0 or 1 or 2 float low_range = ranges[0][0] ,high_range = ranges[0][1];//要计算的直方图区间的上下限 float a = histSize[0]/(high_range - low_range);/* 单位区间内直方图bin的数目*/ float b = -a*low_range; //外循环按行从上往下扫描 for(int y=0; y < img_height ; y++ ) { //指向图像第y行的指针 const uchar* ptr_y = image.ptr

79.  

80. (y); //内循环从左往右扫描 for(int x=0; x < img_width ; x++) { //取出输入图像的第y行第x列第ch通道的像素值 uchar val = ptr_y[x*img_channel + ch]; //灰度值val在输入参数规定的灰度区间内的bin序号的索引值 int idx = cvFloor(val*a + b);//该索引值是数学上线性运算后得到的 float* hist_ptr = hist.ptr

81.  

82. (idx); //指向输出直方图的第idx个bin float weighted_value = kernal(center,Point2f(x,y),hx,hy); //（x,y）位置处的加权值 (*hist_ptr) += weighted_value; //加权累计直方图的每一个bin NormalConst += weighted_value; //加权系数归一化因子 } } float NomalConst_inv = 1.0f/(NormalConst + DBL_EPSILON); //归一化加权直方图 for(int ii=0; ii < histSize[0]; ii++) { float* hist_ptr = hist.ptr

83.  

84. (ii); //指向输出直方图的第ii个bin *hist_ptr *= NomalConst_inv; } return; }

85.  
86.  
87.  
88.  
89.  
90.  
91.  
92.  
93.  
94.  
95.  
96.  
97.  
98.  
99.  
100.  
101.  
102.  
103.  
104.  
105.  

 

 

实验结果：

上图显示了对美女脸部图像的加权直方图的统计结果(实际计算的时候是在灰度图像上计算的），

红色曲线是未加权直方图，它的中间部分有一个小高突起，是因为美女的脸部灰度值在100到150区间内，

蓝色曲线是高斯加权直方图的结果，可以看到中间的小峰比红色曲线高，这是由于用高斯函数进行了加权，把脸部中间的直方图bin的灰度频率值提高了

绿色曲线是Epanechikov核函数的加权结果，可以看到中间的小峰比蓝色曲线还要高，这是因为Epanechikov核函数比高斯核函数在中心区域更加陡峭所致，所以加权效果也比高斯函数要大的多