相关滤波目标追踪一：鼻祖Mosse算法的matlab解析

一些自己的见解，如有不同观点，可以一起讨论。

github地址：https://github.com/Ronales/Mosse_Tracking_matlab

补充：说一下mosse的更新策略：

1. "论文目标就是找到一个滤波器h，使其在输入图像上，跟踪的目标上的位置响应值最大。f表示训练图像，g表示响应输出，h表示滤波器，F，G，H对应其频域值(均为傅里叶变换后的值)。"这句话谁都能懂，落到实处就是：

（1）初始状态（第一帧），我们首先框选一副图像的某一个人的位置，则按照人体目标框选位置裁剪出来的局部图像就是f有了。

（2）由于mosse是一个在线更新的阶段，即响应输出g得到的坐标位置应当作为下一帧更新的输入提供位置裁切参考。这样一来初始状态没有同时具有f和g的条件。（初值状态只有手动框选的一个目标位置），则就无法初始化滤波器h。

这里需要将裁剪出来的局部图像f进行仿射变换，得到128张图像fi。模拟不稳定状态下的请情况，此时这128次仿射变化得到的响应图g最终结果应当为手动框选的那一张局部图f产生的响应值。这样在初始状态下，就能极大情况的模拟各种情况滤波器的情况，构造一个较强的h，以便后续帧进行在线更新。

############################## 响应图的解释：##############################

（3）按照公式，我们还需要一个响应输出G（傅里叶变换g后得到）。这个G怎么来呢。在作者的代码中的实现为：

center = [rect(2)+rect(4)/2 rect(1)+rect(3)/2];%获得矩形框的中心点纵，横坐标（高 宽）即第一帧图像的中心点
sigma = 100;


gsize = size(im);   %获取第一帧图片的尺寸:高 宽 3通道
[R,C] = ndgrid(1:gsize(1), 1:gsize(2)); %ndgrid 生成gsize(1)xgsize(2)的高维矩阵，参数1代表纵向/参数二代表横向
%ndgrid用法：https://blog.youkuaiyun.com/u012183487/article/details/76149279
%产生高斯形状的理想响应g
g = gaussC(R,C, sigma, center);  %调用高斯函数，得到高斯图

怎么解释呢？我们得到的第一帧的图片尺寸，先建立起一个二维矩阵，这个二维矩阵的R和C分别代表横坐标和纵坐标。通过左右两个矩阵的值（比如R矩阵取2，C矩阵取10）.则映射到这一帧图像每个坐标点（的“候选样本”的目标中心点坐标）。（相当于R，C提供索引值确定数组对应位置的值一样。）

（4）于是我们就按照R,C的索引，将得到的第一帧图像输入到高斯函数中，而center作为目标中心点，也传入进去。最终得到了当前帧所有坐标点的一个响应值输出g（因为一张图里有很多个坐标点，所有坐标点得到的响应值即组成了一张响应值图，多个目标则有多个峰值）。（顶峰就代表那一个坐标点对应的其实就是目标中心点的响应值输出g）

############################## 响应图的解释：##############################

（5）根据第一帧初始状态的仿射变换与相关计算，，通过最小化该函数，求得初始化后的H。此外，大写字母（F H G）均为傅里叶变化后的对应(f h g)，目的是为了加速运算。

所以我们就得到了当前帧求得的滤波器H。这样也就方便该滤波器和下一帧图片进行计算，最终得出最大响应值。即为下一帧的目标中心点。获得了下一帧的目标中心点后，我们则可按照上述求解滤波器的方法对H进行更新。依次循环（这是对G/F=H这个公式的理解；它的目的是找到一个滤波器可以通过上一帧找到这一帧的目标中心点）

TIPS：

按照这个思想，mosse算法提出的和这个大同小异。也用到了响应图输出G和F。因为Ai-1和Bi-1在第一帧没有数值，所以就设为0，以此来初始化Ai和Bi。

（6）在后续在线更新过程中，则会根据前一帧的响应值输出g将会求出对应的坐标值（响应值即为目标中心点，宽高不变，通过中心点的变化，动态求的bbox坐标），这也将作为下一帧的输入，用于裁减对应区域的对应图像。持续更新滤波器h。在整个跟踪过程中，通过不断更新滤波器h，不断调整局部区域内目标的中心点位置，进而持续调整下一帧裁切的局部区域，从而达到跟踪的目的。这里需要注意的是，mosse相关滤波并不是在全图进行处理，而是基于宽高不变的初始框选区域，在不断迭代更新的该区域进行目标中心点峰值点的一个查找。

Mosse.m主代码：

% get images from source directory
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'),   %num2str将数字转换为字符数组/
    img_files = num2str((1:seq_len)', [img_path '%04i.jpg']);   %得到所有对应的图片
    %disp(img_files);
else
    error('No image files found in the directory.');
end

% select target from first frame
im = imread(img_files(1,:));    %读取第一帧的图片
%[X,map]=imread('forest.png')  则代表X为图像颜色值，map代表色素(通道)
f = figure('Name', 'Select object to track'); imshow(im);   
%figure表示弹出框名称/imshow表示在弹出窗口显示第一帧图片
rect = getrect; %getrect函数即用鼠标指定矩形
%其中rect返回值为选定矩形以此为：左下角的坐标，宽度，高度
close(f); clear f;
center = [rect(2)+rect(4)/2 rect(1)+rect(3)/2];%获得矩形框的中心点纵，横坐标（高 宽）即第一帧图像的中心点

% plot gaussian
sigma = 100;
gsize = size(im);   %获取第一帧图片的尺寸:高 宽 3通道
[R,C] = ndgrid(1:gsize(1), 1:gsize(2)); %ndgrid 生成gsize(1)xgsize(2)的高维矩阵，参数1代表纵向/参数二代表横向
%ndgrid用法：https://blog.youkuaiyun.com/u012183487/article/details/76149279
%产生高斯形状的理想响应g
g = gaussC(R,C, sigma, center);  %调用高斯函数
g = mat2gray(g);%把一个double类的任意数组转换成值范围在[0,1]的归一化double类数组

% randomly warp original image to create training set
if (size(im,3) == 3) %size(A,n)，size将返回矩阵的行数或列数
    img = rgb2gray(im);  %将RGB图像或彩色图转换为灰度图像
end
%imcrop表示按第一帧选择的矩形大小裁剪图片
img = imcrop(img, rect); %按选定的矩形框裁剪第一帧灰度图大小
%imcrop ： https://blog.youkuaiyun.com/llxue0925/article/details/80431508
g = imcrop(g, rect); %
G = fft2(g);   %二维快速傅里叶变换，将响应输出g变成G
height = size(g,1);   %height=179
width = size(g,2);    %width=130
%f表示输入图像
fi = preprocess(imresize(img, [height width])); 
%重新设置裁剪后的图片尺寸，再调用preprocess函数对图像数据进行标准化处理
Ai = (G.*conj(fft2(fi))); %conj函数用于计算"复数"x的共轭值。即F*--》F的共轭
Bi = (fft2(fi).*conj(fft2(fi)));
%按照求和公式计算第一个H
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.125;   %官方推荐的eta值0.125
fig = figure('Name', 'MOSSE');
mkdir(['results_' dataset]);
for i = 1:size(img_files, 1)  %遍历所有图片
    %size(img_files, 1) 表示图片总数量
    img = imread(img_files(i,:));   %读取第i张图片
    im = img;
    if (size(img,3) == 3)  %若通道为3
        img = rgb2gray(img);  %转成灰度图
    end
    if (i == 1)      %当为第一张图时
        Ai = eta.*Ai;
        Bi = eta.*Bi;
    else
        Hi = Ai./Bi;   %不是第一张图的时候，求解滤波H*
        fi = imcrop(img, rect);   %再按指定矩形大小进行裁剪2-376的对应图片
        fi = preprocess(imresize(fi, [height width]));  %预处理得到输出
        gi = uint8(255*mat2gray(ifft2(Hi.*fft2(fi))));  %得到响应值输出 转换成值范围在[0,1]的归一化double类数组 
        maxval = max(gi(:));  %得到最大响应值
        [P, Q] = find(gi == maxval);  %找到最大响应值所对应的gi
        dx = mean(P)-height/2;
        dy = mean(Q)-width/2;

        
        rect = [rect(1)+dy rect(2)+dx width height];  %更新中心点坐标
        fi = imcrop(img, rect);     %按新更新的rect目标点进行裁剪图片 
        fi = preprocess(imresize(fi, [height width]));
        %最后更新Ai和Bi
        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); %绘制矩形 宽为3
    %imwrite(result, ['results_' dataset num2str(i, '/%04i.jpg')]);
    %imwrite(result, ['results_temp' dataset num2str(i, '/%04i.jpg')]);
    imshow(result);
end