【opencv】目标跟踪之MOSSE算法配合模板匹配实现初始滤波器的自动初始化

因为最近有一个需求是识别视频流中的指定物体，一些轻量级的算法实际测试效果都不太好，所以考虑到目标追踪算法。（是自动识别！！！！！！）

Opencv八种目标追踪算法：

BOOSTING Tracker：和Haar cascades（AdaBoost）背后所用的机器学习算法相同，但是距其诞生已有十多年了。这一追踪器速度较慢，并且表现不好，但是作为元老还是有必要提及的。（最低支持OpenCV 3.0.0）

MIL Tracker：比上一个追踪器更精确，但是失败率比较高。（最低支持OpenCV 3.0.0）

KCF Tracker：比BOOSTING和MIL都快，但是在有遮挡的情况下表现不佳。（最低支持OpenCV 3.1.0）

CSRT Tracker：比KCF稍精确，但速度不如后者。（最低支持OpenCV 3.4.2）

MedianFlow Tracker：在报错方面表现得很好，但是对于快速跳动或快速移动的物体，模型会失效。（最低支持OpenCV 3.0.0）

TLD Tracker：我不确定是不是OpenCV和TLD有什么不兼容的问题，但是TLD的误报非常多，所以不推荐。（最低支持OpenCV 3.0.0）

MOSSE Tracker：速度真心快，但是不如CSRT和KCF的准确率那么高，如果追求速度选它准没错。（最低支持OpenCV 3.4.1）

GOTURN Tracker：这是OpenCV中唯一一深度学习为基础的目标检测器。它需要额外的模型才能运行，本文不详细讲解。（最低支持OpenCV 3.2.0）

建议：

如果追求高准确度，又能忍受慢一些的速度，那么就用CSRT

如果对准确度的要求不苛刻，想追求速度，那么就选KCF

纯粹想节省时间就用MOSSE

因为MOSSE算法可以说是相关滤波目标追踪的开山之作，可以先通过学习MOSSE算法进而学习KCF。

对于MOSSE算法目前我仅限于把源码总体阅读一遍，其原理的讲解我就不班门弄斧了hhhhhhh。

经验参考自: http://www.elecfans.com/d/722414.html.

下面进入正题：

下面是通过对想识别的物体制作模板，以此模板对视频帧进行模板匹配，获取视频帧中该模板的位置，以此来初始化滤波器，从而正式进入MOSSE算法的工作。

先上效果图：
1.模板

2.视频流第一帧进行模板匹配的效果

3.目标追踪效果：（怎么上传视频！！！）

效果还可以

放代码：
tracking.h

#include <opencv2/opencv.hpp>

#include <sstream>

class Tracking{
   
   
public:
    cv::Mat divDFTs( const cv::Mat &src1, const cv::Mat &src2 ) const;
    void preProcess( cv::Mat &window ) const;
    double correlate( const cv::Mat &image_sub, cv::Point &delta_xy ) ;
    cv::Mat randWarp( const cv::Mat& a ) const;
    bool initImpl( const cv::Mat& image, const cv::Rect2d& boundingBox );
    bool updateImpl( const cv::Mat& image, cv::Rect2d& boundingBox );
    cv::Mat VisualDft(cv::Mat input);


public:
    cv::Point2d center; //center of the bounding box
    cv::Size size;      //size of the bounding box
    cv::Mat hanWin;
    cv::Mat g,G;          //goal
    cv::Mat H, A,B;    //state
};

tracking.cpp

#include "tracking.h"
using namespace cv;
using namespace std;

const double eps=0.00001;      // for normalization
const double rate=0.2;         // learning rate
const double psrThreshold=1.7; // no detection, if PSR is smaller than this

Mat Tracking::VisualDft(cv::Mat input) {
   
   
    Mat padded;                 //以0填充输入图像矩阵
    int m = getOptimalDFTSize(input.rows); //getOptimalDFTSize函数返回给定向量尺寸的傅里叶最优尺寸大小。
    int n = getOptimalDFTSize(input.cols);

    //填充输入图像I，输入矩阵为padded，上方和左方不做填充处理
    copyMakeBorder(input, padded, 0, m - input.rows, 0, n - input.cols, BORDER_CONSTANT, Scalar::all(0));//四个方向的常量0扩充

    Mat planes[] = {
   
    Mat_<float>(padded), Mat::zeros(padded.size(),CV_32F) };//两个矩阵
    Mat complexI;
    merge(planes, 2, complexI);     //将planes融合合并成一个多通道数组complexI

    dft(complexI, complexI);        //进行傅里叶变换

    //计算幅值，转换到对数尺度(logarithmic scale)
    //=> log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
    split(complexI, planes);        //planes[0] = Re(DFT(I),planes[1] = Im(DFT(I))
    //即planes[0]为实部,planes[1]为虚部
    magnitude(planes[0], planes[1], planes[0]);     //planes[0] = magnitude,求幅度谱
    Mat magI = planes[0];

    magI += Scalar::all(1);
    log(magI, magI);                //转换到对数尺度(logarithmic scale)

    //如果有奇数行或列，则对频谱进行裁剪
    magI = magI(Rect(0, 0, magI.cols&-2, magI.rows&-2));  //magI.rows&-2得到不大于magI.rows的最大偶数

    //重新排列傅里叶图像中的象限，使得原点位于图像中心
    int cx = magI.cols / 2;
    int cy = magI.rows / 2;

    Mat q0(magI, Rect(0, 0, cx, cy));       //左上角图像划定ROI区域
    Mat q1(magI, Rect(cx, 0, cx, cy));      //右上角图像
    Mat q2(magI, Rect(0, cy, cx, cy));      //左下角图像
    Mat q3(magI, Rect(cx, cy, cx, cy));     //右下角图像

    //变换左上角和右下角象限
    Mat tmp;
    q0.copyTo(tmp);
    q3.copyTo(q0);
    tmp.copyTo(q3);

    //变换右上角和左下角象限
    q1.copyTo(tmp);
    q2.copyTo(q1);
    tmp.copyTo(q2);

    //归一化处理，用0-1之间的浮点数将矩阵变换为可视的图像格式
    normalize(magI, magI, 0, 1, CV_MINMAX);

    return magI;
}//图像的傅利叶变换。

Mat Tracking::divDFTs( const Mat &src1, const Mat &src2 ) const{
   
   
    Mat c1[2],c2[2],a1,a2,s1,s2,denom,re,im;

    // split into re and im per src
    cv::split(src1, c1);
    cv::split(src2, c2);

    // (Re2*Re2 + Im2*Im2) = denom
    //   denom is same for both channels
    cv::multiply(c2[0], c2[0], s1);
    cv::multiply(c2[1], c2[1], s2);
    cv::add(s1, s2, denom);

    // (Re1*Re2 + Im1*Im1)/(Re2*Re2 + Im2*Im2) = Re
    cv::multiply(c1[0], c2[0], a1);
    cv::multiply(c1[1], c2[1], a2);
    cv::divide(a1+a2, denom, re, 1.0 );

    // (Im1*Re2 - Re1*Im2)/(Re2*Re2 + Im2*Im2) = Im
    cv::multiply(c1[1], c2[0], a1);
    cv::multiply(c1[0], c2[1], a2);
    cv::divide(a1+a2, denom, im, -1.0);

    // Merge Re and Im back into a complex matrix
    Mat dst, chn[] = {
   
   re,im};
    cv::merge(chn, 2, dst);
    return dst;
}

void Tracking::preProcess( Mat &window ) const{
   
   
    imshow("灰度图片",window);
    moveWindow("灰度图片",20,20);
    window.convertTo(window, CV_32F);//转化为32位浮点型
    log(window + 1.0f, window);     //对灰度图像进行对数变换，增强低光照情况下的对比度
    //normalize
    Scalar mean,StdDev;//val[4]   各个通道。
    meanStdDev(window, mean, StdDev);       //计算得到均值和标准差
    window = (window-mean[0]) / (StdDev[0]+eps); //图像归一化
    imshow("经过对数变换后",window);
    moveWindow("经过对数变换后",220,20);
    //Gaussain weighting
    window = window.mul(hanWin); //空域对图像加汉宁窗，突出图像中心区域，弱化周围背景
    imshow("加汉宁窗后",window);
    moveWindow("加汉宁窗后",420,20);
}

double Tracking::correlate( const Mat &image_sub, Point &delta_xy )
{
   
   
    Mat IMAGE_SUB, RESPONSE, response;
    // filter in dft space
    dft(image_sub, IMAGE_SUB, DFT_COMPLEX_OUTPUT); //对图像做dft，转化到频域,生成共轭对称矩阵，分别储存实部和虚部
    mulSpectrums(IMAGE_SUB, H, RESPONSE, 0, true ); //时域卷积，频域相乘（矩阵点乘），矩阵中对应元素相乘
    idft(RESPONSE, response, DFT_SCALE|DFT_REAL_OUTPUT); //离散傅里叶反变换，返回到时域
    Mat response2show = VisualDft(response);//求实际响应
    imshow("实际输出响应",response2show);
    moveWindow("实际输出响应",1220,20);
    // update center position
    double maxVal; Point maxLoc;
    minMaxLoc(response, 0, &maxVal, 0, &maxLoc);//寻找输出相关中的响应最大值位置，即为下一帧物体中心
    delta_xy.x = maxLoc.x - int(response.size().width/2); //获得x方向运动偏移
    delta_xy.y = maxLoc.y - int(response.size().height/2);//获得y方向运动偏移
    // normalize response
    Scalar mean,std;
    meanStdDev(response, mean, std);//mean为均值，std为标准差
    cout<<"psr is :"<<(maxVal-mean[0]) / (std[0]+eps)<<endl;
    return (maxVal-mean[0]) / (std[0]+eps); // PSR返回波峰旁瓣比
}

Mat Tracking::randWarp( const Mat& a ) const
{
   
   
    cv::RNG rng(8031965);

    // random rotation
    double C=0.1;
    double ang = rng.uniform(-C,C);
    double c=cos(ang), s=sin(ang);
    // 生成随机仿射矩阵，只产生微小形变，不产生位移。防止单一训练样本带来的过拟合，在初始化时引入训练样本，增加滤波器的鲁棒性。
    Mat_<float> W(2,3);
    W << c + rng.uniform(-C,C), -s + rng.uniform(-C,C), 0,
            s + rng.uniform(-C,C),  c + rng.uniform(-C,C), 0;

    // random translation
    Mat_<float> center_warp(2, 1);
    center_warp << a.cols/2, a.rows/2;
    W.col(2) = center_warp - (W.colRange(0, 2))*center_warp;

    Mat warped;
    warpAffine(a, warped, W, a.size