本人QQ:1154911397,欢饮大家与我交流！！！\(^o^)/~
今天要给大家分享的是：
在彩色图像上进行Mean Shift迭代搜索目标
R\G\B三维加权直方图+巴氏系数+Mean Shift迭代
关于 加权直方图、巴氏系数、Mean Shift迭代 这三者之间的关系请大侠们阅读我的另一篇博文：加权直方图+巴氏系数+Mean Shift的关系
关于一维直方图在灰度图像上的迭代程序请看：一维直方图+巴氏系数+Mean Shift
关于二维直方图在彩色图像上的迭代程序请看：二维直方图+巴氏系数+Mean Shift
请看程序：
// KernalBasedMeanShift_3.cpp : 定义控制台应用程序的入口点。
//
#include "stdafx.h"
#include
#include
#include
using namespace std;
using namespace cv;
Mat image; //母图像
Mat roi_hist; //兴趣区域图像块的直方图
Mat roi_whist ; //兴趣区域图像块的加权直方图
Mat roi_img; //兴趣区域的图像块
Rect roi_rect(153,165,116,146); //ROI图像块在原图的真正区域
Mat ms_hist; //由MeanShift给出的目标区域的直方图
bool selectObject = false; //selectObject的初始值为false,
int calculateHist = 0; //calculateHist的初值是0
bool showHist = false;
Point origin;
Rect selection;
/*选择将要被统计的通道编号*/
vector
select_channels(3) ;
//设定直方图中bin的数目
int histSize[3] = /*{8,8,16}*/ /*{8,8,8}*/ {16,16,16}/*{10,10,10}*/;
//设定取值范围
float range[] = {0,256};//灰度特征的范围都是[0,256)
const float* histRange[] = {range,range,range};
bool uniform = true; //均匀分布，
bool accumu = false;//无累加
//定义一个核函数指针，该函数指针可以指向任意一个标准的核函数
typedef float (*KernalFunc)(Point2f& center,Point2f& xy,int hx,int hy);
int weight_method = 0; //选择核函数 ==0：Epanechnikov；==1：Gaussian
KernalFunc DeriveFuncOfKernal ;//声明一个指向核函数的导函数的函数指针变量
//用于控制是否激活鼠标选取roi区域的功能
bool lock_roi = true; //若果为真，则roi区域被锁定为默认的区域，如果为假，则使用鼠标选中的矩形框
//用鼠标选择目标图像区域
static void onMouse( int event, int x, int y, int, void* );
//该函数用于计算给定矩形图像块的加权直方图
static void CalculateWeightedHistogram(const Mat& img,vector
& channels, Mat& _hist,int weighted_method); //该函数用于计算给定矩形图像块的非加权直方图 static void CalculateHistogram(const Mat& img ,vector
& channels, Mat& _hist); //自己写的计算图像三维直方图的函数 void myCalcHist3D( const Mat& image, vector
& channels,OutputArray _hist, int dims, const int* histSize,const float** ranges); //计算图像二维加权直方图的函数 static void myCalcWeightedHist3D( const Mat& image, vector
& 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); /*核函数的负导数：Derivation of Epanechnikov Kernal*/ static float DerivationOfEpanechnikovKernal(Point2f& center,Point2f& xy,int hx,int hy); /*核函数：Gaussian Kernal */ static float GaussianKernalFunc(Point2f& center,Point2f& xy,int hx,int hy); /*核函数的负导数：Derivation of Gaussian Kernal*/ static float DerivationOfGaussianKernal(Point2f& center,Point2f& xy,int hx,int hy); //计算两个直方图的巴式距离 static float ComputeBhattacharyyaDistance3D(const Mat& hist1 ,const Mat& hist2 ); //根据当前搜索框下的图像块，加权直方图以及目标模板直方图计算更新权值矩阵 static void UpdateWeightsMatrix3D(Mat& image,vector
& channels,Mat& modal_hist,Mat& cur_hist, int* histSize,const float** ranges,OutputArray weights_matrix); //根据均值向量漂移到新的位置(公式26) static void ShiftNewLocationByMeanVector3D(Mat& weights_matrix,KernalFunc DeriveFuncOfKernal, Point& leftup,Point2f& old_center,Point2f& new_center); //自己编写的Mean Shift算法，isJudgeOverShift参数用于控制是否在迭代过程中判断漂移冲过了头； //最后两个参数用来保存迭代路径，用于显示搜索窗口的运动和记录窗口的巴氏系数 int MyMeanShift3D(Mat& image, Mat& modal_hist, Rect& initialWindow,TermCriteria criteria , bool isJudgeOverShift,vector
& iterRects ,vector
& iterBhattaCoeff); int _tmain(int argc, _TCHAR* argv[]) { string filename = "meinv1.jpg"; string winname = "MainWindow"; namedWindow(winname,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(dst); src.copyTo(image); if(src.channels() >2) { select_channels[0]=0; //blue通道 select_channels[1] =1;//green通道 select_channels[2] =2;//red通道 } else { printf("输入图像的通道数必须大于1！\n");getchar(); } /**********控制算法运行行为的主要参数*******************************************************/ //该参数用于控制给定的初始位置框和真正的roi框的重叠度：取值范围(0,1) float overlap = 0.3f; //一般来说，重叠度越高，迭代收敛越快，迭代次数越少 int MaxIterNum = 50; //该参数用于控制Mean Shift的最大迭代次数 bool isJudgeOverShift = true;//用于Mean Shift迭代过程中判断是否冲过头的标志 weight_method = 0; //选择加权核函数 ==0的话选择 Epanechnikov kernal；==1选择Gaussian kernal lock_roi = true; //如果lock_roi为真，则鼠标选取的矩形框被无效化，使用roi_rect的默认位置作为roi区域 /**********控制算法运行行为的主要参数*******************************************************/ for(;;) { ///循环显示读入的图像 imshow(winname,src); if(calculateHist) { roi_img = image(selection);//提取鼠标选中的图像块 imshow("ROI Image",roi_img);//显示roi 图像块 calculateHist = 0; //计算选中的图像块的直方图 CalculateHistogram(roi_img,select_channels,roi_hist); 计算选中图像块的加权直方图 CalculateWeightedHistogram(roi_img,select_channels,roi_whist,weight_method); //调用MeanShift算法在迭代寻找目标图像块 Rect InitialWindow(roi_rect); //初始搜索框位置默认放在真实位置的右下角，与真实目标框有overlap的重叠度 InitialWindow.x += cvRound((1-overlap)*roi_rect.width); InitialWindow.y += cvRound((1-overlap)*roi_rect.height); cout<<"roi区域的精确位置： "<
<
0 && selection.height > 0 ) calculateHist = 1; //当在第一帧用鼠标选定了合适的目标跟踪窗口后，calculateHist的值置为 1 //如果lock_roi为真，则鼠标选取被无效化，使用roi_rect的默认位置作为roi区域 if(lock_roi) selection = roi_rect; cout<<"选中的矩形区域为： "<
<
& channels, Mat& _hist,int weighted_method) { //调用自己写的函数计算选中图像块的加权直方图 KernalFunc kernal;//声明一个KernalFunc类型的函数指针变量 //根据加权参数指定具体的核函数 if(weighted_method == 0) { kernal = EpanechnikovKernalFunc; //调用Epanechnikov核函数 DeriveFuncOfKernal = DerivationOfEpanechnikovKernal;//Epanechnikov的导函数 } else if(weighted_method == 1) { kernal = GaussianKernalFunc; //调用高斯Gaussian核函数 DeriveFuncOfKernal = DerivationOfGaussianKernal;//高斯核函数的导函数 } myCalcWeightedHist3D(img,channels,_hist,3,histSize,histRange,kernal); /// 将直方图归一化到范围 [ 0, 1] //normalize(_hist, _hist, 0, 1, NORM_MINMAX, -1, Mat() );//该函数不支持三维直方图的归一化 } //该函数用于计算给定矩形图像块的直方图 static void CalculateHistogram(const Mat& img,vector
& channels, Mat& _hist) { //调用OpenCV函数计算选中图像块的灰度直方图 //calcHist(&img,1,0,Mat(),_hist,1,&histSize,&histRange,uniform,accumu); //调用自己写的函数计算选中图像块的灰度直方图 myCalcHist3D(img,channels,_hist,3,histSize,histRange); /// 将直方图归一化到范围 [ 0, 1] //normalize(_hist, _hist, 0, 1, CV_MINMAX, -1, Mat() );//该函数不支持三维直方图的归一化 return; } /** 计算给定图像的三维直方图，图像类型必须是CV_8UCx的,x = 3, or 4 or ...., channels规定了通道索引顺序，对于3D直方图,channels数组只能有3个值，且必须都 < x */ void myCalcHist3D( const Mat& image, vector
& 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是输出参数，保存着计算好的直方图 _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 ch1 = channels[dims-3], ch2 = channels[dims-2],ch3 = channels[dims-1]; //对应于直方图的第一维的区间范围，bin个数等 float low_range1 = ranges[0][0] ,high_range1 = ranges[0][1];//要计算的直方图第一维区间的上下限 float a1 = histSize[0]/(high_range1 - low_range1);/* 第一维单位区间内直方图bin的数目*/ float b1 = -a1*low_range1; //对应于直方图的第二维的区间范围，bin个数等 float low_range2= ranges[1][0] ,high_range2 = ranges[1][1];//要计算的直方图第二维区间的上下限 float a2 = histSize[1]/(high_range2 - low_range2);/* 单位区间内直方图bin的数目*/ float b2 = -a2*low_range2; //对应于直方图的第三维的区间范围，bin个数等 float low_range3 = ranges[2][0] ,high_range3 = ranges[2][1];//要计算的直方图第三维区间的上下限 float a3 = histSize[2]/(high_range3 - low_range3);/* 单位区间内直方图bin的数目*/ float b3 = -a3*low_range3; const uchar* hist_data = hist.data; //指向直方图矩阵的数据区的首指针 size_t hist_step0 = hist.step[0]; //直方图第一维的步长 size_t hist_step1 = hist.step[1]; //直方图第er维的步长 size_t hist_step2 = hist.step[2]; //直方图第san维的步长 Mat_
Myx = image;//Mxy是一个uchar类型的三元组，用于取出image中位于(x,y)位置的像素 //外循环按行从上往下扫描 for(int y=0; y < img_height ; y++ ) { //内循环从左往右扫描 for(int x=0; x < img_width ; x++) { //取出输入图像的第y行第x列所有通道的像素值 Vec3b val = Myx(y,x);//val三元组中按照B,G,R的顺序存放着三个通道的uchar值 //计算像素值val在输入参数规定的灰度区间内的bin序号的索引值 int idx1 = cvFloor(val[ch1]*a1 + b1);//输出直方图第一维的bin索引号 int idx2 = cvFloor(val[ch2]*a2 + b2);//输出直方图第二维的bin索引号 int idx3 = cvFloor(val[ch3]*a3 + b3);//输出直方图第三维的bin索引号 //第一维idx1对应于直方图矩阵的面(层)索引,第二维idx2对应于行索引，第三维idx3为列索引 float* hist_ptr =(float*)(hist_data + idx1*hist_step0 + idx2*hist_step1+ idx3*hist_step2); (*hist_ptr)++; //累计(idx1,idx2,idx3)位置处的对应的直方图bin上的数量 } } return; } //计算图像二维加权直方图的函数 static void myCalcWeightedHist3D( const Mat& image, vector
& channels,OutputArray _hist, int dims, const int* histSize,const float** ranges,KernalFunc kernal) { 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]行histSize[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 ch1 = channels[dims-3], ch2 = channels[dims-2],ch3 = channels[dims-1]; //对应于直方图的第一维的区间范围，bin个数等 float low_range1 = ranges[0][0] ,high_range1 = ranges[0][1];//要计算的直方图第一维区间的上下限 float a1 = histSize[0]/(high_range1 - low_range1);/* 第一维单位区间内直方图bin的数目*/ float b1 = -a1*low_range1; //对应于直方图的第二维的区间范围，bin个数等 float low_range2= ranges[1][0] ,high_range2 = ranges[1][1];//要计算的直方图第二维区间的上下限 float a2 = histSize[1]/(high_range2 - low_range2);/* 单位区间内直方图bin的数目*/ float b2 = -a2*low_range2; //对应于直方图的第三维的区间范围，bin个数等 float low_range3 = ranges[2][0] ,high_range3 = ranges[2][1];//要计算的直方图第三维区间的上下限 float a3 = histSize[2]/(high_range3 - low_range3);/* 单位区间内直方图bin的数目*/ float b3 = -a3*low_range3; const uchar* hist_data = hist.data; //指向直方图矩阵的数据区的首指针 size_t hist_step0 = hist.step[0]; //直方图第一维的步长，层(或叫“面”)索引 size_t hist_step1 = hist.step[1]; //直方图第er维的步长，某层图像的行索引 size_t hist_step2 = hist.step[2]; //直方图第san维的步长，某层图像的列索引 Mat_
Myx = image;//Mxy是一个uchar类型的三元组，用于取出image中位于(x,y)位置的像素 //外循环按行从上往下扫描 for(int y=0; y < img_height ; y++ ) { //指向图像第y行的指针 const uchar* ptr_y = image.ptr
(y); //内循环从左往右扫描 for(int x=0; x < img_width ; x++) { //取出输入图像的第y行第x列所有通道的像素值 Vec3b val = Myx(y,x);//val三元组中按照B,G,R的顺序存放着三个通道的uchar值 //计算像素值val在输入参数规定的灰度区间内的bin序号的索引值 int idx1 = cvFloor(val[ch1]*a1 + b1);//输出直方图第一维的bin索引号 int idx2 = cvFloor(val[ch2]*a2 + b2);//输出直方图第二维的bin索引号 int idx3 = cvFloor(val[ch3]*a3 + b3);//输出直方图第三维的bin索引号 //第一维idx1对应于直方图矩阵的面(层)索引,第二维idx2对应于行索引，第三维idx3为列索引 float* hist_ptr =(float*)(hist_data + idx1*hist_step0 + idx2*hist_step1+ idx3*hist_step2); // 调用核函数，计算图像块（x，y）位置处的加权值 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 idx1=0; idx1 < histSize[0]; idx1++) { for(int idx2=0; idx2 < histSize[1]; idx2++) { for(int idx3=0; idx3 < histSize[2]; idx3++) { //第一维idx1对应于直方图矩阵的面(层)索引,第二维idx2对应于行索引，第三维idx3为列索引 float* hist_ptr =(float*)(hist_data + idx1*hist_step0 + idx2*hist_step1+ idx3*hist_step2); *hist_ptr *= NomalConst_inv; //乘以归一化因子 } } } 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; //hx是沿着x轴横向半径 float disty = (xy.y - center.y)/hy; //hy是沿着y轴纵向半径 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; } /*核函数的负导数：Derivation of Epanechnikov Kernal*/ static float DerivationOfEpanechnikovKernal(Point2f& center,Point2f& xy,int hx,int hy) { //计算点xy到中心center的归一化距离 float distx = (xy.x - center.x)/hx; //hx是沿着x轴横向半径 float disty = (xy.y - center.y)/hy; //hy是沿着y轴纵向半径 float dist_square = distx*distx + disty*disty;//距离平方 float result = 0.f; //函数要返回的结果 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; //坐标的维数 //result = 0.5f*Cd_inv*( dimension + 2 ) ; result = 0.63661977237f; //是个常数 = 0.5f*Cd_inv*( dimension + 2 ) } 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; //hx是沿着x轴横向半径 float disty = (xy.y - center.y)/hy; //hy是沿着y轴纵向半径 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; } /*核函数的负导数：Derivation of Gaussian Kernal*/ static float DerivationOfGaussianKernal(Point2f& center,Point2f& xy,int hx,int hy) { //计算点xy到中心center的归一化距离 float distx = (xy.x - center.x)/hx; //hx是沿着x轴横向半径 float disty = (xy.y - center.y)/hy; //hy是沿着y轴纵向半径 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 = 12.566370614f; // ==4*CV_PI float Cd_inv = 1.0f/Cd; int dimension = 2; //坐标的维数 result = Cd_inv*expf( -0.5f*dist_square ); //高斯核函数的导函数 } return result; } //计算两个二维直方图的巴式距离 static float ComputeBhattacharyyaDistance3D(const Mat& hist1 ,const Mat& hist2 ) { //判断两个直方图的维数是否一致 int hist_layer = hist1.size[0] ; //三维直方图矩阵的层数 == histSize[0] int hist_rows = hist1.size[1] ; //一层图像的行数 == histSize[1] int hist_cols = hist1.size[2] ; //一层图像的列数 == histSize[2] if((hist_layer != hist2.size[0])||(hist_rows != hist2.size[1]) ||(hist_cols != hist2.size[2])) { printf("两个直方图的维数不一致,无法计算距离！！！\n");getchar(); } const uchar* hist1_data = hist1.data; //指向直方图矩阵的数据区的首指针 const uchar* hist2_data = hist2.data; //指向直方图矩阵的数据区的首指针 size_t hist_step0 = hist1.step[0]; //直方图第一维的步长 size_t hist_step1 = hist1.step[1]; //直方图第er维的步长 size_t hist_step2 = hist1.step[2]; //直方图第san维的步长 float distance = 0.0f; const float* phist1 =0,*phist2 =0; for(int idx1=0; idx1 < hist_layer; idx1++) { for(int idx2=0; idx2 < hist_rows; idx2++) { for(int idx3=0; idx3 < hist_cols; idx3++) { //第一维idx1对应于直方图矩阵的面(层)索引,第二维idx2对应于行索引，第三维idx3为列索引 phist1 =(float*)(hist1_data + idx1*hist_step0 + idx2*hist_step1+ idx3*hist_step2); phist2 =(float*)(hist2_data + idx1*hist_step0 + idx2*hist_step1+ idx3*hist_step2); distance += sqrtf((*phist1)*(*phist2));//累计距离值 } } } return distance; } //根据当前搜索框下的图像块，加权直方图以及目标模板直方图计算更新权值矩阵 static void UpdateWeightsMatrix3D(Mat& cur_image,vector
& channels,Mat& modal_hist,Mat& cur_hist, int* histSize,const float** ranges,OutputArray weights_matrix) { if((cur_image.channels() < 2 )||(cur_image.depth() != CV_8U)) { printf("该函数只支持CV_8UC3类型的多通道图像\n"); getchar(); } int img_channel = cur_image.channels(); int rows = cur_image.rows; int cols = cur_image.cols; //给权值矩阵开辟空间，权值矩阵的大小与image的大小一样，每一个权值与image的一个像素位置一一对应 weights_matrix.create(rows,cols,CV_32FC1);//创建一个rows行cols列的单通道浮点数矩阵 Mat wei_mat = weights_matrix.getMat(); //从OutputArray类型的参数获取Mat类型的矩阵头索引 //遍历图像块的每一个像素位置，找到该像素点在直方图中的bin索引号，然后求两个直方图对应bin位置处的比值 int dims = 3; //直方图维数 //被统计的通道索引，对应直方图的第一维和第二维,和第三维 int ch1 = channels[dims-3], ch2 = channels[dims-2],ch3 = channels[dims-1]; //对应于直方图的第一维的区间范围，bin个数等 float low_range1 = ranges[0][0] ,high_range1 = ranges[0][1];//要计算的直方图第一维区间的上下限 float a1 = histSize[0]/(high_range1 - low_range1);/* 第一维单位区间内直方图bin的数目*/ float b1 = -a1*low_range1; //对应于直方图的第二维的区间范围，bin个数等 float low_range2= ranges[1][0] ,high_range2 = ranges[1][1];//要计算的直方图第二维区间的上下限 float a2 = histSize[1]/(high_range2 - low_range2);/* 单位区间内直方图bin的数目*/ float b2 = -a2*low_range2; //对应于直方图的第三维的区间范围，bin个数等 float low_range3 = ranges[2][0] ,high_range3 = ranges[2][1];//要计算的直方图第三维区间的上下限 float a3 = histSize[2]/(high_range3 - low_range3);/* 单位区间内直方图bin的数目*/ float b3 = -a3*low_range3; size_t hist_step0 = modal_hist.step[0]; //直方图第一维的步长，层(或叫“面”)索引 size_t hist_step1 = modal_hist.step[1]; //直方图第er维的步长，某层图像的行索引 size_t hist_step2 = modal_hist.step[2]; //直方图第san维的步长，某层图像的列索引 //按照从上往下，从左往右的顺序填充权值矩阵 for(int row =0; row < rows; row++) { //指向image的第row行的指针 uchar* curimg_ptr = cur_image.ptr
(row); //指向权重矩阵第row行的指针 float* weimat_ptr = wei_mat.ptr
(row); for(int col =0; col < cols; col++) { //取出输入图像的第y行第x列第ch1、ch2、ch3 通道的像素值 uchar val1 = curimg_ptr[col*img_channel + ch1]; uchar val2 = curimg_ptr[col*img_channel + ch2]; uchar val3 = curimg_ptr[col*img_channel + ch3]; //计算灰度值val在输入参数规定的灰度区间内的bin序号的索引值 int idx1 = cvFloor(val1*a1 + b1);//直方图第一维的bin索引号 int idx2 = cvFloor(val2*a2 + b2);//直方图第二维的bin索引号 int idx3 = cvFloor(val3*a3 + b3);//直方图第二维的bin索引号 //取出(idx1,idx2)对应位置的直方图的bin值 float* pmodal = (float*)(modal_hist.data + idx1*hist_step0 + idx2*hist_step1+ idx3*hist_step2); float* pcur = (float*)(cur_hist.data + idx1*hist_step0 + idx2*hist_step1+ idx3*hist_step2); //第row行第col列的权重值 float weight = sqrtf( (*pmodal) / (*pcur + DBL_EPSILON) ); *(weimat_ptr + col) = weight; } } return; } //根据均值向量漂移到新的位置(公式26) static void ShiftNewLocationByMeanVector3D(Mat& weights_matrix,KernalFunc DeriveFuncOfKernal, Point& leftup,Point2f& old_center,Point2f& new_center) { //权值矩阵的行列数 int rows = weights_matrix.rows; int cols = weights_matrix.cols; //核函数的横向半径和纵向半径 int hx = cols/2 , hy = rows/2; //计算新的搜索框中心位置 float x_sum=0.f,y_sum=0.f,const_sum=0.f; int x = 0 , y = 0; for(int row =0; row < rows; row++) { //指向权重矩阵第row行的指针 float* weimat_ptr = weights_matrix.ptr
(row); y = leftup.y + row; for(int col =0; col < cols; col++) { x = leftup.x + col; //当前像素点在母图像中的坐标位置 Point2f point(x,y); //调用对应核函数的导函数 float g = DeriveFuncOfKernal(old_center,point,hx,hy); //取出对应位置（row，col）上的权重值 float weight = weimat_ptr[col]; //累加分子与分母 x_sum += x*weight*g ; y_sum += y*weight*g ; const_sum += weight*g ; } } //计算新的中心点 new_center.x = cvRound(x_sum/(const_sum+DBL_EPSILON)); new_center.y = cvRound(y_sum/(const_sum+DBL_EPSILON)); } //最后两个参数用来保存迭代路径，用于显示搜索窗口的运动和记录窗口的巴氏系数 int MyMeanShift3D(Mat& image, Mat& modal_hist, Rect& initialWindow,TermCriteria criteria , bool isJudgeOverShift,vector
& iterRects ,vector
& iterBhattaCoeff) { int i = 0, eps; Rect cur_rect ; //当前搜索窗口 Mat cur_win; //当前搜索窗口下的图像块 Mat cur_hist; //当前搜索窗口下图像块的加权直方图 Mat cur_weights;//当前搜索窗口下的权值矩阵 //初始化迭代路径, iterRects.clear(); iterBhattaCoeff.clear(); Mat mat = image; Rect windowIn(initialWindow); if( mat.channels()<2 )//目前只支持多通道图像 CV_Error( CV_BadNumChannels, "目前只支持多通道图像" ); if( windowIn.height <= 0 || windowIn.width <= 0 ) CV_Error( CV_StsBadArg, "窗口尺寸必须大于0" ); windowIn = Rect(windowIn) & Rect(0, 0, mat.cols, mat.rows); //边界检查 //检查迭代终止条件 criteria = cvCheckTermCriteria( criteria, 1., 100 ); eps = cvRound( criteria.epsilon * criteria.epsilon ); //初始搜索框 cur_rect = windowIn; float BhattaCoeff = 0.0f; //开始MeanShift迭代过程，实际迭代次数小于等于criteria.maxCount for( i = 0; i < criteria.maxCount; i++ ) { int dx, dy, nx, ny; //当前搜索窗口的中心 Point2f cur_center(cur_rect.x +0.5f* cur_rect.width, cur_rect.y + 0.5f* cur_rect.height); //当前搜索窗口的左上角顶点坐标 Point2i cur_leftup(cur_rect.x,cur_rect.y); //从当前帧图像中获取当前搜索窗口下的图像块 cur_win = mat(cur_rect); //计算当前窗口下的图像块的加权直方图 CalculateWeightedHistogram(cur_win,select_channels,cur_hist,weight_method); //计算当前搜索窗口下的直方图与给定的目标模板直方图的巴氏系数 BhattaCoeff = ComputeBhattacharyyaDistance3D(modal_hist,cur_hist); //根据公式25计算当前窗口下的权值矩阵， UpdateWeightsMatrix3D(cur_win,select_channels,modal_hist,cur_hist, histSize,histRange,cur_weights); //公式-25 //根据公式26计算均值漂移向量，得到新的搜索窗口的中心位置 Point2f new_center;//新的搜索窗口的中心位置 ShiftNewLocationByMeanVector3D(cur_weights,DeriveFuncOfKernal, cur_leftup,cur_center,new_center); //公式-26 //新的搜索区域的左上角坐标 nx = cvRound(new_center.x - 0.5f*cur_rect.width); ny = cvRound(new_center.y - 0.5f*cur_rect.height); //x方向的边界检查 if( nx < 0 ) nx = 0; else if( nx + cur_rect.width > mat.cols ) nx = mat.cols - cur_rect.width; //y方向的边界检查 if( ny < 0 ) ny = 0; else if( ny + cur_rect.height > mat.rows ) ny = mat.rows - cur_rect.height; //计算漂移向量(dx,dy) dx = nx - cur_rect.x; dy = ny - cur_rect.y; //记录迭代过程 iterRects.push_back(cur_rect); iterBhattaCoeff.push_back(BhattaCoeff); //形成新的搜索框，移动左上角顶点坐标，长宽不变 Rect new_rect(cur_rect); new_rect.x = nx; new_rect.y = ny; //判断均值漂移是否冲过了头 if(isJudgeOverShift) { Mat new_win = mat(new_rect); //取出新窗口下的图像块 Mat new_hist; //计算新窗口下的图像块的加权直方图 CalculateWeightedHistogram(new_win,select_channels,new_hist,0); //计算新窗口下的直方图与给定的目标模板直方图的巴氏系数 float new_BhattaCoeff = ComputeBhattacharyyaDistance3D(modal_hist,new_hist); //如果新位置上的巴氏系数小于旧位置上的的巴氏系数，说明冲过了头 if(new_BhattaCoeff < BhattaCoeff) { //若果冲过了头，就把两次的位置坐标平均一下，得出新的位置 new_rect.x = cvRound(0.5f*(cur_rect.x + new_rect.x)); new_rect.y = cvRound(0.5f*(cur_rect.y + new_rect.y)); //由于调整了新窗口的位置，所以要重新计算漂移向量(dx,dy) dx = new_rect.x - cur_rect.x; dy = new_rect.y - cur_rect.y; } } /* 检查窗口的移动距离是否小于eps,若小于，则算法收敛，退出循环 */ if( dx*dx + dy*dy < eps ) { cur_rect = new_rect; //如果在没有达到最大迭代次数之前该退出条件已经满足 //则iterRects中的最后一个元素就是最终的迭代结果 iterRects.push_back(cur_rect); iterBhattaCoeff.push_back(BhattaCoeff); break; } //进入下一次迭代 cur_rect = new_rect; //这个搜索框在循环退出时才被加入到迭代轨迹数组中 } //如果超过最大迭代次数的时候，循环退出，把最后一次均值漂移形成的矩形框作为最终的结果 if( i == criteria.maxCount) { iterRects.push_back(cur_rect); iterBhattaCoeff.push_back(BhattaCoeff); } return iterRects.size();//返回实际使用的迭代次数(漂移次数) }
请看运行结果：
扫一扫
8661
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
