本文是CSK与KCF算法推导的第五篇，主要介绍具体实现的细节。

算法流程

  算法的大致流程如上图所示。左边是第一次训练，右边是每次的更新。我认为整个算法最关键的地方在于提取特征图和核函数的计算。作者的代码里有提取hog特征，也有直接用原始的灰度图像数据作为特征图，目前我只看了用灰度图作为特征图的部分。需要注意的是，这里的灰度图是归一化到[-0.5,0.5]的区间内的！！！

fftshift()

  MATLAB里有一个fftshift()的函数。原本DFT的结果中，低频部分是在序列两端，高频在中间；这个函数可以把低频移到中间，高频移到两边。

  在我前面的文章中介绍了，标签值应该是一个峰值在序列两端的高斯函数。既然我们是用这样的标签训练的，我们在检测的时候，如果检测一个和训练样本类似的样本，那么得到的结果也会是一个峰值在两端的结果。这时候可以用类似fftshift()的函数将结果重新排布一下，这样能便于分析。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// Create Gaussian Peak. Function called only in the first frame.
cv::Mat KCFTracker::createGaussianPeak(int sizey, int sizex)
{
    cv::Mat_<float> res(sizey, sizex);

    int syh = (sizey) / 2;
    int sxh = (sizex) / 2;

    float output_sigma = std::sqrt((float) sizex * sizey) / padding * output_sigma_factor;
    float mult = -0.5 / (output_sigma * output_sigma);

    for (int i = 0; i < sizey; i++)
        for (int j = 0; j < sizex; j++)
        {
            int ih = i - syh;
            int jh = j - sxh;
            res(i, j) = std::exp(mult * (float) (ih * ih + jh * jh));
        }
    return FFTTools::fftd(res);
}

  上面是作者采用的生成标签值的代码，可以看到这个生成得到标签值并没有符合峰值在四个角点的要求。但作者的代码确实是能用的，为什么呢？因为作者是在计算高斯核相关的时候做了fftshift。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
// Evaluates a Gaussian kernel with bandwidth SIGMA for all relative shifts between input images X and Y, which must both be MxN. They must    also be periodic (ie., pre-processed with a cosine window).
cv::Mat KCFTracker::gaussianCorrelation(cv::Mat x1, cv::Mat x2)
{
    using namespace FFTTools;
    cv::Mat c = cv::Mat( cv::Size(size_patch[1], size_patch[0]), CV_32F, cv::Scalar(0) );
    // HOG features
    if (_hogfeatures) {
……
    }
    // Gray features
    else {
        cv::mulSpectrums(fftd(x1), fftd(x2), c, 0, true);
        c = fftd(c, true);
        rearrange(c);
        c = real(c);
    }
……

  上面的rearrange()函数就是完成了fftshift()的功能。为什么能够这么做呢？

α^=y^k^xx+λ\hat \alpha = \frac{ {\hat y} }{ { { {\hat k}^{xx} } + \lambda } }

  我们可以对作者代码中采用的计算方式做一个推导。记α′\alpha'是代码中计算的结果，yy是一个峰值在四周的二维高斯函数，y′y'是峰值在中央的二维高斯函数，矩阵的维度是M×N。那么：

α^′=F(y′)F(kxx(x−N/2,  y−M/2))+λ=F(y(x−N/2,  y−M/2))F(kxx(x−N/2,  y−M/2))+λ\hat \alpha ' = \frac{ { {\mathcal F}(y')} }{ { {\mathcal F}\left( {k^{xx}(x - N/2,\;y - M/2)} \right){\rm{ + } }\lambda } }= \frac{ { {\mathcal F}(y(x - N/2,\;y - M/2))} }{ { {\mathcal F}\left( { {k^{xx} }(x - N/2,\;y - M/2)} \right){\rm{ + } }\lambda } }

  二维fftshift()相当于是在矩阵行列方向上分别循环移位一半的长度。自然yy和y′y'也是fftshift()的关系。利用傅立叶变换的时移性质可以得出：

α^′(u,v)=y^WNuN/2WMvM/2k^xxWNuN/2WMvM/2+λ=y^(−1)u+vk^xx(−1)u+v+λ\hat \alpha '(u,v) = \frac{ {\hat yW_N^{uN/2}W_M^{vM/2} } }{ { { {\hat k}^{xx} }W_N^{uN/2}W_M^{vM/2} + \lambda } } = \frac{ {\hat y{ {( - 1)}^{u + v} } } }{ { { {\hat k}^{xx} }{ {( - 1)}^{u + v} } + \lambda } }

=y^(−1)u+vk^xx(−1)u+v+λ(−1)u+v=y^k^xx+λ=α^(u,v)= \frac{ {\hat y{ {( - 1)}^{u + v} } } }{ { { {\hat k}^{xx} }{ {( - 1)}^{u + v} } + \lambda { {( - 1)}^{u + v} } } } = \frac{ {\hat y} }{ { { {\hat k}^{xx} } + \lambda } } = \hat \alpha(u,v)

  所以作者代码中的做法是可行的。

峰值点估计

  因为处理的图像的点都是离散的，所以我们在f(z)f(z)中找到的某一个峰值点，不一定恰好就是真正的峰值，需要对真正的峰值点进行估计。

  用检测到的峰值点和与它临近的两个点可以对真正的峰值点进行估计。用三个点来拟合二次函数，非常简单，设f(x)=ax2+bx+qf(x)=ax^2+bx+q，将三个点代入可以得到：

{a−b+q=lq=ca+b+q=r\left\{ \begin{array}{l} a - b + q = l\\ q = c\\ a + b + q = r \end{array} \right.

x=−b2a=−r−l2(l+r−2c)=r−l2(2c−l−r)x = - \frac{b}{ {2a} } = - \frac{ {r - l} }{ {2(l + r - 2c)} } = \frac{ {r - l} }{ {2(2c - l - r)} }

  下面是作者采用的代码：

1
2
3
4
5
6
7
8
9
10
// Calculate sub-pixel peak for one dimension
float KCFTracker::subPixelPeak(float left, float center, float right)
{   
    float divisor = 2 * center - right - left;

    if (divisor == 0)
        return 0;
    
    return 0.5 * (right - left) / divisor;
}

峰值位置与目标位置

  我们按照这样的样本C(x)C(x)和标签yy训练得到了α\alpha，那么如果接下来给一个下面这样的样本zz：

  那么会得到怎样的响应呢？先把zz写成循环矩阵，这一循环矩阵中的样本分别求出响应值，就得到了f(z)f(z)，峰值点相对于标签值左移了一个单位，而样本则是右移一个单位。

  所以f(z)f(z)中的峰值点相对于原点的偏移表示当前的样本需要进行这样的移位才能和原来的训练样本一样；也就是说当前的样本是朝反方向偏移了。对于二维的图像，f(z)f(z)的峰值点相对于原点偏移了(Δx, Δy)(\Delta x,\ \Delta y)，那么目标实际上是从原来的位置偏移了(−Δx, −Δy)(-\Delta x,\ -\Delta y)。
  那么怎样才能让结果符合我们的直觉呢？只要在f(z)f(z)的行列索引前面加负号，进行倒序就行。根据傅立叶变换的性质，x[−n]=X∗[k]x[-n]=X^*[k]，所以我们其实可以直接求f^∗(z)\hat f^*(z)，而不是求f^(z)\hat f(z)。
  再结合我的上一篇文章中的内容：

f^∗(z)=k^zx⊙α^∗=k^zx⊙real(α^){ {\hat f}^*}(z) = { {\hat k}^{zx} } \odot { {\hat \alpha }^*} = { {\hat k}^{zx} } \odot real(\hat \alpha )

  其中k^zx{\hat k}^{zx}表示在求zz和xx的互相关时，需要对x^\hat x取复共轭，这也是作者的代码中实际采用的做法。