mosse（Minimizing the Output Sum of Squared Error）代码解读

---

mosse论文解读：https://blog.csdn.net/fjswcjswzy/article/details/106172434

---

mosse作为使用相关滤波进行目标跟踪的鼻祖，必须得把这玩意弄明白才能进行接下来的研究，我就从代码的角度，一行一行理解代码，保证弄得一清二楚。此文的代码为matlab，其他语言类似

步骤一： 读取视频帧

首先，要做的事将视频帧读入，代码应该比较好理解

datadir = '../data/';
dataset = 'Surfer';
path = [datadir dataset];
img_path = [path '/img/'];
D = dir([img_path, '*.jpg']);
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

步骤二：在第一帧中选择要跟踪的目标

在把视频帧读入后，就要在第一帧中用矩形框 选出我们要跟踪的目标，并且找出所选矩形框的中心center

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];

这里我们的第一帧和要跟踪的目标如下图：

步骤三：画出高斯图像

画出高斯的代码为：

sigma = 100;
gsize = size(im);
[R,C] = ndgrid(1:gsize(1), 1:gsize(2));
g = gaussC(R,C, sigma, center);
g = mat2gray(g);

这里高斯滤波器的函数为：

function val = gaussC(x, y, sigma, center)
    xc = center(1);
    yc = center(2);
    exponent = ((x-xc).^2 + (y-yc).^2)./(2*sigma);
    val       = (exp(-exponent)); 

• 这里我们使用的图片大小是360x480 ，所以R和C都是大小为360x480的矩阵，

• 将R,C, sigma,center当做参数传入gaussC函数，注意这里的高斯做了一些改变，正常来说二维高斯曲面的公式（x,y代表像素的模板坐标，模板中心位置为原点）为：

  $$\mathbf{g}\left(\mathbf{x},\mathbf{y}\right)=\frac{1}{2{\mathbf{\pi }\mathbf{\sigma }}^2}{\mathbf{e}}^{-\frac{\left({\mathbf{x}}^2+{\mathbf{y}}^2\right)}{2{\mathbf{\sigma }}^2}}$$ 但是为什么代码中，只有${\mathbf{e}}^{-\frac{\left({\mathbf{x}}^2+{\mathbf{y}}^2\right)}{2{\mathbf{\sigma }}^2}}$这一部分呢？------这是因为我们要画出高斯图像 ，使高斯响应在我们所框的框框的中心达到最大的峰值，而$\frac{1}{2{\mathbf{\pi }\mathbf{\sigma }}^2}$里都是常数，对我们找高斯位置并没有影响，所以把其舍去。

• 再者，(x-xc).^2 + (y-yc).^2，为什么要减去xc和yc呢？-------这是为了能让峰值中心移动到我们所框出的框框的中心位置，因为我们如果没有减去xc和yc，那么峰值中心就在原点位置，如下图：
  但是，如果减去xc和yc，那么就如下图：
  

• 这里g = mat2gray(g)实现图像矩阵的归一化操作。所谓”归一化”就是使矩阵的每个元素的值都在0和1之间。

步骤四：对第一帧图像得到一个初始滤波器

if (size(im,3) == 3) 
    img = rgb2gray(im); 
end
img = imcrop(img, rect);
g = imcrop(g, rect);
G = fft2(g);

height = size(g,1);
width = size(g,2);
fi = preprocess(imresize(img, [height width]));

Ai = (G.*conj(fft2(fi)));

Bi = (fft2(fi).*conj(fft2(fi)));

这里预处理函数preprocess为：

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

• 首先将彩色图片转成灰度图像
• 再将原始图片和对原始图片所画出的高斯图像 按我们所框的大小和位置切割下来
• 再对所切割下来的高斯图像做二维傅里叶变换 并将其赋值给G
• 再对所切割下来的原始图片进行预处理，所得的结果赋予fi
• 预处理的作用是：使用对数函数对像素值进行转换，该函数有助于解决低对比度照明情况。 像素值被归一化为平均值为0.0，范数为1.0。 最后，将图像乘以余弦窗，该余弦窗将边缘附近的像素值逐渐减小为零。 这还有一个好处，就是可以将更多的重点放在目标的中心附近。
• 这时候可以得到得到Ai 和 Bi，分别对应论文中的Ai 和 Bi，如下图，从论文公式可以得出，当前滤波器是与前一帧的滤波器有关的，当第一帧时，其实就是求这个滤波器的初始值，$\mathbf{\eta }$为1：
  
• 当然这个初始滤波器是非常理想化的，相当于只有一个样本，对于一个样本，结果是过拟合的，为了去除这种过拟合，我们需要更多的样本，该怎么做呢？-----对图片进行旋转

第五步：随机旋转框内的图像，并且得到训练集

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

这里rand_warp(img)函数为：

function img = rand_warp(img)
    a = -180/16;
    b = 180/16;
    r = a + (b-a).*rand;
    sz = size(img);
    scale = 1-0.1 + 0.2.*rand;
    % trans_scale = randi([-4,4], 1, 1);
    img = imresize(imresize(imrotate(img, r), scale), [sz(1) sz(2)]);

• 这里对所框出的图像进行随机旋转变换（变换幅度比较小），这样就产生了128张训练集样本，部分样本如下图：
• 将新的Ai和Bi加到原始的Ai和Bi上，相当于对原始的Ai和Bi进行一个修正

步骤六：在线训练

在线训练的意思就是根据每一帧得到的结果进行修正反馈更新。代码如下：

eta = 0.125;
fig = figure('Name', 'MOSSE');
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;
    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);
end

• 论文中给出的公式是：
  

• 这里eta就是$\mathbf{\eta }$，取值为0~1，代表着当前帧受前一帧影响的大小，若$\mathbf{\eta }$越小，则表示受前一帧影响就越大

• 当第一帧进入时，我们用第一幅图来训练，所以Ai和Bi不用更新：Ai = eta.*Ai; Bi = eta.*Bi;

• 将第一帧得到的滤波器Hi，用在第二帧上

• gi是什么呢？先知道Gi是对Hi和Fi进行点乘得到的，在对点乘的结果Gi进行反傅里叶变换，就得到gi，在论文中的公式为：
  

• 将Gi变换到gi，是由频域变换到空间域，空间域就是：

• 接下来就要寻找gi的最大值所在的坐标，这里使用mean（P）和mean（Q）是为什么呢？------说明取得最大值的地方有很多，这里取了一个平均值，这里如果是一个比较简单的场景，比如一辆飞机在蓝色的天空中飞过，可能就只有一个最大值，但是如果是一个复杂场景，比如一个人在马路上走，可能就有很多响应的最大值，所以取了一个平均值

• 再在取最大值的点的位置，画出一个与我们所画的矩形框有相同大小的矩形框（通过这个新画出来矩形框继续去运行，直到迭代完全部的视频帧）

• 最后再根据我们的论文，更新Ai和Bi

• 最后再将所更新的矩形框画出来