SIFT源码解析（openMVG-2.1）

本文对openMVG 2.1中的sift相关源码进行了解析，算法原理部分请见上一篇文章万字长文详解SIFT特征提取，本文代码中涉及的一些细节问题可以当作上一篇文章的补充。
调用例程

openMVG的项目源代码中有sift的调用例程，见：
openMVG-2.1\src\openMVG_Samples\features_siftPutativeMatches\siftmatch.cpp
其中核心调用sift的代码：

1. using namespace openMVG::features;
2. std::unique_ptr<Image_describer> image_describer(new SIFT_Anatomy_Image_describer);  //核心类初始化，具体分析见下文
3. std::map<IndexT, std::unique_ptr<features::Regions>> regions_perImage; //用于储存特征
4. image_describer->Describe(imageL, regions_perImage[0]); //主要执行函数
5. image_describer->Describe(imageR, regions_perImage[1]);

复制代码

SIFT_Anatomy_Image_describer类

算法的核心类为 SIFT_Anatomy_Image_describer，位于SIFT_Anatomy_Image_Describer.hpp中
核心参数设置

其参数有一个struct传入，在初始化时会有一个默认的参数值初始化，如下：

1.     Params(
2.       int first_octave = 0, //如果设置为-1，则第一个octave中使用的是扩大一倍的图像，也就是upscale一次
3.       int num_octaves = 6, //最大的octave数量
4.       int num_scales = 3,     // 每层octave输出的尺度数量
5.       float edge_threshold = 10.0f, // 去除边缘时用到的阈值
6.       float peak_threshold = 0.04f,  // 去除对比度小的点时用到的阈值
7.       bool root_sift = true  // 是否使用改进的sift描述符

复制代码

first_octave，在openMVG中，默认的首层octave图像是源图像大小，可以通过这个参数设置为扩大一倍的大小。这里upscale主要是为了能够检测到更加细节的特征，此外，如果原图的分辨率本身较低，upscale一下也有助于检测到小尺度特征。
另，在opencv中，是直接默认首层upscale的。
num_octaves，设置的最大octave数量，此外，在后面的代码中还会根据图片尺寸算一个octave数量，来保证最后一层的分辨率不要太小（不小于32X32），在HierarchicalGaussianScaleSpace类中我们详细介绍。
num_scales每个octave的输出scale数量，这个数量加3，就是每个octave所需要的高斯层数量。
edge_threshold这里判定是否为边缘点的方法和之前原理篇讲的略有不同，但原理基本一致，都是得到局部Hessian矩阵后，判断特征值的大小情况。
在之前的判定中，我们计算出Hessian矩阵M后，采用这个公式计算其R响应：

\[R = \text{Det}(M) - \alpha \text{Tr}(M)^2=\lambda_{1}\lambda_{2}-\alpha(\lambda_{1}+\lambda_{2})^2\]
其中α是一个经验值 0.04 到 0.06。
但在openMVG中，我们假设 $$λ_1 = γλ_2$$
其中γ > 1，该值越大，表示两者之间的比值越大，越有可能是边缘。

我们采用trace平方和determinant的比值来判定：

\[\frac{\mathrm{Tr}(\boldsymbol{M})^2}{\mathrm{Det}(\boldsymbol{M})}=\frac{(\lambda_1 + \lambda_2)^2}{\lambda_1\lambda_2}=\frac{(\gamma\lambda_2+\lambda_2)^2}{\gamma\lambda_2^2}=\frac{(\gamma + 1)^2}{\gamma}\]
这里如果得到的值越大，表示γ的值越大，也就表示越有可能是边缘点。edge_threshold就是这里设置的γ，SIFT原作者给出一个经验值10，在openMVG中也是用的这个经验值。
peak_threshold这个用来去掉对比度小的点，计算方法我们在SIFT_KeypointExtractor类中详述。
root_sift改进的sift描述符，详见Sift_DescriptorExtractor类。
算法流程核心 Describe_SIFT_Anatomy

包含整体的特征提取算法流程。核心流程函数！

1. std::unique_ptr<Regions_type> Describe_SIFT_Anatomy(
2.   const image::Image<unsigned char>& image,   
3.   const image::Image<unsigned char>* mask = nullptr  
4. )
5. //image：传入图像
6. //mask：掩码图像，用于过滤关键点。非零值表示感兴趣区域。通常传一个空即可
7. {
8.   auto regions = std::unique_ptr<Regions_type>(new Regions_type);
10.   if (image.size() == 0)
11.     return regions;
13.   // Convert to float in range [0;1]
14.   //将输入图像转换为浮点数类型，并将像素值归一化到 [0, 1] 范围
15.   const image::Image<float> If(image.GetMat().cast<float>()/255.0f);
17.   // compute sift keypoints
18.   // 计算 SIFT 关键点
19.   {
20.     const int supplementary_images = 3;
21.     // 高斯层数和输出尺度之间的数量差3，详见sift原理文章。
22.     // => in order to ensure each gaussian slice is used in the process 3 extra images are required:
23.     // +1 for dog computation
24.     // +2 for 3d discrete extrema definition
26.     //生成高斯尺度空间，详见后文具体的类分析
27.     HierarchicalGaussianScaleSpace octave_gen(
28.       params_.num_octaves_,
29.       params_.num_scales_,
30.       (params_.first_octave_ == -1) // 这里判定首个octave是否要upscale，如果要upscale，相关尺度要进行调整
31.       ? GaussianScaleSpaceParams(1.6f/2.0f, 1.0f/2.0f, 0.5f, supplementary_images)
32.       : GaussianScaleSpaceParams(1.6f, 1.0f, 0.5f, supplementary_images));
33.     octave_gen.SetImage( If );  //设置输入图像，详见下文相关类分析
35.     std::vector<sift::Keypoint> keypoints; // 存储所有检测到的关键点
36.     keypoints.reserve(5000);  // 预分配空间，提高性能
37.     Octave octave; // 存储当前octave的信息，以后每次迭代都是修改这个里面的信息
38.     while ( octave_gen.NextOctave( octave ) )  // 遍历octave
39.     {
40.       std::vector<sift::Keypoint> keys;
41.       // Find Keypoints 创建一个 SIFT 关键点检测器
42.       sift::SIFT_KeypointExtractor keypointDetector(
43.         params_.peak_threshold_ / octave_gen.NbSlice(),  
44.         params_.edge_threshold_);
45.     //参数1：这里就是原理篇中讲的对比度阈值计算方法，其中peak_threshold_为经验值
46.     //参数2：边缘阈值，详见上文
47.       keypointDetector(octave, keys); // 检测关键点，详细过程见下文
48.       // 创建一个 SIFT 描述符提取器
49.       sift::Sift_DescriptorExtractor descriptorExtractor;
50.       descriptorExtractor(octave, keys); // 计算当前octave中关键点的方向和描述符
52.       // Concatenate the found keypoints
53.       std::move(keys.begin(), keys.end(), std::back_inserter(keypoints));
54.     }
55.     for (const auto & k : keypoints)
56.     {
57.       // Feature masking 特征掩码处理
58.       if (mask)
59.       {
60.         const image::Image<unsigned char> & maskIma = *mask;
61.         if (maskIma(k.y, k.x) == 0) // 如果关键点位于掩码的非感兴趣区域，则跳过该关键点
62.           continue;
63.       }
64.       // Create the SIFT region 创建 SIFT 区域
65.       {
66.         // 将关键点的描述符转换为无符号字符类型，并添加到 Regions 对象中
67.           regions->Descriptors().emplace_back(k.descr.cast<unsigned char>());
68.           // 将关键点的坐标、尺度和方向添加到 Regions 对象中
69.           regions->Features().emplace_back(k.x, k.y, k.sigma, k.theta);
70.       }
71.     }
72.   }
73.   return regions;
74. };

复制代码

在构建高斯金字塔时，往GaussianScaleSpaceParams(1.6f/2.0f, 1.0f/2.0f, 0.5f, supplementary_images)中传了一个0.5，我们默认输入图像已经进行了一次高斯卷积操作，其σ为这个0.5，那么第一层的图像就是从这个0.5卷积后的图像基础上，再进行一次卷积得到的。如果第一层所需的σ为1.6，则再次卷积所需的高斯σ为：

\[\sigma_{extra} = \sqrt{1.6^2-0.5^2}\]
而不是直接对图像做σ为1.6的卷积。
HierarchicalGaussianScaleSpace类

构建高斯金字塔的核心类
设置初始图象 SetImage

设置初始图像，主要是判断是否要upscale，以及计算octave的数量。

1.   virtual void SetImage(const image::Image<float> & img)
2.   {
3.     const double sigma_extra =
4.       sqrt(Square(m_params.sigma_min) - Square(m_params.sigma_in)) / m_params.delta_min; // 计算额外的高斯模糊标准差，用于在图像上应用额外的高斯滤波，见上文分析
5.     if (m_params.delta_min == 1.0f) // 如果最小采样距离为 1.0f，直接对输入图像应用高斯滤波
6.     {
7.       image::ImageGaussianFilter(img, sigma_extra, m_cur_base_octave_image);
8.     }
9.     else  // 如果最小采样距离为 0.5f，需要对图像进行上采样
10.     {
11.       if (m_params.delta_min == 0.5f)
12.       {
13.         image::Image<float> tmp;
14.         ImageUpsample(img, tmp); //上采样
15.         image::ImageGaussianFilter(tmp, sigma_extra, m_cur_base_octave_image);
16.       }
17.       else
18.       {
19.         OPENMVG_LOG_ERROR
20.           << "Upsampling or downsampling with delta equal to: "
21.           << m_params.delta_min << " is not yet implemented";
22.       }
23.     }
24.     //-- Limit the size of the last octave to be at least 32x32 pixels
25.     // 保证最后一个octave至少为 32x32 像素，这个要和传入的最大octave数量做个对比，取小的那个作为octave数量的确定值
26.     const int nbOctaveMax = std::ceil(std::log2( std::min(m_cur_base_octave_image.Width(), m_cur_base_octave_image.Height())/32));
27.     m_nb_octave = std::min(m_nb_octave, nbOctaveMax);
28.   }

复制代码

const int index = (m_params.supplementary_levels == 0) ? 1 : m_params.supplementary_levels; 这里我们用第一个octave举例，如果supplementary_levels为3，则index = 3，那么下采样的层就是第6-3=3层，也就是sigma为\(2σ_{initial}\)的层，把这层作为第二个octave的基础图象。也就是说本来第二个octave的基础图象应该是\(2σ_{initial}\)对原图做卷积，这在第一个octave中已经计算过了，可以直接拿到第二个octave中去用。
SIFT_KeypointExtractor类

特征点提取类。在每一个octave迭代时，包含两个关键步骤，极值提取、极值修正。
默认参数

1.   virtual bool NextOctave(Octave & octave)
2.   {
3.     if (m_cur_octave_id >= m_nb_octave)
4.     {
5.       return false;  // 如果达到最大数量，返回 false 表示没有更多的octave可供计算
6.     }
7.     else
8.     {
9.       octave.octave_level = m_cur_octave_id;
10.       if (m_cur_octave_id == 0) // 如果是第一个octave
11.       {
12.         octave.delta = m_params.delta_min; //设置采样率为最小采样率
13.       }
14.       else
15.       {
16.         octave.delta *= 2.0f; // 对于后续的octave，采样率翻倍（相当于图像缩小）
17.       }
19.       // init the "blur"/sigma scale spaces values
20.       octave.slices.resize(m_nb_slice + m_params.supplementary_levels);
21.       octave.sigmas.resize(m_nb_slice + m_params.supplementary_levels);
22.       for (int s = 0; s < m_nb_slice  + m_params.supplementary_levels; ++s)
23.       {
24.         octave.sigmas[s] =
25.           octave.delta / m_params.delta_min * m_params.sigma_min * pow(2.0,(float)s/(float)m_nb_slice);   //这里计算octave内每一层所需的sigma，其中k = 2^{1/s}，详见原理篇文章
26.       }
28.       // Build the octave iteratively
29.       octave.slices[0] = m_cur_base_octave_image;
30.       for (int s = 1; s < octave.sigmas.size(); ++s)
31.       {
32.         // Iterative blurring the previous image
33.         const image::Image<float> & im_prev = octave.slices[s-1];
34.         image::Image<float> & im_next = octave.slices[s];
35.         const double sig_prev = octave.sigmas[s-1];
36.         const double sig_next = octave.sigmas[s];
37.         const double sigma_extra = sqrt(Square(sig_next) - Square(sig_prev)) / octave.delta;  // 下一层的sigma是由当前层再进行一次卷积得到的，这个额外的卷积核计算方法遵循勾股定理
39.         image::ImageGaussianFilter(im_prev, sigma_extra, im_next); //高斯模糊生成该层图像，每一层都是由上一层进行卷积的到，而非由原始图像卷积得到
40.       }
42.       // Prepare for next octave computation -> Decimation
43.       ++m_cur_octave_id;
44.       if (m_cur_octave_id < m_nb_octave)
45.       {
46.         // Decimate => sigma * 2 for the next iteration
47.         const int index = (m_params.supplementary_levels == 0) ? 1 : m_params.supplementary_levels;
48.         // 对选定的切片进行下采样，作为下一个octave的基础图像
49.         ImageDecimate(octave.slices[octave.sigmas.size()-index], m_cur_base_octave_image);
50.       }
51.       return true;
52.     }
53.   }

复制代码

核心流程

1.     float peak_threshold,     // 对比度阈值，详见上文
2.     float edge_threshold,     // 边缘去除阈值，详见上文
3.     int nb_refinement_step = 5 // 极值修正是一个迭代过程，这里限制了最大迭代次数

复制代码

极值提取 Find_3d_discrete_extrema

在 DoG 空间中查找 3D 离散极值点

1.   // 重载函数调用运算符，计算给定高斯octave的差分高斯（DoG），并在 DoG 空间中查找离散极值点
2.   //对这些关键点的位置进行细化，以得到更精确的关键点位置。
3.   void operator()( const Octave & octave , std::vector<Keypoint> & keypoints )
4.   {
5.     if (!ComputeDogs(octave)) // 差分高斯（DoG），就是两个图像相减，不再赘述
6.       return;
7.     Find_3d_discrete_extrema(keypoints, 0.8f);   //找到离散极值点，这里的0.8是对对比度阈值的一次缩放，可以初步筛出更多的极值点；默认值为1，对原阈值不做修改
8.     Keypoints_refine_position(keypoints);        //对极值点进行修正
9.   }

复制代码

极值修正 Keypoints_refine_position

细化关键点的位置（空间位置和尺度），丢弃无法细化的关键点
[code]//遍历输入的关键点列表，对每个关键点进行位置细化。//通过二次函数插值来优化关键点的位置，直到满足收敛条件或达到最大迭代次数。//对细化后的关键点进行峰值阈值、边缘响应和边界检查，只有通过所有检查的关键点才会被保留。void Keypoints_refine_position(  std::vector & keypoints) const{  std::vector kps;  kps.reserve(keypoints.size());  const float ofstMax = 0.6f; // 细化参数的最大偏移量  // Ratio between two consecutive scales in the slice  // 这个就是k  const float sigma_ratio = m_Dogs.sigmas[1] / m_Dogs.sigmas[0];  // 边缘响应阈值, 计算原理详见上文  const float edge_thres = Square(m_edge_threshold + 1) / m_edge_threshold;  const Octave & octave = m_Dogs;  const int w = octave.slices[0].Width();  const int h = octave.slices[0].Height();  const float delta  = octave.delta;  for (const auto & key : keypoints) // 遍历keypoints  {    float val = key.val;    int ic = key.i; // current discrete value of x coordinate - at each interpolation    int jc = key.j; // current discrete value of y coordinate - at each interpolation    int sc = key.s; // current discrete value of s coordinate - at each interpolation    float ofstX = 0.0f;    float ofstY = 0.0f;    float ofstS = 0.0f;    bool isConv = false;    // While position cannot be refined and the number of refinement is not too much    // 迭代细化极值点    for ( int nIntrp = 0; !isConv && nIntrp < m_nb_refinement_step; ++nIntrp)    {      // Extrema interpolation via a quadratic function      //   only if the detection is not too close to the border (so the discrete 3D Hessian is well defined)      if ( 0 < ic &&  ic < (w-1) && 0 < jc && jc < (h-1) )      {        // 通过taylor展开得到极值点附近的二次项近似，求取极值点偏移，求解过程实际是解一个导数为0的方程，详见算法原理篇        if (Inverse_3D_Taylor_second_order_expansion(octave, ic, jc, sc, &ofstX, &ofstY, &ofstS, &val, ofstMax))          isConv = true;      }      else      {        isConv = false;      }      if (!isConv)      {        // let's explore the neighborhood in        // space...        if (ofstX > +ofstMax && (ic+1) < (w-1) ) {++ic;}        if (ofstX < -ofstMax && (ic-1) >  0    ) {--ic;}        if (ofstY > +ofstMax && (jc+1) < (h-1) ) {++jc;}        if (ofstY < -ofstMax && (jc-1) >  0    ) {--jc;}        // ... and scale.        /*        if (ofstS > +ofstMax && (sc+1) < (ns-1)) {++sc;}        if (ofstS < -ofstMax && (sc-1) >    0  ) {--sc;}        */      }    }    if (isConv)    {      // Peak threshold check      if ( std::abs(val) > m_peak_threshold )  // 对比度高的点保留      {        Keypoint kp = key;        kp.x = (ic + ofstX) * delta;        kp.y = (jc + ofstY) * delta;        kp.i = ic;        kp.j = jc;        kp.s = sc;        kp.sigma = octave.sigmas[sc] * pow(sigma_ratio, ofstS); // logarithmic scale        kp.val = val;        // Edge check 高斯差分对边缘的响应也很高，要把这些边缘点去除掉        if (Compute_edge_response(kp) >=0 && std::abs(kp.edgeResp)