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;
   }
};