数字图像处理实验（九）：meanshift跟踪算法

一、实验原理

MeanShift算法

无参密度估计理论，无参密度估计也叫做非参数估计，属于数理统计的一个分支，和参数密度估计共同构成了概率密度估计方法。
参数密度估计方法：要求特征空间服从一个已知的概率密度函数，在实际的应用中这个条件很难达到。
无参数密度估计方法：对先验知识要求最少，完全依靠训练数据进行估计，并且可以用于任意形状的密度估计。所以依靠无参密度估计方法，即不事先规定概率密度函数的结构形式，在某一连续点处的密度函数值可由该点邻域中的若干样本点估计得出。常用的无参密度估计方法有：直方图法、最近邻域法和核密度估计法。
MeanShift算法正是属于核密度估计法，它不需要任何先验知识而完全依靠特征空间中样本点的计算其密度函数值。对于一组采样数据，直方图法通常把数据的值域分成若干相等的区间，数据按区间分成若干组，每组数据的个数与总参数个数的比率就是每个单元的概率值；核密度估计法的原理相似于直方图法，只是多了一个用于平滑数据的核函数。采用核函数估计法，在采样充分的情况下，能够渐进地收敛于任意的密度函数，即可以对服从任何分布的数据进行密度估计。

（1）均值漂移的基本形式

给定d维空间的n个数据点集X，那么对于空间中的任意点x的mean shift向量基本形式可以表示为：

这个向量就是漂移向量，其中S_k表示的是数据集的点到x的距离小于球半径h的数据点，即：

而漂移的过程，说的简单一点，就是通过计算得漂移向量，然后把球圆心x的位置更新一下，更新公式为：

使得圆心的位置一直处于力的平衡位置：

简而言之：就是求解一个向量，使得圆心一直往数据集密度最大的方向移动。
说的再简单一点，就是每次迭代的时候，都是找到圆里面点的平均位置作为新的圆心位置。

（2）加入核函数的漂移向量

这个说的简单一点就是加入一个高斯权重，最后的漂移向量计算公式为。这里对比一下（1）中的计算公式：

加入一个高斯权重：

每次更新的圆心坐标：

（3）mean shift 聚类流程

假设在一个多维空间中有很多数据点需要进行聚类，Mean Shift的过程如下：
①、在未被标记的数据点中随机选择一个点作为中心center；
②、找出离center距离在bandwidth之内的所有点，记做集合M，认为这些点属于簇c。同时，把这些求内点属于这个类的概率加1，这个参数将用于最后步骤的分类
③、以center为中心点，计算从center开始到集合M中每个元素的向量，将这些向量相加，得到向量shift。
④、center = center+shift。即center沿着shift的方向移动，移动距离是||shift||。
⑤、重复步骤2、3、4，直到shift的大小很小（就是迭代到收敛），记住此时的center。注意，这个迭代过程中遇到的点都应该归类到簇c。
⑥、如果收敛时当前簇c的center与其它已经存在的簇c2中心的距离小于阈值，那么把c2和c合并。否则，把c作为新的聚类，增加1类。
⑦、重复1、2、3、4、5直到所有的点都被标记访问。
⑧、分类：根据每个类，对每个点的访问频率，取访问频率最大的那个类，作为当前点集的所属类。
简单的说，mean shift就是沿着密度上升的方向寻找同属一个簇的数据点。

二、实验代码

实验环境：
（1）OpenCV3.4.3
（2）Ubuntu16.04
（3）VS Code
（4）C++


//
#include <iostream>
#include "opencv2/opencv.hpp"
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/core/core.hpp>
#include <math.h>using namespace cv;
using namespace std;bool selectObject = false;
int trackObject = 0;
cv::Rect selection;
// 被跟踪的图片
cv::Mat tracking_image;
cv::Point origin;// 3 鼠标选择函数
static void onMouse(int event, int x, int y, int, void*) {if (selectObject) {selection.x = MIN(x, origin.x);//选择区域的x坐标选起点与当前点的最小值，保证鼠标不管向右下角还是左上角拉动都正确选择selection.y = MIN(y, origin.y);selection.width = std::abs(x - origin.x);selection.height = std::abs(y - origin.y);selection &= cv::Rect(0, 0, tracking_image.cols, tracking_image.rows);//确保所选矩形在图片范围内}switch (event) {case cv::EVENT_LBUTTONDOWN:  //按下鼠标左键，进行赋初值origin = cv::Point(x, y);selection = cv::Rect(x, y, 0, 0);selectObject = true;break;case cv::EVENT_LBUTTONUP:   //鼠标左键抬起selectObject = false;if (selection.width > 0 && selection.height > 0)trackObject = -1;   // 置-1，主程序开始进行meanshiftbreak;}
}// 实验七：canny算子边缘检测
class Exp7{
public:void initialize(int pic_id){makeGaussTemplate(5);computerGaussTemplateSum();makeSobelTemplate();gradX = cv::Mat::zeros(gray_image_list[pic_id].rows, gray_image_list[pic_id].cols, CV_64F); // 水平梯度gradY = cv::Mat::zeros(gray_image_list[pic_id].rows, gray_image_list[pic_id].cols, CV_64F); // 垂直梯度grad = cv::Mat::zeros(gray_image_list[pic_id].rows, gray_image_list[pic_id].cols, CV_8U);  // 梯度幅值thead = cv::Mat::zeros(gray_image_list[pic_id].rows, gray_image_list[pic_id].cols, CV_64F); // 梯度角度locateGrad = cv::Mat::zeros(gray_image_list[pic_id].rows, gray_image_list[pic_id].cols, CV_64F);  //区域grad_max = 0;}// 0.1、彩色图像转灰度图像cv::Mat color2Gray(cv::Mat src_image) {//创建与原图同类型和同大小的矩阵cv::Mat gray_image(src_image.rows, src_image.cols, CV_8UC1);if (src_image.channels() != 1) {for (int i = 0; i < src_image.rows; i++)for (int j = 0; j < src_image.cols; j++)gray_image.at<uchar>(i, j) = (src_image.at<cv::Vec3b>(i, j)[0] + src_image.at<cv::Vec3b>(i, j)[1] + src_image.at<cv::Vec3b>(i, j)[2]) / 3;}elsegray_image = src_image.clone();return gray_image;}// 0.2 零填充cv::Mat grayZerosPadding(cv::Mat& src, int kernel_size=5){int m = kernel_size/2;cv::Mat dst = cv::Mat::zeros(2*m + src.rows, 2*m + src.cols, CV_8U);dst(cv::Rect(m,m,src.cols,src.rows)) += src;return dst;}// 1 高斯滤波// 1.1 高斯模板void makeGaussTemplate(int size=5, int sigma=1){cv::Mat gaussTemplate = cv::Mat::zeros(size, size, CV_32F);int center=size/2;double min = g(center,center);for(int i=0; i < size; i++)for(int j=0; j < size; j++)gaussTemplate.at<float>(i, j) = g(i-center,j-center)/min;gaussTemplate.convertTo(gaussTemplate, CV_8U);gauss_template.push_back(gaussTemplate);}// 1.2 计算正态分布double g(double x, double y, double sigma=1){return exp(-(x*x + y*y)/(2*sigma*sigma));}// 1.3 高斯模板和void computerGaussTemplateSum(){gauss_template_sum = 0;for(int i=0; i < gauss_template.rows; i++)for(int j=0; j < gauss_template.cols; j++)gauss_template_sum += gauss_template.at<uchar>(i, j);}// 1.4 生成sobel模板void makeSobelTemplate(){sobel_template.push_back((cv::Mat_<float>(3,3) <<  -1,   0,   1, -2,   0,   2, -1,   0,   1));sobel_template.push_back((cv::Mat_<float>(3,3) <<  -1,  -2,  -1, 0,   0,   0, 1,   2,   1));}// 2 canny边缘检测算法void Canny(int pic_id, double lowThres, double hightThres){cv::Mat result = cv::Mat::zeros(gray_image_list[pic_id].rows, gray_image_list[pic_id].cols, CV_8UC1);// 1 高斯滤波grayGaussFiltering(pic_id);// 2 sobel算子,计算梯度强度 与 方向sobelDetection(pic_id);/******************* 3 非极大值抑制 && 4 双阈值检测 *******************/// 4.1 确定双阈值// 将灰度幅值拉成一维向量，便于排序，求解灰度幅值的最大值std::vector<int> caculateValue;convert2vector(grad, caculateValue);qsort_vector(caculateValue, 0, caculateValue.size() - 1);long long highIndex = grad.rows * grad.cols * hightThres;// 确定双阈值:int highValue = caculateValue[highIndex]; int lowValue = highValue * lowThres;for(int i = 1 ; i < grad.rows-1; i++ )for( int j = 1; j < grad.cols-1; j++){// 八个方位double N = grad.at<uchar>(i-1, j);     // 正上（北）double NE = grad.at<uchar>(i-1, j+1);  // 右上方（东北）double E = grad.at<uchar>(i, j+1);     // 右边（东）double SE = grad.at<uchar>(i+1, j+1);  // 右下方（东南）double S = grad.at<uchar>(i+1, j);     // 正下（南）double SW = grad.at<uchar>(i-1, j-1);  // 左下方（西南）double W = grad.at<uchar>(i, j-1);     // 左边（西）double NW = grad.at<uchar>(i -1, j -1);// 左上方（西北）// 区域判断，线性插值处理double tanThead;    // tan角度double Gp1;         // 两个方向的梯度强度double Gp2;// 求角度，绝对值tanThead = abs(tan(thead.at<float>(i,j)));// 求解正向梯度 与 反向梯度switch ((int)locateGrad.at<float>(i,j)) {case 0: Gp1 = (1- tanThead) * E + tanThead * NE;    // 正向梯度Gp2 = (1- tanThead) * W + tanThead * SW;    // 反向梯度break;case 1: Gp1 = (1- tanThead) * N + tanThead * NE;Gp2 = (1- tanThead) * S + tanThead * SW;break;case 2:Gp1 = (1- tanThead) * N + tanThead * NW;Gp2 = (1- tanThead) * S + tanThead * SE;break;case 3:Gp1 = (1- tanThead) * W + tanThead *NW;Gp2 = (1- tanThead) * E + tanThead *SE;break;default:break;}// NMS -非极大值抑制和双阈值检测// 3.1 非极大值抑制if(abs(Gp1) > 255)Gp1 = 255;if(abs(Gp2) > 255)Gp2 = 255;Gp1 = (Gp1 > 0) ? Gp1 : -Gp1;Gp2 = (Gp2 > 0) ? Gp2 : -Gp2;if(grad.at<uchar>(i, j) >= Gp1  && grad.at<uchar>(i, j) >= Gp2){  // 4.2 双阈值检测if(grad.at<uchar>(i, j) >= highValue){ // 高于高阈值： 表明是边缘grad.at<uchar>(i, j) = highValue;result.at<uchar>(i, j) = 255;}else if(grad.at<uchar>(i, j) < lowValue)grad.at<uchar>(i, j) = 0;          // 低于低阈值：抑制掉这个点，它不可能是边缘elsegrad.at<uchar>(i, j) = lowValue;   // 介于两个阈值之间：可能是边缘 }elsegrad.at<uchar>(i, j) = 0;}int m = 1;//5、抑制孤立的弱边缘for(int i = m ; i < grad.rows-m; i++ )for( int j = m; j < grad.cols-m; j++)if(grad.at<uchar>(i, j) == lowValue){result.at<uchar>(i,j) = find8StrongGradient(grad(cv::Rect(j - m, i - m, 3, 3)), highValue);}// 保存结果result_list.push_back(result);}// 3 灰度滤波void grayGaussFiltering(int pic_id){int size = gauss_template.rows;int m = size/2;// 这里的灰度已经领填充过了,但是彩色图像没有cv::Mat gray_image = gray_image_list[pic_id];cv::Mat gauss_image = cv::Mat::zeros(gray_image_list[pic_id].rows, gray_image_list[pic_id].cols, gray_image.type());// 计算高斯滤波for(int i = m; i < gray_image.rows - m; i++)for(int j = m; j < gray_image.cols - m; j++){ cv::Mat sub_matrix = gray_image(cv::Rect(j - m, i - m, size, size));gauss_image.at<uchar>(i - m, j - m) = sub_matrix.dot(gauss_template)/gauss_template_sum;}gray_gauss_image_list.push_back(gauss_image);}// 4 sobel算子计算void sobelDetection(int pic_id){int size = sobel_template[0].rows;int m = size/2;// 经过了高斯滤波,小了一圈cv::Mat gray_image = gray_gauss_image_list[pic_id];for(int i = m; i < gray_image.rows - m; i++)for(int j = m; j < gray_image.cols - m; j++){ cv::Mat image_block = gray_image(cv::Rect(j - m, i - m, size, size));computerSobelResult(image_block, i, j);}}// 4.1 计算梯度强度,方向,区域,角度void computerSobelResult(cv::Mat image_block, int i, int j){float Gx = 0.0, Gy = 0.0, g =0.0;for(int n = 0; n < image_block.rows; n++)for(int m = 0; m < image_block.cols; m++){Gx += float(image_block.at<uchar>(n, m)) * sobel_template[0].at<float>(n, m);Gy += float(image_block.at<uchar>(n, m)) * sobel_template[1].at<float>(n, m);}g = abs(Gx) + abs(Gy);if(g > 255)g = 255;if(g > grad_max)grad_max = g;grad.at<uchar>(i,j) = int(g);gradX.at<float>(i,j) = Gx;gradY.at<float>(i,j) = Gy;thead.at<float>(i,j) = Gy/Gx;if(0 <= thead.at<float>(i, j) && thead.at<float>(i, j) <= (M_PI/4.0))locateGrad.at<float>(i, j) = 0;else if(M_PI/4.0 < thead.at<float>(i,j) && thead.at<float>(i,j) <= (M_PI/2.0))locateGrad.at<float>(i, j) = 1;else if(-M_PI/2.0 <= thead.at<float>(i,j) && thead.at<float>(i,j) <= (-M_PI/4.0))locateGrad.at<float>(i, j) = 2;else if(-M_PI/4.0 < thead.at<float>(i,j) && thead.at<float>(i,j) < 0)locateGrad.at<float>(i, j) = 3;}// 5.1 拉成一维向量void convert2vector(cv::Mat& src, std::vector<int>& nums){for(int i = 0; i < src.rows; i++)for(int j = 0; j < src.cols; j++)nums.push_back(src.at<uchar>(i, j));}// 5.2 向量排序void qsort_vector(std::vector<int>& nums, int begin, int end){if(begin >= end)return;int i = partition(nums, begin, end);//std::cout<<"i :"<<i<<"\n";qsort_vector(nums, begin, i-1);qsort_vector(nums, i+1, end);}// 5.2.1 交换void exchange(std::vector<int>& nums, int a,int b){int c = nums[a];nums[a] = nums[b];nums[b] = c;return;}// 5.2.2 划分int partition(std::vector<int>& nums, int begin, int end){int i = begin, j = end + 1;int x = nums[begin];while (true) {while (nums[++i] < x) {// 向右扫描if (i == end)break;}while (nums[--j] > x) {// 向左扫描if (j == begin)break;}if (i >= j) // 指针相遇，切分位置确定break;exchange(nums, i, j);// 交换左右逆序元素}// 运行到最后：i指向从左往右第一个大于x的元素，j指向从右往左第一个小于x的元素。exchange(nums, begin, j);// 将切分元素放在切分位置return j;}// 6 寻找8邻域的强梯度double find8StrongGradient(cv::Mat image_block, int highValue){for(int i=0; i < 3; i++)for(int j=0; j < 3; j++)if(image_block.at<double>(i,j) == highValue)return 255;return 0;    }// 7 提供外部接口Mat getEdgeDetectionImage(int pic_id, Mat& src) {gray_image_list.push_back(color2Gray(src));if(pic_id==0)initialize(pic_id);Canny(pic_id, 0.9, 0.9);return result_list[pic_id];}private:std::vector<cv::Mat> color_image_list;// 对原始图像进行零填充std::vector<cv::Mat> gray_image_list;std::vector<cv::Mat> result_list;// 高斯滤波cv::Mat gauss_template;float gauss_template_sum;std::vector<cv::Mat> gray_gauss_image_list;// sobel算子计算梯度强度 与方向std::vector<cv::Mat> sobel_template;// 沿X方向,Y方向,(x,y)处的梯度,夹角,所在区域cv::Mat gradX;cv::Mat gradY;cv::Mat grad;cv::Mat thead;cv::Mat locateGrad;int grad_max;
};// 实验九：综合实验，meanshift算法、边缘提取、滤波
class Exp9 {
public:// 1 跟踪代码算法void videoTracking() {cv::VideoCapture cap(0);//exp4 = new Exp4;exp7 = new Exp7();//直方图特征子区间bin的数目，分成16个区间hsize = 16;   //色度 hue 范围float hranges[] = {0,180};const float* phranges = hranges;//摄像头未打开的处理方式if (!cap.isOpened()) {std::cout << "***Could not initialize capturing...***\n";}cv::namedWindow("meanShift", 0);cv::namedWindow("Histogram", 0);cv::namedWindow("Edge", 0);cv::setMouseCallback("meanShift", onMouse, 0);cv::Mat histimg = cv::Mat::zeros(200, 320, CV_8UC3);bool paused = false;int i = 0;while (true) {if (!paused) {cap >> frame;frame.copyTo(tracking_image);image_edge = exp7->getEdgeDetectionImage(i, tracking_image);cvtColor(tracking_image, hsv, cv::COLOR_BGR2HSV);if (trackObject) {int from_to[] = { 0, 0 };hue.create(hsv.size(), hsv.depth());//提取h分量，色度cv::mixChannels(&hsv, 1, &hue, 1, from_to, 1);//画直方图if (trackObject < 0)drawHist(histimg, phranges);//进行meanshiftcalcBackProject(&hue, 1, 0, hist, backproj, &phranges);cv::TermCriteria criteria(cv::TermCriteria::EPS | cv::TermCriteria::COUNT, 10, 1);cv::meanShift(backproj, trackWindow, criteria);cv::rectangle(tracking_image, trackWindow, cv::Scalar(255, 0, 0), 3);cv::rectangle(image_edge, trackWindow, cv::Scalar(255, 0, 0), 3);}}if (selectObject && selection.width > 0 && selection.height > 0) {cv::Mat roi(tracking_image, selection);cv::bitwise_not(roi, roi);}cv::imshow("meanShift", tracking_image);cv::imshow("Histogram", histimg);cv::imshow("Edge", image_edge);i++;if (cv::waitKey(10) == 27)break;}}// 2 绘制跟踪目标直方图void drawHist(Mat& histimg, const float* phranges) {cv::Mat roi(hue, selection);cv::calcHist(&roi, 1, 0, cv::Mat(), hist, 1, &hsize, &phranges);// 直方图归一化到0-255cv::normalize(hist, hist, 0, 255, cv::NORM_MINMAX);trackWindow = selection;// 置1，除非重新框选目标，否则只执行一次trackObject = 1;// 绘制背景底色histimg = cv::Scalar::all(0);//320/16；每个区间可以占用histimg20个列（column）int binW = histimg.cols / hsize;//定义一个缓冲单bin矩阵，1行16列cv::Mat buf(1, hsize, CV_8UC3);//Vec3b为3个char值的向量，存放颜色数据for (int i = 0; i < hsize; i++)buf.at<cv::Vec3b>(i) = cv::Vec3b(cv::saturate_cast<uchar>(i*180. / hsize), 255, 255);cv::cvtColor(buf, buf, cv::COLOR_HSV2BGR);// 绘制直方图for (int i = 0; i < hsize; i++) {//获取直方图高度，根据histimg图像的高换算成0到histimg.rows的值int val = cv::saturate_cast<int>(hist.at<float>(i)*histimg.rows / 255);//画出直方图，左上角坐标，右下角坐标，高度，颜色，大小，线型//rows-va1是因为矩阵是左上角开始数的，y轴向下，而要画的图是y轴向上cv::rectangle(histimg, cv::Point(i*binW, histimg.rows),cv::Point((i + 1)*binW, histimg.rows - val),cv::Scalar(buf.at<cv::Vec3b>(i)), -1, 8);}}private:// 当前帧Mat frame;// hsv图像Mat hsv;// 色调Mat hue;// 直方图信息Mat hist;// 反向投影Mat backproj;Mat image_noise;Mat image_process;Mat image_edge;cv::Rect trackWindow;Exp7* exp7;int hsize;};int main()
{Exp9 a;a.videoTracking();
}

三、实验结果