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大小都相同。从而提高滤波器模板的鲁棒性。
同时，作者为了让滤波器对与形变、光照等外界影响具有更好的鲁棒性，采取了如下的模板更新策略,作者将滤波器的模型公式分为分子和分母两个部分，每个部分都分别的进行更新:。
其中 At和 At−1分别表示的是当前帧和上一帧的分子。
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:N
fi = 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']);
else
error('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:N
fi = 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);%灰度图
end
if (i == 1)
Ai = eta.*Ai;
Bi = eta.*Bi;%似乎没有意义，在i=2的时候同时约了eta
else
%预测位置
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
% visualization
text_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
1292
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
