KCFコードリーディング - takminの書きっぱなし備忘録 @はてなブログ

はじめに

この記事はOpenCV Advent Calender 2016の第22日目の記事です。

OpenCVのextraモジュールの中にTracking APIというものがあり、それを使用すると動画中の物体を追跡するための様々なアルゴリズムを使用することが出来ます。

Tracking APIに使い方については、以下の記事を参考にしてください。

OpenCV Tracking API について

こんなに簡単!? トラッキング

ここではその中のKernelized Correlation Filter (KCF)というアルゴリズムについて、コードリーディングしてみたので、それについて解説します。

今回、コードリーディングをしようと思った理由は２つです。

Trackingの使い方の記事は見かけるけど、KCFのアルゴリズムについて解説した記事は見当たらない。論文自体が結構面白いのでぜひ紹介したい。
Trackingはextraモジュールに含まれているが、extraモジュールはバージョンが変わると突然使えなくなったりするので、そうなった時のためにも自分でコードの中身まで理解しておきたい。

以下では、まずアルゴリズムについて解説し、次に実際のコードを見ていきます。

アルゴリズムの解説

KCFは物体を追跡しながら物体テンプレートを随時学習する手法です[1]。

この動画のように、第１フレームで指定した物体を追跡し続けます。
この時、学習には１フレームだけあれば良く、後は随時物体を追跡しながら学習を更新していきます。

物体テンプレートを随時学習しながら追跡を行う他の手法として、代表的なものにTLD（またはP-N Learning）という手法があります[2]。TLDでは追跡している物体から正例を、その物体周辺から負例画像を抽出し、識別器（randomized fern）を都度学習することで追跡を行っています。（TLDについてはこちらの記事を参照。尚Tracking APIにも含まれています。）

KCFでは物体検出を、入力画像全体から物体位置を抽出する回帰問題として捉えます。その際画像をシフトさせることで、１枚の画像から学習画像を仮想的に増やして学習を行います。

KCFの面白い点は、シフトさせた画像を巡回行列ととらえ、この回帰問題をリッジ回帰の問題にすることで、最終的に周波数領域におけるスペクトル同士の乗算/除算の形で簡単に解けることを示したことです。これは結果としてCorrelation Filterと呼ばれるものと同じ手法になりますが、KCFでは更にカーネルを適用した場合でも同様に周波数空間でのスペクトル同士の乗算/除算に落とし込めることを示しました。

コードを眺める前に、もう少しアルゴリズムについて論文[1]を元に解説します。尚、式の導出については論文を参照して下さい。

リッジ回帰は、以下の式で表せます。

ただし、

です。

ここで、 $\bf{x}_i$ は $i$ ピクセルずつシフトした学習画像、 $y_i$ はそれに対応する物体の位置（教師ラベル）、 $\lambda$ は過学習を防ぐための正則化項への重みです。ここでの目的は入力画像にかけると物体位置を返すようなフィルター $\bf{w}$ を求めることです。

この問題は解析的に解くことが可能で、フィルター $\bf{w}$ は以下の式で求まります。（式番号は論文[1]と合わせています）

(3)式は都合により複素空間にまで解を広げています。ここで、 $\bf{X}$ は少しずつシフトした画像を行ごとに並べた行列で、結果として(6)式のような巡回行列になります。

ここで、 $\bf{X}^H$ は $\bf{X}$ のエルミート転置と呼ばれるもので、 $\bf{X}$ の複素共役の転置を取ったものです。

ここで、 $\bf{X}$ が巡回行列であることを利用すると、 $\bf{w}$ は以下のように導出できます。

ここで $\hat{\bf{w}}$ 、 $\hat{\bf{x}}$ 、 $\hat{\bf{x}}^*$ 、 $\hat{\bf{y}}$ は、それぞれ $\bf{w}$ 、 $\bf{x}$ 、 $\bf{x}^*$ ( $\bf{x}$ の複素共役)、 $\bf{y}$ の離散フーリエ変換により求めた周波数スペクトルです。 $\odot$ はベクトルの要素同士の乗算を表し、またここでは分母と分子を要素同士で割り算しています。

このようにフィルター $\bf{w}$ を周波数スペクトル同士の乗算/除算の形で表すことができました。
あとはこのフィルターで画像内を畳み込み、最も高い値を返す場所を物体の個所とみなすことができます。畳み込み演算は周波数領域ではやはりスペクトル同士の乗算の形で表されるので、ほぼ全ての計算が周波数領域で完結します。

さらに入力 $\bf{X}$ にカーネルを適用した場合を考えます。

カーネル法とは、入力 $\bf{x}$ をなんらかの非線形な関数 $\varphi(\bf{x})$ によって特徴を抽出することで識別性能を高めるようなアプローチです。この際、非線形変換したベクトル同士の内積を $\varphi^T(\bf{x})\varphi(\bf{x}')=\kappa(\bf{x},\bf{x}')$ とおくことで、色々な計算を簡略化できます。この $\kappa(\bf{x},\bf{x}')$ のことをカーネル関数と呼びます。

入力画像およびそれをシフトさせた画像に対し、非線形変換を行った $\varphi(\bf{x}_i)$ を新たな入力とした場合、やはり周波数領域でほとんどの計算を行うことが出来、以下のように計算できます。

ここで、 $\bf{k^{xx'}}$ はベクトル $\varphi(\bf{x'})$ と $\varphi(\bf{x}_i)$ との内積を要素に持つベクトル（この論文ではカーネル相関と呼んでいる）で、その $i$ 番目の要素 $k^{xx'}_i$ は(18)式で表せます。

$P^{i-1}\bf{x}$ は入力画像 $\bf{x}$ をi-1回シフトさせた画像です。

このカーネル相関ベクトル $\bf{k^{xx'}}$ を離散フーリエ変換した $\hat{\bf{k}}^{\bf{xx'}}$ を用いることで、周波数領域での乗算/除算によって未知の画像に対する物体位置の応答を求めることができました。

ただし、(17)式が成り立つためには巡回行列 $\bf{X}$ に対するグラム行列 $\bf{K}$ が、やはり巡回行列とならなければいけないという制約があります。このようなカーネル関数としては、RBFカーネルや多項式カーネル、ガウスカーネル、χ2乗カーネルなどがあります。

尚、OpenCVの実装ではガウスカーネルを使用しています。その場合 $k_i^{xx'}$ は(30)式で表せます。

ここで $\it{F}^{-1}$ は逆フーリエ変換を表します。

また、(30)式はカラー画像やあるいはHOG特徴量などのマルチチャネルへと拡張できます。

ここでcがチャネルを表します。実際にOpenCVの実装ではカラー画像や特徴量などのマルチチャネルを扱うことができます。

コードの解説

前置きが長くなりましたが、実際にコードを覗いていきます。
（とはいえ、前置き部分ですでに力が尽きかけているので、さっと流していきます。。。）

尚、バージョンは3.1のリリース版を使用しています。

https://github.com/opencv/opencv_contrib/releases

このうち、KCFの実装部分はほとんどが

modules/tracking/src/trackerKCF.cpp

というファイルに記述されています。

KCFはcv::Trackerクラスによってラップされており、基本的にこのクラス経由で使用します。

使い方についてはこちらのページが参考になります（サンプルコードあり）

http://docs.opencv.org/3.1.0/d2/d0a/tutorial_introduction_to_tracker.html

KCFトラッカーはまず、

Ptr<Tracker> tracker = Tracker::create( "KCF" );

とすることで、KCFトラッカーを生成し、次にトラッカーの初期化を行います。

tracker->init(frame,roi);

ここで、frameは入力画像(cv::Mat)、roiは追跡したい物体の位置(cv::Rect)です。

以後は、カメラ/動画からのフレームを１枚１枚取得して、テンプレートを更新しながら追跡をしていきます。

 tracker->update(frame,roi);

ここでroiは入力として前のフレームにおける物体位置を与え、出力として入力画像における物体位置を返します。

Tracker::init()とTracker::update()関数内では、TrackerKCFImpl::initImpl()とTrackerKCFImpl::updateImpl()関数がそれぞれ呼び出されるため、ここでは主にこの２つの関数について順を追ってみていきます。

TrackerKCFImpl::initImpl()

  /*
   * Initialization:
   * - creating hann window filter
   * - ROI padding
   * - creating a gaussian response for the training ground-truth
   * - perform FFT to the gaussian response
   */
  bool TrackerKCFImpl::initImpl( const Mat& /*image*/, const Rect2d& boundingBox ){
    frame=0;
    roi = boundingBox;

    //calclulate output sigma
    output_sigma=sqrt(roi.width*roi.height)*params.output_sigma_factor;
    output_sigma=-0.5/(output_sigma*output_sigma);

    //resize the ROI whenever needed
    if(params.resize && roi.width*roi.height>params.max_patch_size){
      resizeImage=true;
      roi.x/=2.0;
      roi.y/=2.0;
      roi.width/=2.0;
      roi.height/=2.0;
    }

    // add padding to the roi
    roi.x-=roi.width/2;
    roi.y-=roi.height/2;
    roi.width*=2;
    roi.height*=2;

    // initialize the hann window filter
    createHanningWindow(hann, roi.size(), CV_64F);

    // hann window filter for CN feature
    Mat _layer[] = {hann, hann, hann, hann, hann, hann, hann, hann, hann, hann};
    merge(_layer, 10, hann_cn);

    // create gaussian response
    y=Mat::zeros((int)roi.height,(int)roi.width,CV_64F);
    for(unsigned i=0;i<roi.height;i++){
      for(unsigned j=0;j<roi.width;j++){
        y.at<double>(i,j)=(i-roi.height/2+1)*(i-roi.height/2+1)+(j-roi.width/2+1)*(j-roi.width/2+1);
      }
    }

    y*=(double)output_sigma;
    cv::exp(y,y);

    // perform fourier transfor to the gaussian response
    fft2(y,yf);

    model=Ptr<TrackerKCFModel>(new TrackerKCFModel(params));

    // record the non-compressed descriptors
    if((params.desc_npca & GRAY) == GRAY)descriptors_npca.push_back(GRAY);
    if((params.desc_npca & CN) == CN)descriptors_npca.push_back(CN);
    if(use_custom_extractor_npca)descriptors_npca.push_back(CUSTOM);
    features_npca.resize(descriptors_npca.size());

    // record the compressed descriptors
    if((params.desc_pca & GRAY) == GRAY)descriptors_pca.push_back(GRAY);
    if((params.desc_pca & CN) == CN)descriptors_pca.push_back(CN);
    if(use_custom_extractor_pca)descriptors_pca.push_back(CUSTOM);
    features_pca.resize(descriptors_pca.size());

    // accept only the available descriptor modes
    CV_Assert(
      (params.desc_pca & GRAY) == GRAY
      || (params.desc_npca & GRAY) == GRAY
      || (params.desc_pca & CN) == CN
      || (params.desc_npca & CN) == CN
      || use_custom_extractor_pca
      || use_custom_extractor_npca
    );

    // TODO: return true only if roi inside the image
    return true;
  }

initImple()では、入力として与えられたboundingBoxから(1)式の教師ラベル $\bf{y}$ および周波数スペクトル $\hat{\bf{y}}$ を生成します。また、入力画像に関してはここでは特に使用しません。
テンプレートの学習はTracker::update()を呼び出したときに行われます。

教師ラベルは、下図のようにROIの中心位置が１の値を持つガウス分布で表します。

    //calclulate output sigma
    output_sigma=sqrt(roi.width*roi.height)*params.output_sigma_factor;
    output_sigma=-0.5/(output_sigma*output_sigma);

ここでは、ガウス分布のσを計算しています。

    // add padding to the roi
    roi.x-=roi.width/2;
    roi.y-=roi.height/2;
    roi.width*=2;
    roi.height*=2;

アルゴリズム解説では学習データとして入力画像全体を使っていましたが、実際の実装ではROIのサイズを2倍に拡張した領域についてのみ処理しています。
ここでは、そのために処理領域を２倍に拡張しています。（そのため、1フレームでROIのサイズ以上物体が動くと追跡が失敗するかもしれません）

    // initialize the hann window filter
    createHanningWindow(hann, roi.size(), CV_64F);

    // hann window filter for CN feature
    Mat _layer[] = {hann, hann, hann, hann, hann, hann, hann, hann, hann, hann};
    merge(_layer, 10, hann_cn);

ハニング窓を生成します。
フーリエ変換は本来周期的な信号を入力とした変換のため、入力の境界部分の両端が同じ値を持つ必要があります。そこで入力画像に対してハニング窓をかけることで両端の値が連続するように変換します。

尚、カラー画像に対しては３チャネルから１０チャネルへ変換して処理が行われるため、カラー画像用のハニング窓を１０個用意しています。

    // create gaussian response
    y=Mat::zeros((int)roi.height,(int)roi.width,CV_64F);
    for(unsigned i=0;i<roi.height;i++){
      for(unsigned j=0;j<roi.width;j++){
        y.at<double>(i,j)=(i-roi.height/2+1)*(i-roi.height/2+1)+(j-roi.width/2+1)*(j-roi.width/2+1);
      }
    }

    y*=(double)output_sigma;
    cv::exp(y,y)

応答となるガウス分布を計算しています。ROIの中心がもっとも値が大きくなります。

    // perform fourier transfor to the gaussian response
    fft2(y,yf);

教師ラベル $\bf{y}$ を高速フーリエ変換によって周波数スペクトル $\hat{\bf{y}}$ を求めています。

    // record the non-compressed descriptors
    if((params.desc_npca & GRAY) == GRAY)descriptors_npca.push_back(GRAY);
    if((params.desc_npca & CN) == CN)descriptors_npca.push_back(CN);
    if(use_custom_extractor_npca)descriptors_npca.push_back(CUSTOM);
    features_npca.resize(descriptors_npca.size());

    // record the compressed descriptors
    if((params.desc_pca & GRAY) == GRAY)descriptors_pca.push_back(GRAY);
    if((params.desc_pca & CN) == CN)descriptors_pca.push_back(CN);
    if(use_custom_extractor_pca)descriptors_pca.push_back(CUSTOM);
    features_pca.resize(descriptors_pca.size());

ここでは入力のタイプを設定しています。
入力にはグレースケール（GRAY）、カラー画像（CN)、カスタム(CUSTOM)があり、またGRAYとCNに関しては、それぞれPCAによる次元圧縮を行うか行わないかを選択できます。

TrackerKCFImpl::updateImpl()

  /*
   * Main part of the KCF algorithm
   */
  bool TrackerKCFImpl::updateImpl( const Mat& image, Rect2d& boundingBox ){
    double minVal, maxVal;	// min-max response
    Point minLoc,maxLoc;	// min-max location

    Mat img=image.clone();
    // check the channels of the input image, grayscale is preferred
    CV_Assert(img.channels() == 1 || img.channels() == 3);

    // resize the image whenever needed
    if(resizeImage)resize(img,img,Size(img.cols/2,img.rows/2));

    // detection part
    if(frame>0){

      // extract and pre-process the patch
      // get non compressed descriptors
      for(unsigned i=0;i<descriptors_npca.size()-extractor_npca.size();i++){
        if(!getSubWindow(img,roi, features_npca[i], img_Patch, descriptors_npca[i]))return false;
      }
      //get non-compressed custom descriptors
      for(unsigned i=0,j=(unsigned)(descriptors_npca.size()-extractor_npca.size());i<extractor_npca.size();i++,j++){
        if(!getSubWindow(img,roi, features_npca[j], extractor_npca[i]))return false;
      }
      if(features_npca.size()>0)merge(features_npca,X[1]);

      // get compressed descriptors
      for(unsigned i=0;i<descriptors_pca.size()-extractor_pca.size();i++){
        if(!getSubWindow(img,roi, features_pca[i], img_Patch, descriptors_pca[i]))return false;
      }
      //get compressed custom descriptors
      for(unsigned i=0,j=(unsigned)(descriptors_pca.size()-extractor_pca.size());i<extractor_pca.size();i++,j++){
        if(!getSubWindow(img,roi, features_pca[j], extractor_pca[i]))return false;
      }
      if(features_pca.size()>0)merge(features_pca,X[0]);

      //compress the features and the KRSL model
      if(params.desc_pca !=0){
        compress(proj_mtx,X[0],X[0],data_temp,compress_data);
        compress(proj_mtx,Z[0],Zc[0],data_temp,compress_data);
      }

      // copy the compressed KRLS model
      Zc[1] = Z[1];

      // merge all features
      if(features_npca.size()==0){
        x = X[0];
        z = Zc[0];
      }else if(features_pca.size()==0){
        x = X[1];
        z = Z[1];
      }else{
        merge(X,2,x);
        merge(Zc,2,z);
      }

      //compute the gaussian kernel
      denseGaussKernel(params.sigma,x,z,k,layers,vxf,vyf,vxyf,xy_data,xyf_data);

      // compute the fourier transform of the kernel
      fft2(k,kf);
      if(frame==1)spec2=Mat_<Vec2d >(kf.rows, kf.cols);

      // calculate filter response
      if(params.split_coeff)
        calcResponse(alphaf,alphaf_den,kf,response, spec, spec2);
      else
        calcResponse(alphaf,kf,response, spec);

      // extract the maximum response
      minMaxLoc( response, &minVal, &maxVal, &minLoc, &maxLoc );
      roi.x+=(maxLoc.x-roi.width/2+1);
      roi.y+=(maxLoc.y-roi.height/2+1);

      // update the bounding box
      boundingBox.x=(resizeImage?roi.x*2:roi.x)+boundingBox.width/2;
      boundingBox.y=(resizeImage?roi.y*2:roi.y)+boundingBox.height/2;
    }

    // extract the patch for learning purpose
    // get non compressed descriptors
    for(unsigned i=0;i<descriptors_npca.size()-extractor_npca.size();i++){
      if(!getSubWindow(img,roi, features_npca[i], img_Patch, descriptors_npca[i]))return false;
    }
    //get non-compressed custom descriptors
    for(unsigned i=0,j=(unsigned)(descriptors_npca.size()-extractor_npca.size());i<extractor_npca.size();i++,j++){
      if(!getSubWindow(img,roi, features_npca[j], extractor_npca[i]))return false;
    }
    if(features_npca.size()>0)merge(features_npca,X[1]);

    // get compressed descriptors
    for(unsigned i=0;i<descriptors_pca.size()-extractor_pca.size();i++){
      if(!getSubWindow(img,roi, features_pca[i], img_Patch, descriptors_pca[i]))return false;
    }
    //get compressed custom descriptors
    for(unsigned i=0,j=(unsigned)(descriptors_pca.size()-extractor_pca.size());i<extractor_pca.size();i++,j++){
      if(!getSubWindow(img,roi, features_pca[j], extractor_pca[i]))return false;
    }
    if(features_pca.size()>0)merge(features_pca,X[0]);

    //update the training data
    if(frame==0){
      Z[0] = X[0].clone();
      Z[1] = X[1].clone();
    }else{
      Z[0]=(1.0-params.interp_factor)*Z[0]+params.interp_factor*X[0];
      Z[1]=(1.0-params.interp_factor)*Z[1]+params.interp_factor*X[1];
    }

    if(params.desc_pca !=0 || use_custom_extractor_pca){
      // initialize the vector of Mat variables
      if(frame==0){
        layers_pca_data.resize(Z[0].channels());
        average_data.resize(Z[0].channels());
      }

      // feature compression
      updateProjectionMatrix(Z[0],old_cov_mtx,proj_mtx,params.pca_learning_rate,params.compressed_size,layers_pca_data,average_data,data_pca, new_covar,w_data,u_data,vt_data);
      compress(proj_mtx,X[0],X[0],data_temp,compress_data);
    }

    // merge all features
    if(features_npca.size()==0)
      x = X[0];
    else if(features_pca.size()==0)
      x = X[1];
    else
      merge(X,2,x);

    // initialize some required Mat variables
    if(frame==0){
      layers.resize(x.channels());
      vxf.resize(x.channels());
      vyf.resize(x.channels());
      vxyf.resize(vyf.size());
      new_alphaf=Mat_<Vec2d >(yf.rows, yf.cols);
    }

    // Kernel Regularized Least-Squares, calculate alphas
    denseGaussKernel(params.sigma,x,x,k,layers,vxf,vyf,vxyf,xy_data,xyf_data);

    // compute the fourier transform of the kernel and add a small value
    fft2(k,kf);
    kf_lambda=kf+params.lambda;

    double den;
    if(params.split_coeff){
      mulSpectrums(yf,kf,new_alphaf,0);
      mulSpectrums(kf,kf_lambda,new_alphaf_den,0);
    }else{
      for(int i=0;i<yf.rows;i++){
        for(int j=0;j<yf.cols;j++){
          den = 1.0/(kf_lambda.at<Vec2d>(i,j)[0]*kf_lambda.at<Vec2d>(i,j)[0]+kf_lambda.at<Vec2d>(i,j)[1]*kf_lambda.at<Vec2d>(i,j)[1]);

          new_alphaf.at<Vec2d>(i,j)[0]=
          (yf.at<Vec2d>(i,j)[0]*kf_lambda.at<Vec2d>(i,j)[0]+yf.at<Vec2d>(i,j)[1]*kf_lambda.at<Vec2d>(i,j)[1])*den;
          new_alphaf.at<Vec2d>(i,j)[1]=
          (yf.at<Vec2d>(i,j)[1]*kf_lambda.at<Vec2d>(i,j)[0]-yf.at<Vec2d>(i,j)[0]*kf_lambda.at<Vec2d>(i,j)[1])*den;
        }
      }
    }

    // update the RLS model
    if(frame==0){
      alphaf=new_alphaf.clone();
      if(params.split_coeff)alphaf_den=new_alphaf_den.clone();
    }else{
      alphaf=(1.0-params.interp_factor)*alphaf+params.interp_factor*new_alphaf;
      if(params.split_coeff)alphaf_den=(1.0-params.interp_factor)*alphaf_den+params.interp_factor*new_alphaf_den;
    }

    frame++;
    return true;
  }

updateImpl()では、物体テンプレートの学習と追跡を両方行います。

    // detection part
    if(frame>0){
     .
     .
     .
    }

では、物体の検出を行います。そのため、init()直後に呼び出した場合は、（まだテンプレートが学習されていないため）スキップされます。
この検出の手順を具体的に見てみます。

      // extract and pre-process the patch
      // get non compressed descriptors
      for(unsigned i=0;i<descriptors_npca.size()-extractor_npca.size();i++){
        if(!getSubWindow(img,roi, features_npca[i], img_Patch, descriptors_npca[i]))return false;
      }
      //get non-compressed custom descriptors
      for(unsigned i=0,j=(unsigned)(descriptors_npca.size()-extractor_npca.size());i<extractor_npca.size();i++,j++){
        if(!getSubWindow(img,roi, features_npca[j], extractor_npca[i]))return false;
      }
      if(features_npca.size()>0)merge(features_npca,X[1]);

      // get compressed descriptors
      for(unsigned i=0;i<descriptors_pca.size()-extractor_pca.size();i++){
        if(!getSubWindow(img,roi, features_pca[i], img_Patch, descriptors_pca[i]))return false;
      }
      //get compressed custom descriptors
      for(unsigned i=0,j=(unsigned)(descriptors_pca.size()-extractor_pca.size());i<extractor_pca.size();i++,j++){
        if(!getSubWindow(img,roi, features_pca[j], extractor_pca[i]))return false;
      }
      if(features_pca.size()>0)merge(features_pca,X[0]);

ここでは各入力の処理タイプ（GRAY、CN、CUSTOM）によって入力画像の変換を行い、最終的にPCAによる圧縮を行う信号をX[0]へ、行わないものをX[1]へ格納しています。

getSubWindow()関数は入力画像から探索領域（init()で設定したROIの２倍）を切り取って正規化を行う関数です。その際、入力の処理タイプによって異なる正規化を行います。例えばGRAYの場合は-0.5〜0.5の範囲に正規化してハミング窓をかけます。CNの場合は3チャネルをLook Up Tableを用いて10チャネルに変換した上でハミング窓をかけています。

この10チャネルへの変換は文献[3]に従っています。

      //compress the features and the KRSL model
      if(params.desc_pca !=0){
        compress(proj_mtx,X[0],X[0],data_temp,compress_data);
        compress(proj_mtx,Z[0],Zc[0],data_temp,compress_data);
      }

入力をPCAを用いて次元削減を行っています。

      //compute the gaussian kernel
      denseGaussKernel(params.sigma,x,z,k,layers,vxf,vyf,vxyf,xy_data,xyf_data);

入力xとテンプレートzを用いて(31)式を計算し、値をkに格納します。これは、(22)式の $\bf{k^{xz}$ にあたります。

      // compute the fourier transform of the kernel
      fft2(k,kf);
      if(frame==1)spec2=Mat_<Vec2d >(kf.rows, kf.cols);

$\bf{k^{xz}$ から $\hat{\bf{k}}^{\bf{xz}}$ を求めています。

      // calculate filter response
      if(params.split_coeff)
        calcResponse(alphaf,alphaf_den,kf,response, spec, spec2);
      else
        calcResponse(alphaf,kf,response, spec);

      // extract the maximum response
      minMaxLoc( response, &minVal, &maxVal, &minLoc, &maxLoc );
      roi.x+=(maxLoc.x-roi.width/2+1);
      roi.y+=(maxLoc.y-roi.height/2+1);

calcResponse()関数は(23)式の周波数空間における応答を計算し、それを逆フーリエ変換します。
従って、minMaxLoc()関数で最も反応の強かった位置を検出すれば、それが物体の位置になります。

以降のコードでは学習を行います。まず検出の時と同様、getSubWindow()関数による探索領域の切り出しやハミング窓の適用を行います。この時、もし検出をこの入力画像に対してすでに行われていた場合、検出結果をROIとして新たに設定します。
また、検出時と同様にX[0]とX[1]にそれぞれ次元圧縮用ベクトルと非次元圧縮用ベクトルを格納します。

    //update the training data
    if(frame==0){
      Z[0] = X[0].clone();
      Z[1] = X[1].clone();
    }else{
      Z[0]=(1.0-params.interp_factor)*Z[0]+params.interp_factor*X[0];
      Z[1]=(1.0-params.interp_factor)*Z[1]+params.interp_factor*X[1];
    }

もし一番最初のフレームの場合、テンプレート(Z)に入力画像(X)のROI領域をコピーします。２フレーム以降の場合は、検出した物体の画像領域と過去のテンプレート情報の重み付き平均を用いてテンプレートを更新します。

    // Kernel Regularized Least-Squares, calculate alphas
    denseGaussKernel(params.sigma,x,x,k,layers,vxf,vyf,vxyf,xy_data,xyf_data);

    // compute the fourier transform of the kernel and add a small value
    fft2(k,kf);
    kf_lambda=kf+params.lambda;

(31)式を用いて $\hat{\bf{k}}^{\bf{xx}}$ および $\hat{\bf{k}}^{\bf{xx}}+\lambda$ を計算します。

    double den;
    if(params.split_coeff){
      mulSpectrums(yf,kf,new_alphaf,0);
      mulSpectrums(kf,kf_lambda,new_alphaf_den,0);
    }else{
      for(int i=0;i<yf.rows;i++){
        for(int j=0;j<yf.cols;j++){
          den = 1.0/(kf_lambda.at<Vec2d>(i,j)[0]*kf_lambda.at<Vec2d>(i,j)[0]+kf_lambda.at<Vec2d>(i,j)[1]*kf_lambda.at<Vec2d>(i,j)[1]);

          new_alphaf.at<Vec2d>(i,j)[0]=
          (yf.at<Vec2d>(i,j)[0]*kf_lambda.at<Vec2d>(i,j)[0]+yf.at<Vec2d>(i,j)[1]*kf_lambda.at<Vec2d>(i,j)[1])*den;
          new_alphaf.at<Vec2d>(i,j)[1]=
          (yf.at<Vec2d>(i,j)[1]*kf_lambda.at<Vec2d>(i,j)[0]-yf.at<Vec2d>(i,j)[0]*kf_lambda.at<Vec2d>(i,j)[1])*den;
        }
      }
    }

(17)式を計算して $\hat{\bf{\alpha}}$ を求めます。
params.split_coeffというパラメータ（デフォルトでtrue）によって若干求め方が変わります。
trueの時は、(17)式の分母と分子を別々に計算しています。(17)式と異なるのは、こちらのコードでは分母が $\hat{\bf{k}}^{\bf{xx}}\odot(\hat{\bf{k}}^{\bf{xx}}+\lambda)$ 、分子が $\hat{\bf{k}}^{\bf{xx}}\odot\hat{\bf{y}}$ と、なぜか分母分子ともに $\hat{\bf{k}}^{\bf{xx}}$ の要素の掛け算を行っています。

そもそもなぜparams.split_coeffというパラメータを用意したのか謎です。

    // update the RLS model
    if(frame==0){
      alphaf=new_alphaf.clone();
      if(params.split_coeff)alphaf_den=new_alphaf_den.clone();
    }else{
      alphaf=(1.0-params.interp_factor)*alphaf+params.interp_factor*new_alphaf;
      if(params.split_coeff)alphaf_den=(1.0-params.interp_factor)*alphaf_den+params.interp_factor*new_alphaf_den;
    }

最初のフレームの場合、求めた $\hat{\bf{\alpha}}$ を保存し、それ以降のフレームの場合は重みを付けて $\hat{\bf{\alpha}}$ を更新します。

余談

今回コードリーディングをして気になったのですが、KCFに関してTracker::init()にはなんの学習機能も実装されておらず、代わりにupdate()によってはじめて学習が実行されます。ということは、以下のサイトのサンプルコードは実は問題あると思われます。

http://docs.opencv.org/3.1.0/d2/d0a/tutorial_introduction_to_tracker.html

ここでは、

1. create()

2. 画像取得＋ROI指定

3. init()

4. 画像取得

5. update()

6. 4に戻る

という処理を行っていますが、この場合最初のupdate()で学習されるテンプレートはinit()時にROIを指定した画像とは別のものが与えられています。

従って正しくは、

1. create()

2. 画像取得＋ROI指定

3. init()

4. update()

5. 画像取得

6. update()

7. 5に戻る

にすべきではないかと思います。

最後に

というわけで、KCFのアルゴリズムとコードについて解説しました。

これは自分に課していた宿題でもあったので、これを機会に実現できてよかったです。

ただ、１回のエントリーとしては長すぎたので、２回に分けるべきでした。

このように論文とコードの間を行ったり来たりすると、より相互の理解も深まる気がします。

かなりマイナーな記事だと思いますが、どなたかのお役立てればうれしいです。

参考文献

[1]Henriques, J., Caseiro, R., Martins, P., & Batista, J. High-Speed Tracking with Kernelized Correlation Filters. IEEE Transactions on Pattern Analysis and Machine Intelligence (TPAMI) 37 (3), 583-596

[2]Kalal, Z., Jiri, M., & Krystian, M. (2010). P-N Learning : Bootstrapping Binary Classifiers by Structural Constraints. In IEEE Conference on Computer Vision and Pattern Recognition (CVPR).

[3]Danelljan, M., Khan, F. S., Felsberg, M., & Weijer, J. Van De. (2014). Adaptive Color Attributes for Real-Time Visual Tracking. In IEEE Conference on Computer Vision and Pattern Recognition (CVPR).