BRIEF算法学习

http://blog.sina.com.cn/s/blog_916b71bb0100wbdh.html    算法的中文简介

http://cvlab.epfl.ch/research/detect/brief/   算法的homepage

这个算法比较简单，

1.利用已有的算法（比如说fast-9算法）找到角点，在作者的示例程序中使用surf算法来搜寻角点。

2.在找到的角点的周围规定一个S*S的patch，对patch所在的S*S区域进行高斯模糊。

3.论文中介绍了5种算法对patch中的找到nd个（x，y)个像素对（其中x and y代表一个像素点的位置）进行取样。其中第二种算法，作者认为是最好的，并且在给中的代码中仅给出了这种算法。

翻译：（X,Y)均服从均值为0，方差为1/25*s*s的正太分布，实验是在各向同性的高斯分布下抽样的。我们发现：在方差为1/25*s*s参数下可以取得最好的识别率。

4.抽样完毕后根据：

计算出每个tao，然后根据：

计算出特征描述子。特征描述子是一个binary序列。

5.在进行匹配的时候使用hamming距离计算相似度。

论文代码分析：

TestSampler< TestLocT >::sample(tests_, DESC_LEN, INIT_PATCH_SZ, 1);

下面这个函数产生高斯取样。

double randNormal(const double std, const double mu){double x1, x2, w, y1;static double y2;static bool use_last = false;if (use_last)        /* use value from previous call */{y1 = y2;use_last = false;}else{do {x1 = 2.0 * cvRandReal(&LocalStorage::getInst().rng_obj) - 1.0;x2 = 2.0 * cvRandReal(&LocalStorage::getInst().rng_obj) - 1.0;w = x1 * x1 + x2 * x2;} while ( w >= 1.0 );w = sqrt( (-2.0 * log( w ) ) / w );y1 = x1 * w;y2 = x2 * w;use_last = true;}return mu + y1 * std;}} // namespace utils

下面的代码产生brief描述子。把描述子放在res中。

template< typename KPT_T >void BRIEF::getBRIEF(const std::vector< KPT_T >& kpts,                     std::vector< std::bitset<DESC_LEN> >& res,                     const IplImage* precomp_int_img){   assert(kpts[0].image->nChannels == 1);   // Make sure all points are from the same image   const IplImage* img = kpts[0].image;   assert(img);   #if DEBUG      for (int i=1; i<(int)kpts.size(); ++i)         assert(kpts[i].image == img);   #endif   // Pre-compute integral image//   #if USE_INTEGRAL_IMAGE//      if (!precomp_int_img && cur_int_img_.src != img)      // Update needed for our integral image?//      {//         if (!cur_int_img_.src)//            cur_int_img_.data = new float[img->width*img->height];//         else if (cur_int_img_.src->width*cur_int_img_.src->height < img->width*img->height)//         {//            delete [] cur_int_img_.data;//            cur_int_img_.data = new float[img->width*img->height];//         }////         cur_int_img_.src = img;//         utils::computeIntegralImage(img, cur_int_img_.data);//      }//   #endif   // Recompute bounding box, if needed   if (act_patch_sz_ <= 0)      recomputeActPatchSz();   // Compute all descriptors   std::bitset<DESC_LEN> some_desc;   res.reserve(kpts.size());   for (int i=0; i<(int)kpts.size(); ++i)   {      const KPT_T &kp = kpts[i];      #if DEBUG         if (!checkAllTestInsideImage(kp)) {            some_desc.reset();            res.push_back(some_desc);            continue;         }      #endif      // Rotate and scale tests, if information in keypoints valid      // Note: Mind the INIT_PATCH_SZ and BOX_SZ params or this will segfault      #if SOS_SCALE != SOS_DO_NOT_INTERPRET && SOS_ORI != SOS_DO_NOT_INTERPRET         rotateAndScaleTests(interpretOri(kp.ori), interpretScale(kp.scale), true);         #define KERNEL_SZ act_kernel_sz_      #elif SOS_SCALE != SOS_DO_NOT_INTERPRET         scaleTests(interpretScale(kp.scale), true);         #define KERNEL_SZ act_kernel_sz_      #elif SOS_ORI != SOS_DO_NOT_INTERPRET         rotateTests(interpretOri(kp.ori), true);         #define KERNEL_SZ INIT_KERNEL_SZ      #else         #define KERNEL_SZ INIT_KERNEL_SZ      #endif      // Smooth patch, compute descriptor      #if USE_INTEGRAL_IMAGE         assert(!precomp_int_img || precomp_int_img->depth == IPL_DEPTH_64F);         const float* ii_data = cur_int_img_.data;         const int ii_step = cur_int_img_.src->width;         #if SMOOTH_ENTIRE_PATCH    // smooth entire patch, then apply tests            #if DEBUG               if (!checkSmoothAreaInsideIntImage(kp)) {                  some_desc.reset();                  res.push_back(some_desc);                  continue;               }            #endif            // Compute UNNORMALIZED smoothed patch around (kp.x, kp.y), ok since only            // relative values are of importance            // Using "area = B + D - A - C" where arrangement is  B  A            // 对patch区域进行高斯模糊。                                                   C  D            float *p_smooth = smooth_;            for (int v=0; v<act_patch_sz_; ++v)            {               // cix: smoothing window central pixel offset               int cix = (kp.y-act_patch_sz_/2+v)*ii_step + kp.x-act_patch_sz_/2;               const float *p_a = &ii_data[cix - KERNEL_SZ/2*ii_step + KERNEL_SZ/2];               const float *p_b = &ii_data[cix - KERNEL_SZ/2*ii_step - KERNEL_SZ/2];               const float *p_c = &ii_data[cix + KERNEL_SZ/2*ii_step - KERNEL_SZ/2];               const float *p_d = &ii_data[cix + KERNEL_SZ/2*ii_step + KERNEL_SZ/2];               for (int u=0; u<act_patch_sz_; ++u, ++p_smooth, ++p_a, ++p_b, ++p_c, ++p_d)                  *p_smooth = ((*p_b + *p_d - *p_a - *p_c));            }            // Ready for tests            res.push_back(some_desc);            applyTests(smooth_, res[res.size()-1], tests_);         #else    // smooth around test locations only            #if DEBUG               if (!checkSmoothAreaForTestInsideIntImage(kp)) {                  some_desc.reset();                  res.push_back(some_desc);                  continue;               }            #endif            res.push_back(some_desc);            std::bitset<DESC_LEN> &bits = res[res.size()-1];            const TestLocT *pxy = tests_;//获得高斯采样的的数组，这样可以获得高的速度，一次计算后，以后查表。            for (int j=0; j<DESC_LEN; ++j, pxy+=4)            {               // cix: central pixel offset               const int cix_1 = (kp.y + pxy[1])*ii_step + kp.x+pxy[0];               const int cix_2 = (kp.y + pxy[3])*ii_step + kp.x+pxy[2];               // Compute tests按照公式（1）计算。               bits[j] = ( ii_data[cix_1 - KERNEL_SZ/2*ii_step - KERNEL_SZ/2] +                           ii_data[cix_1 + KERNEL_SZ/2*ii_step + KERNEL_SZ/2] -                           ii_data[cix_1 - KERNEL_SZ/2*ii_step + KERNEL_SZ/2] -                           ii_data[cix_1 + KERNEL_SZ/2*ii_step - KERNEL_SZ/2] )                         >                         ( ii_data[cix_2 - KERNEL_SZ/2*ii_step - KERNEL_SZ/2] +                           ii_data[cix_2 + KERNEL_SZ/2*ii_step + KERNEL_SZ/2] -                           ii_data[cix_2 - KERNEL_SZ/2*ii_step + KERNEL_SZ/2] -                           ii_data[cix_2 + KERNEL_SZ/2*ii_step - KERNEL_SZ/2] );            }         #endif      #else    // not using integral image but cvSmooth instead         #ifdef DEBUG            if (!checkPatchInsideImage(kp)) {               some_desc.reset();               res.push_back(some_desc);               continue;            }         #endif         // Need to copy patch data before smoothing         assert(img->depth == (int)IPL_DEPTH_8U || img->depth == (int)IPL_DEPTH_8S);  // or memcpy will fail         const int step = img->widthStep;         uchar* start = (uchar*)&img->imageData[(kp.y-act_patch_sz_/2)*step + kp.x-act_patch_sz_/2];         for (int j=0; j<act_patch_sz_; ++j, start += step)            memcpy(patch_8u_->imageData + j*act_patch_sz_, start, act_patch_sz_);         cvSmooth(patch_8u_, patch_8u_2_, CV_BLUR, KERNEL_SZ, KERNEL_SZ);         // Ready for tests         res.push_back(some_desc);         applyTests((uchar*)patch_8u_2_->imageData, res[res.size()-1], tests_);      #endif   }   #undef KERNEL_SZ}

下面的代码对计算出的描述子进行match，以下代码本人没有看懂。

template< typename KPT_T, int DESC_LEN >struct BRIEFMatcher{   // Matching with L/R check   void matchLeftRight(const std::vector< std::bitset< DESC_LEN > > feat_left,                       const std::vector< std::bitset< DESC_LEN > > feat_right,                       const std::vector< KPT_T > kpts_left,                       const std::vector< KPT_T > kpts_right,                       std::insert_iterator< std::vector< KPT_T > > lmatch_ins,                       std::insert_iterator< std::vector< KPT_T > > rmatch_ins)   {      int *dist = new int[MAX(feat_left.size(), feat_right.size())];      int *left_nn = new int[feat_left.size()];      int *right_nn = new int[feat_right.size()];      // Matching with left-right check      int min_ix;      for (int i=0; i<(int)feat_left.size(); ++i)      {         for (int k=0; k<(int)feat_right.size(); ++k)            dist[k] = (int)(feat_left[i] ^ feat_right[k]).count();         utils::minVect(dist, (int)feat_right.size(), &min_ix);         left_nn[i] = min_ix;      }      for (int i=0; i<(int)feat_right.size(); ++i)      {         for (int k=0; k<(int)feat_left.size(); ++k)            dist[k] = (int)(feat_right[i] ^ feat_left[k]).count();         utils::minVect(dist, (int)feat_left.size(), &min_ix);         right_nn[i] = min_ix;      }      for (int i=0; i<(int)feat_left.size(); ++i)         if (right_nn[left_nn[i]] == i)         {            *lmatch_ins = kpts_left[i];            *rmatch_ins = kpts_right[left_nn[i]];            ++lmatch_ins; ++rmatch_ins;         }      delete [] dist;      delete [] left_nn;      delete [] right_nn;   }};