optimization

modules/tracking/include/opencv2/tracking/tracker.hpp

@@ -1204,7 +1204,11 @@ class CV_EXPORTS_W TrackerKCF : public Tracker
     -   "GRAY" -- Use grayscale values as the feature
     -   "CN" -- Color-names feature
    */
-  enum MODE {GRAY, CN, CN2};
+  enum MODE {
+    GRAY = (1u << 0),
+    CN = (1u << 1),
+    CUSTOM = (1u<<2)
+  };
 
   struct CV_EXPORTS Params
   {
@@ -1234,9 +1238,12 @@ class CV_EXPORTS_W TrackerKCF : public Tracker
     bool compress_feature;        //!<  activate the pca method to compress the features
     int max_patch_size;           //!<  threshold for the ROI size
     int compressed_size;          //!<  feature size after compression
-    MODE descriptor;              //!<  descriptor type
+    unsigned int desc_pca;        //!<  compressed descriptors of TrackerKCF::MODE
+    unsigned int desc_npca;       //!<  non-compressed descriptors of TrackerKCF::MODE
   };
 
+  virtual void setFeatureExtractor(void (*)(const Mat, const Rect, Mat&), bool pca_func = false);
+
   /** @brief Constructor
       @param parameters KCF parameters TrackerKCF::Params
   */
@@ -1355,8 +1362,8 @@ Rect2d CV_EXPORTS_W selectROI(Mat img, bool fromCenter = true);
 Rect2d CV_EXPORTS_W selectROI(const std::string& windowName, Mat img, bool showCrossair = true, bool fromCenter = true);
 void CV_EXPORTS_W selectROI(const std::string& windowName, Mat img, std::vector<Rect2d> & boundingBox, bool fromCenter = true);
 
-//! @}
-
 } /* namespace cv */
 
+//! @}
+
 #endif

modules/tracking/src/trackerKCF.cpp

@@ -73,6 +73,7 @@ namespace cv{
     TrackerKCFImpl( const TrackerKCF::Params &parameters = TrackerKCF::Params() );
     void read( const FileNode& /*fn*/ );
     void write( FileStorage& /*fs*/ ) const;
+    void setFeatureExtractor(void (*f)(const Mat, const Rect, Mat&), bool pca_func = false);
 
   protected:
      /*
@@ -86,19 +87,22 @@ namespace cv{
     /*
     * KCF functions and vars
     */
-    void createHanningWindow(OutputArray _dst, const cv::Size winSize, const int type) const;
+    void createHanningWindow(OutputArray dest, const cv::Size winSize, const int type) const;
-    void inline fft2(const Mat src, std::vector<Mat> & dest) const;
+    void inline fft2(const Mat src, std::vector<Mat> & dest, std::vector<Mat> & layers_data) const;
     void inline fft2(const Mat src, Mat & dest) const;
     void inline ifft2(const Mat src, Mat & dest) const;
     void inline pixelWiseMult(const std::vector<Mat> src1, const std::vector<Mat>  src2, std::vector<Mat>  & dest, const int flags, const bool conjB=false) const;
     void inline sumChannels(std::vector<Mat> src, Mat & dest) const;
-    void inline updateProjectionMatrix(const Mat src, Mat & old_cov,Mat &  _proj_mtx,double pca_rate, int compressed_sz) const;
+    void inline updateProjectionMatrix(const Mat src, Mat & old_cov,Mat &  proj_matrix,double pca_rate, int compressed_sz,
-    void inline compress(const Mat _proj_mtx, const Mat src, Mat & dest) const;
+                                       std::vector<Mat> & layers_pca,std::vector<Scalar> & average, Mat pca_data, Mat new_cov, Mat w, Mat u, Mat v) const;
-    bool getSubWindow(const Mat img, const Rect roi, Mat& patch) const;
+    void inline compress(const Mat proj_matrix, const Mat src, Mat & dest, Mat & data, Mat & compressed) const;
-    void extractCN(Mat _patch, Mat & cnFeatures) const;
+    bool getSubWindow(const Mat img, const Rect roi, Mat& feat, Mat& patch, TrackerKCF::MODE desc = GRAY) const;
-    void denseGaussKernel(const double sigma, const Mat _x, const Mat _y, Mat & _k) const;
+    bool getSubWindow(const Mat img, const Rect roi, Mat& feat, void (*f)(const Mat, const Rect, Mat& )) const;
-    void calcResponse(const Mat _alphaf, const Mat _k, Mat & _response) const;
+    void extractCN(Mat patch_data, Mat & cnFeatures) const;
-    void calcResponse(const Mat _alphaf, const Mat _alphaf_den, const Mat _k, Mat & _response) const;
+    void denseGaussKernel(const double sigma, const Mat , const Mat y_data, Mat & k_data,
+                          std::vector<Mat> & layers_data,std::vector<Mat> & xf_data,std::vector<Mat> & yf_data, std::vector<Mat> xyf_v, Mat xy, Mat xyf ) const;
+    void calcResponse(const Mat alphaf_data, const Mat kf_data, Mat & response_data, Mat & spec_data) const;
+    void calcResponse(const Mat alphaf_data, const Mat alphaf_den_data, const Mat kf_data, Mat & response_data, Mat & spec_data, Mat & spec2_data) const;
 
     void shiftRows(Mat& mat) const;
     void shiftRows(Mat& mat, int n) const;
@@ -108,17 +112,46 @@ namespace cv{
     double output_sigma;
     Rect2d roi;
     Mat hann; 	//hann window filter
+    Mat hann_cn; //10 dimensional hann-window filter for CN features,
 
     Mat y,yf; 	// training response and its FFT
-    Mat x,xf; 	// observation and its FFT
+    Mat x; 	// observation and its FFT
     Mat k,kf;	// dense gaussian kernel and its FFT
     Mat kf_lambda; // kf+lambda
     Mat new_alphaf, alphaf;	// training coefficients
     Mat new_alphaf_den, alphaf_den; // for splitted training coefficients
-    Mat z, new_z; // model
+    Mat z; // model
     Mat response; // detection result
     Mat old_cov_mtx, proj_mtx; // for feature compression
 
+    // pre-defined Mat variables for optimization of private functions
+    Mat spec, spec2;
+    std::vector<Mat> layers;
+    std::vector<Mat> vxf,vyf,vxyf;
+    Mat xy_data,xyf_data;
+    Mat data_temp, compress_data;
+    std::vector<Mat> layers_pca_data;
+    std::vector<Scalar> average_data;
+    Mat img_Patch;
+
+    // storage for the extracted features, KRLS model, KRLS compressed model
+    Mat X[2],Z[2],Zc[2];
+
+    // storage of the extracted features
+    std::vector<Mat> features_pca;
+    std::vector<Mat> features_npca;
+    std::vector<MODE> descriptors_pca;
+    std::vector<MODE> descriptors_npca;
+
+    // optimization variables for updateProjectionMatrix
+    Mat data_pca, new_covar,w_data,u_data,vt_data;
+
+    // custom feature extractor
+    bool use_custom_extractor_pca;
+    bool use_custom_extractor_npca;
+    std::vector<void(*)(const Mat img, const Rect roi, Mat& output)> extractor_pca;
+    std::vector<void(*)(const Mat img, const Rect roi, Mat& output)> extractor_npca;
+
     bool resizeImage; // resize the image whenever needed and the patch size is large
 
     int frame;
@@ -135,8 +168,9 @@ namespace cv{
   {
     isInit = false;
     resizeImage = false;
+    use_custom_extractor_pca = false;
+    use_custom_extractor_npca = false;
 
-    CV_Assert(params.descriptor == GRAY || params.descriptor == CN /*|| params.descriptor == CN2*/);
   }
 
   void TrackerKCFImpl::read( const cv::FileNode& fn ){
@@ -179,10 +213,10 @@ namespace cv{
 
     // initialize the hann window filter
     createHanningWindow(hann, roi.size(), CV_64F);
-    if(params.descriptor==CN){
+
-      Mat layers[] = {hann, hann, hann, hann, hann, hann, hann, hann, hann, hann};
+    // hann window filter for CN feature
-      merge(layers, 10, hann);
+    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);
@@ -200,6 +234,28 @@ namespace cv{
 
     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;
   }
@@ -210,33 +266,71 @@ namespace cv{
   bool TrackerKCFImpl::updateImpl( const Mat& image, Rect2d& boundingBox ){
     double minVal, maxVal;	// min-max response
     Point minLoc,maxLoc;	// min-max location
-    Mat zc;
 
     Mat img=image.clone();
     // check the channels of the input image, grayscale is preferred
-    CV_Assert(image.channels() == 1 || image.channels() == 3);
+    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));
 
-    // extract and pre-process the patch
-    if(!getSubWindow(img,roi, x))return false;
-
     // 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
-      if(params.compress_feature){
+      denseGaussKernel(params.sigma,x,z,k,layers,vxf,vyf,vxyf,xy_data,xyf_data);
-        compress(proj_mtx,x,x);
+
-        compress(proj_mtx,z,zc);
+      // compute the fourier transform of the kernel
-        denseGaussKernel(params.sigma,x,zc,k);
+      fft2(k,kf);
-      }else
+      if(frame==1)spec2=Mat_<Vec2d >(kf.rows, kf.cols);
-        denseGaussKernel(params.sigma,x,z,k);
 
       // calculate filter response
       if(params.split_coeff)
-        calcResponse(alphaf,alphaf_den,k,response);
+        calcResponse(alphaf,alphaf_den,kf,response, spec, spec2);
       else
-        calcResponse(alphaf,k,response);
+        calcResponse(alphaf,kf,response, spec);
 
       // extract the maximum response
       minMaxLoc( response, &minVal, &maxVal, &minLoc, &maxLoc );
@@ -249,43 +343,84 @@ namespace cv{
     }
 
     // extract the patch for learning purpose
-    if(!getSubWindow(img,roi, x))return false;
+    // 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
-    new_z=x.clone();
+    if(frame==0){
-    if(frame==0)
+      Z[0] = X[0].clone();
-      z=x.clone();
+      Z[1] = X[1].clone();
-    else
+    }else{
-      z=(1.0-params.interp_factor)*z+params.interp_factor*new_z;
+      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());
+      }
 
-    if(params.compress_feature){
       // feature compression
-      updateProjectionMatrix(z,old_cov_mtx,proj_mtx,params.pca_learning_rate,params.compressed_size);
+      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,x);
+      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);
+    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;
 
-    /* TODO: optimize this element-wise division
+    double den;
-     * new_alphaf=yf./kf
-     * z=(a+bi)/(c+di)[(ac+bd)+i(bc-ad)]/(c^2+d^2)
-     */
-    new_alphaf=Mat_<Vec2d >(yf.rows, yf.cols);
-    std::complex<double> temp;
-
     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++){
-          temp=std::complex<double>(yf.at<Vec2d>(i,j)[0],yf.at<Vec2d>(i,j)[1])/(std::complex<double>(kf_lambda.at<Vec2d>(i,j)[0],kf_lambda.at<Vec2d>(i,j)[1])/*+std::complex<double>(0.0000000001,0.0000000001)*/);
+          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]=temp.real();
+
-          new_alphaf.at<Vec2d >(i,j)[1]=temp.imag();
+          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;
         }
       }
     }
@@ -311,11 +446,11 @@ namespace cv{
   /*
    * hann window filter
    */
-  void TrackerKCFImpl::createHanningWindow(OutputArray _dst, const cv::Size winSize, const int type) const {
+  void TrackerKCFImpl::createHanningWindow(OutputArray dest, const cv::Size winSize, const int type) const {
       CV_Assert( type == CV_32FC1 || type == CV_64FC1 );
 
-      _dst.create(winSize, type);
+      dest.create(winSize, type);
-      Mat dst = _dst.getMat();
+      Mat dst = dest.getMat();
 
       int rows = dst.rows, cols = dst.cols;
 
@@ -350,27 +485,14 @@ namespace cv{
    * simplification of fourier transform function in opencv
    */
   void inline TrackerKCFImpl::fft2(const Mat src, Mat & dest) const {
-    std::vector<Mat> layers(src.channels());
+    dft(src,dest,DFT_COMPLEX_OUTPUT);
-    std::vector<Mat> outputs(src.channels());
-
-    split(src, layers);
-
-    for(int i=0;i<src.channels();i++){
-      dft(layers[i],outputs[i],DFT_COMPLEX_OUTPUT);
-    }
-
-    merge(outputs,dest);
   }
 
-  void inline TrackerKCFImpl::fft2(const Mat src, std::vector<Mat> & dest) const {
+  void inline TrackerKCFImpl::fft2(const Mat src, std::vector<Mat> & dest, std::vector<Mat> & layers_data) const {
-    std::vector<Mat> layers(src.channels());
+    split(src, layers_data);
-    dest.clear();
-    dest.resize(src.channels());
-
-    split(src, layers);
 
     for(int i=0;i<src.channels();i++){
-      dft(layers[i],dest[i],DFT_COMPLEX_OUTPUT);
+      dft(layers_data[i],dest[i],DFT_COMPLEX_OUTPUT);
     }
   }
 
@@ -385,9 +507,6 @@ namespace cv{
    * Point-wise multiplication of two Multichannel Mat data
    */
   void inline TrackerKCFImpl::pixelWiseMult(const std::vector<Mat> src1, const std::vector<Mat>  src2, std::vector<Mat>  & dest, const int flags, const bool conjB) const {
-    dest.clear();
-    dest.resize(src1.size());
-
     for(unsigned i=0;i<src1.size();i++){
       mulSpectrums(src1[i], src2[i], dest[i],flags,conjB);
     }
@@ -406,55 +525,51 @@ namespace cv{
   /*
    * obtains the projection matrix using PCA
    */
-  void inline TrackerKCFImpl::updateProjectionMatrix(const Mat src, Mat & old_cov,Mat &  _proj_mtx, double pca_rate, int compressed_sz) const {
+  void inline TrackerKCFImpl::updateProjectionMatrix(const Mat src, Mat & old_cov,Mat &  proj_matrix, double pca_rate, int compressed_sz,
+                                                     std::vector<Mat> & layers_pca,std::vector<Scalar> & average, Mat pca_data, Mat new_cov, Mat w, Mat u, Mat vt) const {
     CV_Assert(compressed_sz<=src.channels());
 
-    // compute average
+    split(src,layers_pca);
-    std::vector<Mat> layers(src.channels());
-    std::vector<Scalar> average(src.channels());
-    split(src,layers);
 
     for (int i=0;i<src.channels();i++){
-      average[i]=mean(layers[i]);
+      average[i]=mean(layers_pca[i]);
-      layers[i]-=average[i];
+      layers_pca[i]-=average[i];
     }
 
     // calc covariance matrix
-    Mat data,new_cov;
+    merge(layers_pca,pca_data);
-    merge(layers,data);
+    pca_data=pca_data.reshape(1,src.rows*src.cols);
-    data=data.reshape(1,src.rows*src.cols);
 
-    new_cov=1.0/(double)(src.rows*src.cols-1)*(data.t()*data);
+    new_cov=1.0/(double)(src.rows*src.cols-1)*(pca_data.t()*pca_data);
     if(old_cov.rows==0)old_cov=new_cov.clone();
 
     // calc PCA
-    Mat w, u, vt;
     SVD::compute((1.0-pca_rate)*old_cov+pca_rate*new_cov, w, u, vt);
 
     // extract the projection matrix
-    _proj_mtx=u(Rect(0,0,compressed_sz,src.channels())).clone();
+    proj_matrix=u(Rect(0,0,compressed_sz,src.channels())).clone();
-    Mat proj_vars=Mat::eye(compressed_sz,compressed_sz,_proj_mtx.type());
+    Mat proj_vars=Mat::eye(compressed_sz,compressed_sz,proj_matrix.type());
     for(int i=0;i<compressed_sz;i++){
       proj_vars.at<double>(i,i)=w.at<double>(i);
     }
 
     // update the covariance matrix
-    old_cov=(1.0-pca_rate)*old_cov+pca_rate*_proj_mtx*proj_vars*_proj_mtx.t();
+    old_cov=(1.0-pca_rate)*old_cov+pca_rate*proj_matrix*proj_vars*proj_matrix.t();
   }
 
   /*
    * compress the features
    */
-  void inline TrackerKCFImpl::compress(const Mat _proj_mtx, const Mat src, Mat & dest) const {
+  void inline TrackerKCFImpl::compress(const Mat proj_matrix, const Mat src, Mat & dest, Mat & data, Mat & compressed) const {
-    Mat data=src.reshape(1,src.rows*src.cols);
+    data=src.reshape(1,src.rows*src.cols);
-    Mat compressed=data*_proj_mtx;
+    compressed=data*proj_matrix;
-    dest=compressed.reshape(_proj_mtx.cols,src.rows).clone();
+    dest=compressed.reshape(proj_matrix.cols,src.rows).clone();
   }
 
   /*
    * obtain the patch and apply hann window filter to it
    */
-  bool TrackerKCFImpl::getSubWindow(const Mat img, const Rect _roi, Mat& patch) const {
+  bool TrackerKCFImpl::getSubWindow(const Mat img, const Rect _roi, Mat& feat, Mat& patch, TrackerKCF::MODE desc) const {
 
     Rect region=_roi;
 
@@ -486,77 +601,108 @@ namespace cv{
     if(patch.rows==0 || patch.cols==0)return false;
 
     // extract the desired descriptors
-    switch(params.descriptor){
+    switch(desc){
-      case GRAY:
-        if(img.channels()>1)cvtColor(patch,patch, CV_BGR2GRAY);
-        patch.convertTo(patch,CV_64F);
-        patch=patch/255.0-0.5; // normalize to range -0.5 .. 0.5
-        break;
       case CN:
         CV_Assert(img.channels() == 3);
-        extractCN(patch,patch);
+        extractCN(patch,feat);
+        feat=feat.mul(hann_cn); // hann window filter
         break;
-      case CN2:
+      default: // GRAY
-        if(patch.channels()>1)cvtColor(patch,patch, CV_BGR2GRAY);
+        if(img.channels()>1)
+          cvtColor(patch,feat, CV_BGR2GRAY);
+        else
+          feat=patch;
+        feat.convertTo(feat,CV_64F);
+        feat=feat/255.0-0.5; // normalize to range -0.5 .. 0.5
+        feat=feat.mul(hann); // hann window filter
         break;
     }
 
-    patch=patch.mul(hann); // hann window filter
-
     return true;
 
   }
 
+  /*
+   * get feature using external function
+   */
+  bool TrackerKCFImpl::getSubWindow(const Mat img, const Rect _roi, Mat& feat, void (*f)(const Mat, const Rect, Mat& )) const{
+
+    // return false if roi is outside the image
+    if((_roi.x+_roi.width<0)
+      ||(_roi.y+_roi.height<0)
+      ||(_roi.x>=img.cols)
+      ||(_roi.y>=img.rows)
+    )return false;
+
+    f(img, _roi, feat);
+
+    if(_roi.width != feat.cols || _roi.height != feat.rows){
+      printf("error in customized function of features extractor!\n");
+      printf("Rules: roi.width==feat.cols && roi.height = feat.rows \n");
+    }
+
+    Mat hann_win;
+    std::vector<Mat> _layers;
+
+    for(int i=0;i<feat.channels();i++)
+      _layers.push_back(hann);
+
+    merge(_layers, hann_win);
+
+    feat=feat.mul(hann_win); // hann window filter
+
+    return true;
+  }
+
   /* Convert BGR to ColorNames
    */
-  void TrackerKCFImpl::extractCN(Mat _patch, Mat & cnFeatures) const {
+  void TrackerKCFImpl::extractCN(Mat patch_data, Mat & cnFeatures) const {
-    Vec3b & pixel = _patch.at<Vec3b>(0,0);
+    Vec3b & pixel = patch_data.at<Vec3b>(0,0);
     unsigned index;
 
-    Mat temp = Mat::zeros(_patch.rows,_patch.cols,CV_64FC(10));
+    if(cnFeatures.type() != CV_64FC(10))
+      cnFeatures = Mat::zeros(patch_data.rows,patch_data.cols,CV_64FC(10));
 
-    for(int i=0;i<_patch.rows;i++){
+    for(int i=0;i<patch_data.rows;i++){
-      for(int j=0;j<_patch.cols;j++){
+      for(int j=0;j<patch_data.cols;j++){
-        pixel=_patch.at<Vec3b>(i,j);
+        pixel=patch_data.at<Vec3b>(i,j);
         index=(unsigned)(floor(pixel[2]/8)+32*floor(pixel[1]/8)+32*32*floor(pixel[0]/8));
 
         //copy the values
         for(int _k=0;_k<10;_k++){
-          temp.at<Vec<double,10> >(i,j)[_k]=ColorNames[index][_k];
+          cnFeatures.at<Vec<double,10> >(i,j)[_k]=ColorNames[index][_k];
         }
       }
     }
 
-    cnFeatures=temp.clone();
   }
 
   /*
    *  dense gauss kernel function
    */
-  void TrackerKCFImpl::denseGaussKernel(const double sigma, const Mat _x, const Mat _y, Mat & _k) const {
+  void TrackerKCFImpl::denseGaussKernel(const double sigma, const Mat x_data, const Mat y_data, Mat & k_data,
-    std::vector<Mat> _xf,_yf,xyf_v;
+                                        std::vector<Mat> & layers_data,std::vector<Mat> & xf_data,std::vector<Mat> & yf_data, std::vector<Mat> xyf_v, Mat xy, Mat xyf ) const {
-    Mat xy,xyf;
     double normX, normY;
 
-    fft2(_x,_xf);
+    fft2(x_data,xf_data,layers_data);
-    fft2(_y,_yf);
+    fft2(y_data,yf_data,layers_data);
 
-    normX=norm(_x);
+    normX=norm(x_data);
     normX*=normX;
-    normY=norm(_y);
+    normY=norm(y_data);
     normY*=normY;
 
-    pixelWiseMult(_xf,_yf,xyf_v,0,true);
+    pixelWiseMult(xf_data,yf_data,xyf_v,0,true);
     sumChannels(xyf_v,xyf);
     ifft2(xyf,xyf);
 
     if(params.wrap_kernel){
-      shiftRows(xyf, _x.rows/2);
+      shiftRows(xyf, x_data.rows/2);
-      shiftCols(xyf, _x.cols/2);
+      shiftCols(xyf, x_data.cols/2);
     }
 
     //(xx + yy - 2 * xy) / numel(x)
-    xy=(normX+normY-2*xyf)/(_x.rows*_x.cols*_x.channels());
+    xy=(normX+normY-2*xyf)/(x_data.rows*x_data.cols*x_data.channels());
 
     // TODO: check wether we really need thresholding or not
     //threshold(xy,xy,0.0,0.0,THRESH_TOZERO);//max(0, (xx + yy - 2 * xy) / numel(x))
@@ -568,7 +714,7 @@ namespace cv{
 
     double sig=-1.0/(sigma*sigma);
     xy=sig*xy;
-    exp(xy,_k);
+    exp(xy,k_data);
 
   }
 
@@ -626,36 +772,42 @@ namespace cv{
   /*
    * calculate the detection response
    */
-  void TrackerKCFImpl::calcResponse(const Mat _alphaf, const Mat _k, Mat & _response) const {
+  void TrackerKCFImpl::calcResponse(const Mat alphaf_data, const Mat kf_data, Mat & response_data, Mat & spec_data) const {
     //alpha f--> 2channels ; k --> 1 channel;
-    Mat _kf;
+    mulSpectrums(alphaf_data,kf_data,spec_data,0,false);
-    fft2(_k,_kf);
+    ifft2(spec_data,response_data);
-    Mat spec;
-    mulSpectrums(_alphaf,_kf,spec,0,false);
-    ifft2(spec,_response);
   }
 
   /*
    * calculate the detection response for splitted form
    */
-  void TrackerKCFImpl::calcResponse(const Mat _alphaf, const Mat _alphaf_den, const Mat _k, Mat & _response) const {
+  void TrackerKCFImpl::calcResponse(const Mat alphaf_data, const Mat _alphaf_den, const Mat kf_data, Mat & response_data, Mat & spec_data, Mat & spec2_data) const {
-    Mat _kf;
-    fft2(_k,_kf);
-    Mat spec;
-    Mat spec2=Mat_<Vec2d >(_k.rows, _k.cols);
-    std::complex<double> temp;
 
-    mulSpectrums(_alphaf,_kf,spec,0,false);
+    mulSpectrums(alphaf_data,kf_data,spec_data,0,false);
 
-    for(int i=0;i<_k.rows;i++){
+    //z=(a+bi)/(c+di)=[(ac+bd)+i(bc-ad)]/(c^2+d^2)
-      for(int j=0;j<_k.cols;j++){
+    double den;
-        temp=std::complex<double>(spec.at<Vec2d>(i,j)[0],spec.at<Vec2d>(i,j)[1])/(std::complex<double>(_alphaf_den.at<Vec2d>(i,j)[0],_alphaf_den.at<Vec2d>(i,j)[1])/*+std::complex<double>(0.0000000001,0.0000000001)*/);
+    for(int i=0;i<kf_data.rows;i++){
-        spec2.at<Vec2d >(i,j)[0]=temp.real();
+      for(int j=0;j<kf_data.cols;j++){
-        spec2.at<Vec2d >(i,j)[1]=temp.imag();
+        den=1.0/(_alphaf_den.at<Vec2d>(i,j)[0]*_alphaf_den.at<Vec2d>(i,j)[0]+_alphaf_den.at<Vec2d>(i,j)[1]*_alphaf_den.at<Vec2d>(i,j)[1]);
+        spec2_data.at<Vec2d>(i,j)[0]=
+          (spec_data.at<Vec2d>(i,j)[0]*_alphaf_den.at<Vec2d>(i,j)[0]+spec_data.at<Vec2d>(i,j)[1]*_alphaf_den.at<Vec2d>(i,j)[1])*den;
+        spec2_data.at<Vec2d>(i,j)[1]=
+          (spec_data.at<Vec2d>(i,j)[1]*_alphaf_den.at<Vec2d>(i,j)[0]-spec_data.at<Vec2d>(i,j)[0]*_alphaf_den.at<Vec2d>(i,j)[1])*den;
       }
     }
 
-    ifft2(spec2,_response);
+    ifft2(spec2_data,response_data);
+  }
+
+  void TrackerKCFImpl::setFeatureExtractor(void (*f)(const Mat, const Rect, Mat&), bool pca_func){
+    if(pca_func){
+      extractor_pca.push_back(f);
+      use_custom_extractor_pca = true;
+    }else{
+      extractor_npca.push_back(f);
+      use_custom_extractor_npca = true;
+    }
   }
   /*----------------------------------------------------------------------*/
 
@@ -669,9 +821,10 @@ namespace cv{
       output_sigma_factor=1.0/16.0;
       resize=true;
       max_patch_size=80*80;
-      descriptor=CN;
       split_coeff=true;
       wrap_kernel=false;
+      desc_npca = GRAY;
+      desc_pca = CN;
 
       //feature compression
       compress_feature=true;
@@ -683,4 +836,6 @@ namespace cv{
 
   void TrackerKCF::Params::write( cv::FileStorage& /*fs*/ ) const{}
 
+  void TrackerKCF::setFeatureExtractor(void (*)(const Mat, const Rect, Mat&), bool ){};
+
 } /* namespace cv */