1原理

https://www.bilibili.com/video/av74302620/?spm_id_from=333.788.videocard.0

https://blog.csdn.net/fzp95/article/details/78385795?utm_medium=distribute.pc_relevant.none-task-blog-baidujs_title-4&spm=1001.2101.3001.4242

相关和卷积操作》
https://blog.csdn.net/qq_17783559/article/details/82254996

1.1 相关性

在这里插入图片描述
其中 f∗表示 f的复共轭。correlation的直观解释就是衡量两个函数在某个时刻相似程度。

1.2 MOSSE 滤波器H公式推导

作者提出的滤波器叫做Minimum Output Sum of Squared Error filter（MOSSE）（误差最小平方和滤波器）。

MOSSE是一种从较少的训练图像中生成滤波器的算法。

1.2.1 结论

先放结论：
整个滤波器的模型公式：在这里插入图片描述
分子是输入和期望输出之间的相关性，分母是输入的能谱。
其中gi表示响应输出，fi表示输入图像，对多组输入和对应的输出求均值，得到h，h表示滤波模板。

Q1: 对于框定的一个目标框，怎么得到多个训练样本输入？
A1: 可以通过任意旋转裁剪

Q2: 训练时候对应的响应是什么？
A2: 一般来说，gi 可以是任何形状。在这种情况下，gi是由ground truth生成的，它在训练图像fi中有一个以目标为中心的紧凑的(σ = 2.0)二维高斯峰.

Q3：经过样本训练得到了H ，之后怎么用？
A3：利用在这里插入图片描述得到相关性最大的，即响应最大的地方，作为下一次跟踪的目标位置。

ps

在实际跟踪的过程中我们要考虑到目标的外观变换等因素的影响，所以需要同时考虑目标的m个图像作为参考，对跟踪框（groundtruth）进行随机仿射变换，获取一系列的训练样本 fi，

而 gi则是由高斯函数产生，并且其峰值位置是在 fi的中心位置。

获得了一系列的训练样本和结果之后，就可以计算滤波器h的值。注意：这里的f,g,h的size大小都相同。从而提高滤波器模板的鲁棒性。

同时，作者为了让滤波器对与形变、光照等外界影响具有更好的鲁棒性，采取了如下的模板更新策略,作者将滤波器的模型公式分为分子和分母两个部分，每个部分都分别的进行更新:。
在这里插入图片描述
其中 At和 At−1分别表示的是当前帧和上一帧的分子。

1.2.2 推导（可以不管）

我们需要找到一个滤波器，使其在目标上的响应最大，则如下公式：
在这里插入图片描述
其中g表示响应输出，f表示输入图像，h表示滤波模板。

显然，我们要是想获得比较获得响应输出，只需确定滤波器模板h即可。上式的计算是进行卷积计算，这在计算机中的计算消耗时很大的，因此作者对上式进行快速傅里叶变换（FFT）这样卷积操作经过FFT后就变成了点乘操作，极大的减少了计算量。上式变成如下形式：

在这里插入图片描述
为了方便描述，将上式写成如下形式：

后面跟踪的任务就是找到 H∗了

因为上式的操作都是元素级别的，因此要想找到，只要使其中的每个元素的MOSSE都最小即可。因此上式可转换为如下形式：
在这里插入图片描述（w和v是H中每个元素的索引）

要想得到最小的 H*wv,只需要对上式求偏导，并使偏导为0即可。即
在这里插入图片描述
即每个位置的值为：

结论

整个滤波器的模型公式：在这里插入图片描述
分子是输入和期望输出之间的相关性，分母是输入的能谱。

ps

但是在实际跟踪的过程中我们要考虑到目标的外观变换等因素的影响，所以需要同时考虑目标的m个图像作为参考，对跟踪框（groundtruth）进行随机仿射变换，获取一系列的训练样本 fi，

而 gi则是由高斯函数产生，并且其峰值位置是在 fi的中心位置。

获得了一系列的训练样本和结果之后，就可以计算滤波器h的值。注意：这里的f,g,h的size大小都相同。从而提高滤波器模板的鲁棒性。

2 具体计算过程

首先，它需要一组训练图像fiand训练输出gi。

2.1 输入目标区域

rect

%% 选择目标区域
% select target from first frame
im = imread(img_files(1,:));
f = figure('Name', 'Select object to track'); 
imshow(im);
rect = getrect;
close(f); clear f;
center = [rect(2)+rect(4)/2 rect(1)+rect(3)/2];

2.2 获得高斯核

%% 高斯核
% plot gaussian
sigma = 100;
gsize = size(im);%获取图像尺寸
[R,C] = ndgrid(1:gsize(1), 1:gsize(2));
%两个矩阵R和C，都是m行n列
g = gaussC(R,C, sigma, center);%通过R和C产生size等同于im的高斯滤波函数
g = mat2gray(g);
g = imcrop(g, rect);
G = fft2(g);    % 得到目标区域的高斯滤波核的傅里叶变换

2.3 计算输入的fi

其中fi是框选的目标区域

img = imcrop(img, rect);

人为地连接图像fi的边界会引入一个影响相关输出的伪影。

解决方法：

首先，使用日志函数对像素值进行转换，这有助于低对比度的照明情况。像素值被归一化为平均值0.0和平均值1.0。

最后，将图像乘以一个余弦窗口，使边缘附近的像素值逐渐减少到零。这也有一个好处，它把重点放在了靠近目标中心的地方。
原因：由于傅里叶变换是周期性的，它不尊重图像的边界。非周期图像的相对边之间的巨大不连续将导致有噪声的傅里叶表示。

%将高斯滤波函数变换到频域
height = size(g,1);
width = size(g,2);
%imresize(img, [height width])将图片调整成滤波器的大小
fi = preprocess(imresize(img, [height width]));

preprocess.m处理函数

function img = preprocess(img)
[r,c] = size(img);
win = window2(r,c,@hann);
save win
eps = 1e-5;
img = log(double(img)+1);
img = (img-mean(img(:)))/(std(img(:))+eps);
img = img.*win;
end

2.4 计算初始的Hi(Ai,Bi)及生成多个训练模板

Ai = (G.*conj(fft2(fi)));%相关性
Bi = (fft2(fi).*conj(fft2(fi)));%fi能量谱N = 128;
for i = 1:Nfi = preprocess(rand_warp(img));Ai = Ai + (G.*conj(fft2(fi)));Bi = Bi + (fft2(fi).*conj(fft2(fi)));
end

2.5 利用初始Hi进行预测响应位置，并在线更新

%预测位置Hi = Ai./Bi;fi = imcrop(img, rect);            fi = preprocess(imresize(fi, [height width])); gi = uint8(255*mat2gray(ifft2(Hi.*fft2(fi))));maxval = max(gi(:))[P, Q] = find(gi == maxval);dx = mean(P)-height/2;dy = mean(Q)-width/2;rect = [rect(1)+dy rect(2)+dx width height];%更新滤波器fi = imcrop(img, rect); fi = preprocess(imresize(fi, [height width]));Ai = eta.*(G.*conj(fft2(fi))) + (1-eta).*Ai;Bi = eta.*(fft2(fi).*conj(fft2(fi))) + (1-eta).*Bi;

code

matlab

%get images from source directory
%% 获取图像
%此处仅仅用于得到图片序列所在地址
datadir = 'H:/code/1doctor_code/matlab_code/tracking/mosse-tracker-master/data/';
dataset = 'Surfer';
path = [datadir dataset];
img_path = [path '/img/'];
D = dir([img_path, '*.jpg']);%在img_path下的所有jpg后缀文件的地址放入D中seq_len = length(D(not([D.isdir])));%得到图片总数
if exist([img_path num2str(1, '%04i.jpg')], 'file')img_files = num2str((1:seq_len)', [img_path '%04i.jpg']);
elseerror('No image files found in the directory.');
end
%% 选择目标区域
% select target from first frame
im = imread(img_files(1,:));
f = figure('Name', 'Select object to track'); 
imshow(im);
rect = getrect;
close(f); clear f;
center = [rect(2)+rect(4)/2 rect(1)+rect(3)/2];
%% 高斯核
% plot gaussian
sigma = 100;
gsize = size(im);%获取图像尺寸
[R,C] = ndgrid(1:gsize(1), 1:gsize(2));
%两个矩阵R和C，都是m行n列
g = gaussC(R,C, sigma, center);%通过R和C产生size等同于im的高斯滤波函数
g = mat2gray(g);
g = imcrop(g, rect);
G = fft2(g);    % 得到目标区域的高斯滤波核的傅里叶变换
%% 随机裁剪原始目标区域，得到N个训练样本
%randomly warp original image to create training set
if (size(im,3) == 3) img = rgb2gray(im); 
end
img = imcrop(img, rect);%将高斯滤波函数变换到频域
height = size(g,1);
width = size(g,2);
%imresize(img, [height width])将图片调整成滤波器的大小
fi = preprocess(imresize(img, [height width]));
Ai = (G.*conj(fft2(fi)));
Bi = (fft2(fi).*conj(fft2(fi)));N = 128;
for i = 1:Nfi = preprocess(rand_warp(img));Ai = Ai + (G.*conj(fft2(fi)));Bi = Bi + (fft2(fi).*conj(fft2(fi)));
end% MOSSE online training regimen
eta = 0.25;
fig = figure('Name', 'MOSSE');
t = figure;
mkdir(['results_' dataset]);
for i = 1:size(img_files, 1)img = imread(img_files(i,:));im = img;if (size(img,3) == 3)img = rgb2gray(img);%灰度图endif (i == 1)Ai = eta.*Ai;Bi = eta.*Bi;%似乎没有意义，在i=2的时候同时约了etaelse%预测位置Hi = Ai./Bi;fi = imcrop(img, rect);            fi = preprocess(imresize(fi, [height width])); gi = uint8(255*mat2gray(ifft2(Hi.*fft2(fi))));maxval = max(gi(:))[P, Q] = find(gi == maxval);dx = mean(P)-height/2;dy = mean(Q)-width/2;rect = [rect(1)+dy rect(2)+dx width height];%更新滤波器fi = imcrop(img, rect); fi = preprocess(imresize(fi, [height width]));Ai = eta.*(G.*conj(fft2(fi))) + (1-eta).*Ai;Bi = eta.*(fft2(fi).*conj(fft2(fi))) + (1-eta).*Bi;end% visualizationtext_str = ['Frame: ' num2str(i)];box_color = 'green';position=[1 1];result = insertText(im, position,text_str,'FontSize',15,'BoxColor',...box_color,'BoxOpacity',0.4,'TextColor','white');result = insertShape(result, 'Rectangle', rect, 'LineWidth', 3);imwrite(result, ['results_' dataset num2str(i, '/%04i.jpg')]);imshow(result);drawnow;rect
end

c

// This file ispart of the OpenCV project. 
// It is subject to the license terms in the LICENSEfile found in the top-level directory 
// of this distribution and athttp://opencv.org/license.html. // 
//[1] David S. Bolme et al. "Visual Object Trackingusing Adaptive Correlation Filters" 
//   http://www.cs.colostate.edu/~draper/papers/bolme_cvpr10.pdf 
// // 
// credits: 
// Kun-Hsin Chen: for initial c++ code 
// Cracki: for the idea of only converting the usedpatch to gray 
// #include "opencv2/tracking.hpp" namespace cv { 
namespace tracking { struct DummyModel :TrackerModel 
{ virtual void modelUpdateImpl(){} virtual void modelEstimationImpl( const std::vector<Mat>& ){} 
}; const double eps=0.00001;      // fornormalization 
const double rate=0.2;         //learning rate 
const double psrThreshold=5.7; //no detection, if PSR is smaller than this struct MosseImpl :TrackerMOSSE 
{ 
protected: Point2d center;//center of the bounding box Size size;      //size ofthe bounding box Mat hanWin; Mat G;          //goal Mat H, A,B;    //state//  Element-wisedivision of complex numbers in src1 and src2 Mat 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 forboth 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 preProcess( Mat &window) const { window.convertTo(window, CV_32F); log(window+ 1.0f, window); //normalize Scalarmean,StdDev; meanStdDev(window, mean, StdDev); window =(window-mean[0]) /(StdDev[0]+eps); //Gaussain weighting window =window.mul(hanWin); } double correlate( const Mat &image_sub, Point&delta_xy ) const//计算相对位移 { MatIMAGE_SUB, RESPONSE, response; // filter in dft space dft(image_sub, IMAGE_SUB, DFT_COMPLEX_OUTPUT); mulSpectrums(IMAGE_SUB, H, RESPONSE, 0, true ); idft(RESPONSE, response, DFT_SCALE|DFT_REAL_OUTPUT); // update center position double maxVal; Point maxLoc; minMaxLoc(response, 0, &maxVal, 0, &maxLoc); delta_xy.x= maxLoc.x - int(response.size().width/2); delta_xy.y= maxLoc.y - int(response.size().height/2); // normalize response Scalarmean,std; meanStdDev(response, mean, std); return (maxVal-mean[0]) / (std[0]+eps); // PSR } Mat randWarp( const Mat& a ) const { cv::RNGrng(8031965); // random rotation double C=0.1; double ang = rng.uniform(-C,C); double c=cos(ang), s=sin(ang); // affine warp matrix 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(), BORDER_REFLECT); return warped; } virtual bool initImpl( const Mat& image, const Rect2d& boundingBox ){ model =makePtr<DummyModel>(); Mat img; if (image.channels() == 1) img =image; else cvtColor(image, img, COLOR_BGR2GRAY); int w = getOptimalDFTSize(int(boundingBox.width)); int h = getOptimalDFTSize(int(boundingBox.height)); //Get the center position int x1 = int(floor((2*boundingBox.x+boundingBox.width-w)/2)); int y1 = int(floor((2*boundingBox.y+boundingBox.height-h)/2)); center.x =x1 + (w)/2; center.y =y1 + (h)/2; size.width= w; size.height= h; Mat window;getRectSubPix(img, size, center, window); createHanningWindow(hanWin, size, CV_32F); // goal Matg=Mat::zeros(size,CV_32F); g.at<float>(h/2, w/2) = 1;GaussianBlur(g, g, Size(-1,-1),2.0); double maxVal; minMaxLoc(g, 0,&maxVal); g = g /maxVal; dft(g, G,DFT_COMPLEX_OUTPUT); // initial A,B and H A =Mat::zeros(G.size(), G.type()); B =Mat::zeros(G.size(), G.type()); for(int i=0; i<8; i++) { Matwindow_warp = randWarp(window); preProcess(window_warp); MatWINDOW_WARP, A_i, B_i; dft(window_warp, WINDOW_WARP, DFT_COMPLEX_OUTPUT); mulSpectrums(G          ,WINDOW_WARP, A_i, 0, true); mulSpectrums(WINDOW_WARP, WINDOW_WARP, B_i, 0, true); A+=A_i;B+=B_i;} H =divDFTs(A,B); return true; } virtual bool updateImpl( const Mat& image,Rect2d& boundingBox ) { if (H.empty()) // not initialized return false; Matimage_sub; getRectSubPix(image, size, center, image_sub); if (image_sub.channels() != 1) cvtColor(image_sub, image_sub, COLOR_BGR2GRAY); preProcess(image_sub); Pointdelta_xy; double PSR =correlate(image_sub, delta_xy); if (PSR < psrThreshold) return false; //update location center.x +=delta_xy.x; center.y +=delta_xy.y; Matimg_sub_new; getRectSubPix(image, size, center, img_sub_new); if (img_sub_new.channels() !=1) cvtColor(img_sub_new, img_sub_new, COLOR_BGR2GRAY); preProcess(img_sub_new); // new state for A and B Mat F, A_new,B_new; dft(img_sub_new, F, DFT_COMPLEX_OUTPUT); mulSpectrums(G, F, A_new, 0, true ); mulSpectrums(F, F, B_new, 0, true ); // update A ,B, and H A = A*(1-rate) + A_new*rate; B = B*(1-rate) + B_new*rate; H =divDFTs(A, B); // return tracked rect double x=center.x, y=center.y; int w = size.width,h=size.height; boundingBox= Rect2d(Point2d(x-0.5*w,y-0.5*h), Point2d(x+0.5*w, y+0.5*h)); return true; } public: MosseImpl() {isInit = 0; } // dummy implementation. virtual void read( const FileNode& ){} virtual void write( FileStorage& ) const{} }; // MosseImpl } // tracking Ptr<TrackerMOSSE> TrackerMOSSE::create() 
{ returnmakePtr<tracking::MosseImpl>(); 
} } // cv