KCF代码阅读笔记（matlab版）

因为自己是个憨憨，所以写的比较琐碎，很多没见过的matlab函数也解释了一下，代码是按运行顺序排的。tracker.m中调用了较多子函数，它们的解释在后面。

1. 主函数run_tracker

提取HOG特征用的函数和CPP文件作者已经给出并编译
什么也不输入时，默认为使用HOG特征和高斯核函数（还可以选择使用灰度特征和线性核）
要把base_path改成自己的benchamark文件夹所在路径（第43行），比如：
base_path = ‘F:\track\Benchmark\’;

主函数中给出的相关参数：
padding 目标搜索范围系数，决定了每帧进行检测的范围
Lambada 正则化参数λ
output_sigma_factor 空间带宽：与目标大小成比例
interp_factor 时间更新系数
sigma 高斯核函数的带宽/方差
poly_a/poly_b
cell_size 求HOG特征时cell的尺寸

1-1 选择序列choose_video

输入base_path即benchmark序列库的路径，输出你选择的序列的名称（如’boy’，’Car4’等）
利用listdlg函数生成了一个列表GUI，方便进行选择

1-2 提取序列中目标的参数load_video_info

输入base_path和通过choose_video选择得到的序列名video_name，得到所选序列的存储路径。

video_path = [base_path video '/'];

并且打开benchmark中预留的groundtruth_rect.txt文件，将序列中目标的参数存储到ground_truth中。

filename = [video_path 'groundtruth_rect' suffix '.txt'];
f = fopen(filename);

ground_truth是一个矩阵，行数为序列的帧数。开始时列数为4，第三列和第四列为目标跟踪框的宽和高，第一列和第二列为目标跟踪框左上角的列坐标和行坐标_（这样我猜大概是因为imcrop函数中行和列的位置是反过来的吧）_。跟踪框左上角坐标再加上其二分之一的边长后由4列变为2列，表示目标的中心位置轨迹。
将第一帧的参数作为初始化时的参数target_sz和pos*（有4个序列作者应该是出于跟踪精度的需要修改了参数，暂且不表）*。

target_sz = [ground_truth(1,4), ground_truth(1,3)];%目标大小(高，宽）
pos = [ground_truth(1,2), ground_truth(1,1)] + floor(target_sz/2);%中心点坐标
ground_truth = ground_truth(:,[2,1]) + ground_truth(:,[4,3]) / 2;%4列变为2列

img_files按顺序存储了所有的序列图片的文件名，如’0001.jpg’（图片不必设置为0001.jpg，0002.jpg…0030.jpg这种格式，即使是从一个完整序列中截取的片段，也可以顺利读取，比如0031.jpg,0033.jpg…0080.jpg），长度为图片数。
如果要加入自己采集的序列，最好将名称都改成“xxxx.jpg”的格式，并且要放在一个叫img的子文件夹下，比如这个路径：“F:\track\Benchmark\IRA\img\0001.jpg”。

video_path = [video_path 'img/'];%图片应放在video_path路径下的img子文件夹中
img_files = dir([video_path '*.jpg']);

可以自己仿照benchmark里的groundtruth_rect.txt进行标注，
比如“F:\track\Benchmark\IRA\groundtruth_rect.txt”。不标注的话就需要对代码做一些修改，使其可以顺利读入无标注的序列，并通过imcrop的方式初始化跟踪窗口。

2.跟踪的核心tracker

下面进行逐行详解，其中频域的变量都以f结尾
41——45尺寸调整，如果目标尺寸过大就适当缩小

49，确定检测范围window_sz，作者将其设为跟踪框尺寸的1+padding倍_（如果有快速运动，目标在两帧间的移动超过了这个范围，会导致跟踪丢失？）_

window_sz = floor(target_sz * (1 + padding)); 

58，用output_sigma_factor，cell_sz和跟踪框尺寸计算高斯标签的带宽output_sigma
59，先用gaussian_shaped_labels生成回归标签，然后进行傅里叶变换转换到频域上的yf

output_sigma = sqrt(prod(target_sz)) * output_sigma_factor / cell_size;
yf = fft2(gaussian_shaped_labels(output_sigma, floor(window_sz / cell_size)));

62，生成汉宁窗cos_window，尺寸与yf相同，即floor(window_sz / cell_size)，检测范围window_sz中cell的个数。
对信号进行傅里叶变换时，为了减少频谱泄漏，通常在采样后对信号加窗。

cos_window = hann(size(yf,1)) * hann(size(yf,2))';

65——67，利用show_video将跟踪结果可视化
73，position记录每帧目标中心位置，time记录用时，图片读取和显示的时间不计入

time = 0;  %to calculate FPS
positions = zeros(numel(img_files), 2);  %to calculate precision

75，开始逐帧进行跟踪
76——86，图像预处理，图片读取为im，转为灰度图，图片读取的时间不计入总时间

    for frame = 1:numel(img_files)
        im = imread([video_path img_files{frame}]);
        if size(im,3) > 1
            im = rgb2gray(im);
        end

        tic()

        if frame > 1

87–118行为第一帧后的跟踪过程，跳过来到119行，利用第一帧图像进行初始化
119，由benchamrk中的ground_truth或者自己框选出的第一帧目标位置及之前计算出的检测范围window_sz，利用get_subwindow获得图像中第一帧的检测的区域patch，如果超过图像尺寸会加以修正。作者的修正方法是认为超出部分的值都与边界的值相同。
120，利用get_feature获得第一帧patch的特征矩阵，再经过傅里叶变换到频率域得到xf
125，利用gaussian_correlation得到频率域上的高斯响应kf
131，岭回归计算，得到分类器参数alphaf
133——135，将以上得到的alphaf和xf作为第一帧时分类器训练的结果，初始化结束

 patch = get_subwindow(im, pos, window_sz);

xf = fft2(get_features(patch,cell_size, cos_window));

kf = gaussian_correlation(xf, xf, kernel.sigma);

alphaf = yf ./ (kf + lambda);

if frame == 1  %first frame, train with a single image
    model_alphaf = alphaf;
    model_xf = xf;

①初始化完成后，从第二帧开始正式跟踪
90，在上一帧的跟踪结果pos的基础上，根据window_sz在这一帧图像上提取检测区域patch
这一帧的训练，与下一帧的检测，使用的是同一块patch
91，利用get_feature得到检测区域patch的特征zf
96，由zf和model_xf计算高斯响应kzf
102，model_alphaf与kzf点乘后进行傅里叶反变换，回到时域。保留实部，得到实数响应图response map，并且响应均为归一化值
109——114
115，将响应最大的位置作为目标中心位置，这部分中有一段没有看懂

if frame > 1

    patch = get_subwindow(im, pos, window_sz);

    zf = fft2(get_features(patch , cell_size, cos_window));

    kzf = gaussian_correlation(zf, model_xf, kernel.sigma);

    response = real(ifft2(model_alphaf .* kzf));

    %vert_delta, horiz_delta为response map上峰值最大处
    [vert_delta, horiz_delta] = find(response == max(response(:)), 1);

    %没有理解这一步的意义
    if vert_delta > size(zf,1) / 2  
        vert_delta = vert_delta - size(zf,1);
    end
    if horiz_delta > size(zf,2) / 2 
        horiz_delta = horiz_delta - size(zf,2);
    end

    pos = pos + cell_size * [vert_delta - 1, horiz_delta - 1];

②这一帧的跟踪完成后，用新的目标训练分类器，部分代码与初始化时的训练相同
119，得到这一帧的patch位置
120，得到patch的特征xf
125，计算高斯响应kf
131，岭回归计算，得到新的分类器参数alphaf
133——139，model_alphaf和model_xf随着时间更新，更新系数为主函数中设置的interp_factor

 patch = get_subwindow(im, pos, window_sz);

xf = fft2(get_features(patch,cell_size, cos_window));

kf = gaussian_correlation(xf, xf, kernel.sigma);

alphaf = yf ./ (kf + lambda);

if frame == 1 
    model_alphaf = alphaf;
    model_xf = xf;
else
     model_alphaf = (1 - interp_factor) * model_alphaf + interp_factor * alphaf;
     model_xf = (1 - interp_factor) * model_xf + interp_factor * xf;
 end

①，②检测和训练交替进行，从而实现跟踪
这一帧的训练，与下一帧的检测，使用的是同一块patch
KCF算法采用固定间隔的更新策略，即每一帧都进行分类器的训练，并以固定的系数进行融合

143,144，用position记录跟踪的目标中心位置，time累计时间
147——160，跟踪结果可视化

2-1 高斯形状和回归标签计算gaussian_shaped_labels

输入高斯带宽output_sigma，以及检测范围内cell的个数=floor(window_sz / cell_size)
借助ndgrid函数生成了回归标签labels
ndgrid函数 [X,Y] = ndgrid(x,y)：
输入x和y分别为长度为Lx和Ly的行向量
输出X和Y是维度为Lx×Ly的矩阵，矩阵X的每一列都是向量xT，矩阵Y每一行都是向量y

[rs, cs] = ndgrid((1:sz(1)) - floor(sz(1)/2), (1:sz(2)) - floor(sz(2)/2));
labels = exp(-0.5 / sigma^2 * (rs.^2 + cs.^2));

此时生成的回归标签的峰值在中心，再用circshift将labels的峰值移到左上角

labels = circshift(labels, -floor(sz(1:2) / 2) + 1);
assert(labels(1,1) == 1)%检验峰值是否移动到label的左上角处，即label(1,1)位置

2-2跟踪结果可视化show_video

2-3检测区域提取get_subwindow

输入该帧的图像，目标中心位置及开始计算出的检测范围window_sz，提取图像中的检测区域patch，如果范围超过图像尺寸会加以修正。作者的修正方法是认为超出部分的灰度值都与边界像素的灰度值相同。

%以目标中心位置为中心，扩展出大小等于window_sz的矩形窗口
xs = floor(pos(2)) + (1:sz(2)) - floor(sz(2)/2);
ys = floor(pos(1)) + (1:sz(1)) - floor(sz(1)/2);

%当patch超出图片边缘时，认为超出部分的灰度值与边缘相同
xs(xs < 1) = 1;
ys(ys < 1) = 1;
xs(xs > size(im,2)) = size(im,2);
ys(ys > size(im,1)) = size(im,1);

%在图像中提取出检测区域patch
out = im(ys, xs, :);

2-4特征提取get_features

输入从图像中提取的patch，主函数中已设置好的存在features结构体中的特征类型与参数和HOG的cell尺寸cell_size，以及之前根据回归标签yf生成的汉宁窗 cos_window。
HOG特征利用fhog计算，输出的值是一个三维的矩阵，第三维度的最后值全为0，因此置空，第三维的大小变为31。之后与汉宁窗cos_window相乘。
最终输出的特征矩阵x，前两维尺寸与回归标签yf、汉宁窗尺寸相同，均为window_sz / cell_size，第三维为31。
灰度特征就是图像像素灰度归一化之后减去均值
可以在此处加入自己设计的特征，比如深度卷积特征，HCF算法就是如此

x = double(fhog(single(im) / 255, cell_size));
x(:,:,end) = [];%将第三维度的最后置空

%为patch的特征x加窗，做傅里叶变换的准备
if ~isempty(cos_window)
    x = bsxfun(@times, x, cos_window);
end

2-5HOG特征计算fhog

输入归一化的图像patch，HOG的cell_size和主函数里选择的features.hog_orientations即特征参数，利用了一个工具箱里的函数gradientMex，输出一个三维的矩阵，数据类型为单精度，前两维的尺寸等于patch尺寸的四分之一，与汉宁窗cos_window的尺寸相同，第三维大小为32。

2-6响应度计算gaussian_correlation

输入两个频率域上的尺寸完全相同的两个三维矩阵xf，yf以及主函数中选择的高斯方差gaussian.sigam。 N为前两维尺寸之积，即“面积”。
xx和yy等于xf和yf中各元素的模的平方和再除以N，都是实数**（Matlab的转置对复矩阵自带共轭效果！！！！！）**
xyf等于xf与yf的共轭矩阵点乘
xyf经过傅里叶反变换回到时域，取实部并将第三维度上的值求和，最后得到尺寸与输入xf和yf相同的二维矩阵xy
用xx+yy-2xy再除以xf中的元素数，然后乘以高斯核函数，最后经过傅里叶变换得到所有位置的高斯响应_（我觉得N和numel(xf)其实是相等的？）_

N = size(xf,1) * size(xf,2);
xx = xf(:)' * xf(:) / N;  %squared norm of x
yy = yf(:)' * yf(:) / N;  %squared norm of y

%cross-correlation term in Fourier domain
xyf = xf .* conj(yf);
xy = sum(real(ifft2(xyf)), 3);  %to spatial domain

%calculate gaussian response for all positions, then go back to the
%Fourier domain
kf = fft2(exp(-1 / sigma^2 * max(0, (xx + yy - 2 * xy) / numel(xf))));

本文转自 https://blog.csdn.net/bitopyx/article/details/81948875，如有侵权，请联系删除。