SIFT源码分析

        SIFT的原理以及逻辑过程我就不细说了，网上有很多的教程大家可以参考，今天我主要是对SIFT的源码进行细致的分析，包括代码中的各种细节也都会一一讲解。

        我是先贴代码然后做解释

  std::unique_ptr<Regions> Describe(
    const image::Image<unsigned char>& image,
    const image::Image<unsigned char>* mask = nullptr
  ) override
  {
    return DescribeSIFT(image, mask);
  }

        那我们就从MVG调到SIFT的这段代码 return DescribeSIFT这里开始讲起....

        首先进入到 DescribeSIFT这个函数，先看它的的参数，

std::unique_ptr<Regions_type> DescribeSIFT(
    const image::Image<unsigned char>& image,
    const image::Image<unsigned char>* mask = nullptr
)

        它的输入参数有两个，一个是我们要处理的图像，另一个就是掩膜，我们一般用不到值为nullptr。

const int w = image.Width(), h = image.Height();
//Convert to float
const image::Image<float> If(image.GetMat().cast<float>());

VlSiftFilt *filt = vl_sift_new(w, h,
  _params._num_octaves, _params._num_scales, _params._first_octave);
if (_params._edge_threshold >= 0)
  vl_sift_set_edge_thresh(filt, _params._edge_threshold);
if (_params._peak_threshold >= 0)
  vl_sift_set_peak_thresh(filt, 255*_params._peak_threshold/_params._num_scales);

Descriptor<vl_sift_pix, 128> descr;
Descriptor<unsigned char, 128> descriptor;

        这里定义的w和h是图像的宽和长，然后把它转换成float型，便于更精确的计算。接下来我们需要创建一个SIFT过滤器，过滤器的作用就是检测图像中的关键点，并计算这些关键点的SIFT描述符。

        vl_sift_new这个函数有五个参数，w和h是图像的宽度和高度，_params._num_octaves是要计算的尺度空间金字塔的八度数量，_params._num_scales是在每个八度中，除了第一个和最后一个尺度外，要计算的中间尺度的数量（因为第一个和最后一个尺度是不计算极值的）。每个八度实际上包含_num_scales+2个尺度（包括第一个和最后一个通过高斯模糊得到的尺度）。_params._first_octave是尺度空间金字塔中第一个八度的索引。

        下面那两个if语句分别用来设置边缘阈值和峰值阈值，那为什么要设置这两个阈值，我们要知道边缘阈值用于在特征检测过程中排除边缘响应过强的点，因为边缘上的点通常不是稳定的特征点。峰值阈值则是一个直接应用于特征点响应值的阈值，用于确定哪些点应该被视为关键点。

Descriptor<vl_sift_pix, 128> descr;
Descriptor<unsigned char, 128> descriptor;

// Process SIFT computation
vl_sift_process_first_octave(filt, If.data());

        在这段代码中，首先使用了两种不同的数据类型来定义SIFT描述符的容器，然后就进入SIFT特征计算了。接下来我们进入到vl_sift_process_first_octave函数（主要用于处理第一个八度）。

vl_sift_process_first_octave (VlSiftFilt *f, vl_sift_pix const *im)
{
  int o, s, h, w ;
  double sa, sb ;
  vl_sift_pix *octave ;

  /* shortcuts */
  vl_sift_pix *temp   = f-> temp ;
  int width           = f-> width ;
  int height          = f-> height ;
  int o_min           = f-> o_min ;
  int s_min           = f-> s_min ;
  int s_max           = f-> s_max ;
  double sigma0       = f-> sigma0 ;
  double sigmak       = f-> sigmak ;
  double sigman       = f-> sigman ;
  double dsigma0      = f-> dsigma0 ;

  /* restart from the first */
  f->o_cur = o_min ;
  f->nkeys = 0 ;
  w = f-> octave_width  = VL_SHIFT_LEFT(f->width,  - f->o_cur) ;
  h = f-> octave_height = VL_SHIFT_LEFT(f->height, - f->o_cur) ;

  /* is there at least one octave? */
  if (f->O == 0)
    return VL_ERR_EOF ;

        首先这个函数有两个输入参数，VlSiftFilt *f 是指向SIFT滤波器结构体的指针，该结构体包含了SIFT计算所需的所有参数和状态信息。vl_sift_pix const *im 则是指向输入图像的指针。接下来定义了许多局部变量用于循环、计算尺度和处理图像大小。

• temp：指向SIFT滤波器中临时存储空间的指针。
• width和height：输入图像的宽度和高度。
• o_min、s_min、s_max：分别表示最小尺度组（octave）索引、最小尺度层级（scale level）索引和最大尺度层级索引。
• sigma0、sigmak、sigman、dsigma0：与高斯模糊相关的参数，用于控制尺度空间的构建。

        然后是重启操作，主要就是将一些变量归归位， f->o_cur = o_min 表示将当前尺度组索引f->o_cur重置为最小尺度组索引o_min。重置已检测到的关键点数量f->nkeys为0。根据当前尺度组索引和输入图像大小，计算当前尺度组的宽度w和高度h。VL_SHIFT_LEFT这个宏就不细说了，有兴趣的可以看看，主要就是按位左移实现快速的缩放计算。最后，如果判断没有足够的数据来构建至少一个完整的尺度组，函数就会返回VL_ERR_EOF错误码。

  octave = vl_sift_get_octave (f, s_min) ;

  if (o_min < 0) {
    /* double once */
    copy_and_upsample_rows (temp,   im,   width,      height) ;
    copy_and_upsample_rows (octave, temp, height, 2 * width ) ;

    /* double more */
    for(o = -1 ; o > o_min ; --o) {
      copy_and_upsample_rows (temp, octave,
                              width << -o,      height << -o ) ;
      copy_and_upsample_rows (octave, temp,
                              width << -o, 2 * (height << -o)) ;
    }
  }
  else if (o_min > 0) {
    /* downsample */
    copy_and_downsample (octave, im, width, height, o_min) ;
  }
  else {
    /* direct copy */
    memcpy(octave, im, sizeof(vl_sift_pix) * width * height) ;
  }

         在这段代码中，首先调用vl_sift_get_octave，获取第一组尺度空间的图像数据。vl_sift_get_octave代码如下：

vl_sift_get_octave (VlSiftFilt const *f, int s)
{
  int w = vl_sift_get_octave_width  (f) ;
  int h = vl_sift_get_octave_height (f) ;
  return f->octave + w * h * (s - f->s_min) ;
}

        很简单，获取长和宽后，主要计算f->octave + w * h * (s - f->s_min)，这一步其实是个偏移操作，找到当前尺度层级图像像素的起始位置，举个例子octave_data[0...63]是第一组第一层的图像，octave_data[64,127]是第一组第二层的图像，知道像素是怎么保存的这就很好理解。s - f->s_min是获取当前是第几个层级，因为层级不一定是从1开始，后面会说到。

         接下来回到vl_sift_process_first_octave中，根据o_min的值不同会有三种处理方法，如果o_min < 0，表示需要对输入图像进行上采样（放大）以构建第一组尺度空间。这是因为当o_min为负数时，第一组尺度空间的尺寸会大于原始图像尺寸。如果o_min < 0，表示需要对输入图像进行上采样（放大）以构建第一组尺度空间。这是因为当o_min为负数时，第一组尺度空间的尺寸会大于原始图像尺寸。如果o_min == 0，表示第一组尺度空间的尺寸与原始图像相同，因此直接将输入图像im复制到octave中。

        一般我们都是把原图像作为第一组第一层，也就是o_min = 0这种情况。

  sa = sigma0 * pow (sigmak,   s_min) ;
  sb = sigman * pow (2.0,    - o_min) ;

  if (sa > sb) {
    double sd = sqrt (sa*sa - sb*sb) ;
    _vl_sift_smooth (f, octave, temp, octave, w, h, sd) ;
  }

  /* -----------------------------------------------------------------
   *                                          Compute the first octave
   * -------------------------------------------------------------- */

  for(s = s_min + 1 ; s <= s_max ; ++s) {
    double sd = dsigma0 * pow (sigmak, s) ;
    _vl_sift_smooth (f, vl_sift_get_octave(f, s), temp,
                     vl_sift_get_octave(f, s - 1), w, h, sd) ;
  }

  return VL_ERR_OK ;

         这段代码根据一些初始参数（如sigma0、sigmak、sigman、o_min、s_min、s_max）计算出初始的平滑标准差（sa和sb），然后根据这些参数和计算出的标准差来平滑图像，以生成尺度空间中的各个尺度层级。

• sa 是基于 sigma0（初始标准差）和 sigmak（尺度因子）计算出的，用于第一组尺度空间（octave）中第一个尺度层级的平滑。
• sb 是基于 sigman（与尺度组大小相关的标准差）和 o_min（最小尺度组索引）计算出的，用于确定是否需要对原始图像进行额外的平滑以匹配第一组尺度空间的尺度。

        接下来判断如果 sa > sb，说明需要对原始图像进行额外的平滑，以匹配第一组尺度空间的第一个尺度层级的尺度。这里使用 _vl_sift_smooth 函数进行平滑，其中 sd 是根据 sa 和 sb 计算出的平滑标准差。

        然后就是计算第一个组的所有层级，代码遍历从 s_min + 1 到 s_max 的所有尺度层级索引 s，对于每个 s，计算该尺度层级的平滑标准差 sd。然后，使用 _vl_sift_smooth 函数对前一尺度层级的图像进行平滑，以生成当前尺度层级的图像。

        到此为止就构建出了图像金字塔的第一个组。那接下来我们会回到DescribeSIFT中。

    // Build alias to cached data
    auto regions = std::unique_ptr<Regions_type>(new Regions_type);

        这段代码中创建了regions,regions是一个std::unique_ptr<Regions_type>类型的智能指针，它用于管理和存储与SIFT（尺度不变特征变换）算法检测到的关键点（keypoints）和描述符（descriptors）相关的信息。

    // reserve some memory for faster keypoint saving
    regions->Features().reserve(2000);
    regions->Descriptors().reserve(2000);

        这两行代码调用了regions指向的Regions_type对象的成员函数，并使用reserve方法为即将存储的关键点和描述符预留了足够的内存空间。这样做可以提高后续添加元素的效率，因为容器不需要在每次添加新元素时都重新分配内存。这里假设大约会有2000个关键点和相应的描述符被检测并存储。

while (true) {
  vl_sift_detect(filt);

  VlSiftKeypoint const *keys  = vl_sift_get_keypoints(filt);
  const int nkeys = vl_sift_get_nkeypoints(filt);

  // Update gradient before launching parallel extraction
  vl_sift_update_gradient(filt);

  #ifdef OPENMVG_USE_OPENMP
  #pragma omp parallel for private(descr, descriptor)
  #endif
  for (int i = 0; i < nkeys; ++i) {

    // Feature masking
    if (mask)
    {
      const image::Image<unsigned char> & maskIma = *mask;
      if (maskIma(keys[i].y, keys[i].x) == 0)
        continue;
    }

        接下来进入到while true循环，首先调用SIFT检测函数来在当前图像或视频帧中检测关键点，通过vl_sift_get_keypoints(filt)和vl_sift_get_nkeypoints(filt)获取检测到的关键点和它们的数量。关键点存储在VlSiftKeypoint类型的数组中，数量存储在nkeys变量中。vl_sift_update_gradient(filt); 在进行并行描述符提取之前更新图像的梯度。这是因为SIFT描述符的计算依赖于图像的梯度信息。

        如果定义了OPENMVG_USE_OPENMP宏（这通常意味着代码可以在支持OpenMP的编译器上并行编译和执行），则使用#pragma omp parallel for指令来并行化关键点的处理。这可以提高处理大量关键点时的效率。

        然后遍历每个关键点，如果提供了掩码（mask非空），则检查当前关键点是否位于掩码图像的未掩蔽区域。如果关键点位于掩蔽区域（即掩码图像在该点的值为0），则跳过该关键点的后续处理。

      double angles [4] = {0.0, 0.0, 0.0, 0.0};
      int nangles = 1; // by default (1 upright feature)
      if (_bOrientation)
      { // compute from 1 to 4 orientations
        nangles = vl_sift_calc_keypoint_orientations(filt, angles, keys+i);
      }

      for (int q=0 ; q < nangles ; ++q) {
        vl_sift_calc_keypoint_descriptor(filt, &descr[0], keys+i, angles[q]);
        const SIOPointFeature fp(keys[i].x, keys[i].y,
          keys[i].sigma, static_cast<float>(angles[q]));

        siftDescToUChar(&descr[0], descriptor, _params._root_sift);
        #ifdef OPENMVG_USE_OPENMP
        #pragma omp critical
        #endif
        {
          regions->Descriptors().push_back(descriptor);
          regions->Features().push_back(fp);
        }
      }
    }
    if (vl_sift_process_next_octave(filt))
      break; // Last octave
  }
  vl_sift_delete(filt);

  return regions;
}

        如果启用了方向计算（_bOrientation为真），则调用vl_sift_calc_keypoint_orientations函数来计算关键点的一个或多个方向。这个函数会填充angles数组，并返回计算出的方向数量（nangles）。方向计算是可选的，因为某些情况下只需要关键点而不需要它们的方向。这里可能计算出四个方向，主要是为了提高鲁棒性。然后，通过一个for循环遍历每个角度/方向（nangles表示角度的数量），这通常用于计算每个关键点的多个方向（通常是主方向及其周围的一些次要方向）的描述符。对于每个角度q，调用vl_sift_calc_keypoint_descriptor函数计算关键点的描述符。这个函数需要几个参数：滤波后的图像filt、用于存储描述符的数组descr（这里从&descr[0]开始），当前关键点的信息keys+i（i是一个外部循环的索引，用于遍历关键点），以及当前的角度angles[q]。然后，根据当前关键点的位置x、y、尺度sigma和当前的角度angles[q]，创建一个SIOPointFeature对象fp。使用siftDescToUChar函数将计算出的描述符（存储在descr中）转换为另一种格式（是为了存储或传输的便利性），并存储在descriptor变量中。在每次迭代结束时，检查是否还有下一个组（通过调用vl_sift_process_next_octave(filt)），如果有，则继续循环；否则，跳出循环（一般是设置6个组）最后，使用vl_sift_delete函数释放与SIFT滤波器相关的资源。然后继续遍历下一个方向，如果为最后一个则遍历下一个关键点......

        整个流程说完，但是最重要的几个函数还没有讲，vl_sift_detect，vl_sift_get_keypoints，vl_sift_update_gradient，vl_sift_process_next_octave这几个关键的函数，我会在接下来几个博客中详细讲解，谢谢大家！