SIFT源码解析（openMVG-2.1）
本文对openMVG 2.1中的sift相关源码进行了解析，算法原理部分请见上一篇文章万字长文详解SIFT特征提取，本文代码中涉及的一些细节问题可以当作上一篇文章的补充。调用例程
openMVG的项目源代码中有sift的调用例程，见：
openMVG-2.1\src\openMVG_Samples\features_siftPutativeMatches\siftmatch.cpp
其中核心调用sift的代码：
using namespace openMVG::features;
std::unique_ptr<Image_describer> image_describer(new SIFT_Anatomy_Image_describer);//核心类初始化，具体分析见下文
std::map<IndexT, std::unique_ptr<features::Regions>> regions_perImage; //用于储存特征
image_describer->Describe(imageL, regions_perImage); //主要执行函数
image_describer->Describe(imageR, regions_perImage);SIFT_Anatomy_Image_describer类
算法的核心类为 SIFT_Anatomy_Image_describer，位于SIFT_Anatomy_Image_Describer.hpp中
核心参数设置
其参数有一个struct传入，在初始化时会有一个默认的参数值初始化，如下：
Params(
int first_octave = 0, //如果设置为-1，则第一个octave中使用的是扩大一倍的图像，也就是upscale一次
int num_octaves = 6, //最大的octave数量
int num_scales = 3, // 每层octave输出的尺度数量
float edge_threshold = 10.0f, // 去除边缘时用到的阈值
float peak_threshold = 0.04f,// 去除对比度小的点时用到的阈值
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响应：
\
其中α是一个经验值 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
包含整体的特征提取算法流程。核心流程函数！
std::unique_ptr<Regions_type> Describe_SIFT_Anatomy(
const image::Image<unsigned char>& image,
const image::Image<unsigned char>* mask = nullptr
)
//image：传入图像
//mask：掩码图像，用于过滤关键点。非零值表示感兴趣区域。通常传一个空即可
{
auto regions = std::unique_ptr<Regions_type>(new Regions_type);
if (image.size() == 0)
return regions;
// Convert to float in range
//将输入图像转换为浮点数类型，并将像素值归一化到 范围
const image::Image<float> If(image.GetMat().cast<float>()/255.0f);
// compute sift keypoints
// 计算 SIFT 关键点
{
const int supplementary_images = 3;
// 高斯层数和输出尺度之间的数量差3，详见sift原理文章。
// => in order to ensure each gaussian slice is used in the process 3 extra images are required:
// +1 for dog computation
// +2 for 3d discrete extrema definition
//生成高斯尺度空间，详见后文具体的类分析
HierarchicalGaussianScaleSpace octave_gen(
params_.num_octaves_,
params_.num_scales_,
(params_.first_octave_ == -1) // 这里判定首个octave是否要upscale，如果要upscale，相关尺度要进行调整
? GaussianScaleSpaceParams(1.6f/2.0f, 1.0f/2.0f, 0.5f, supplementary_images)
: GaussianScaleSpaceParams(1.6f, 1.0f, 0.5f, supplementary_images));
octave_gen.SetImage( If );//设置输入图像，详见下文相关类分析
std::vector<sift::Keypoint> keypoints; // 存储所有检测到的关键点
keypoints.reserve(5000);// 预分配空间，提高性能
Octave octave; // 存储当前octave的信息，以后每次迭代都是修改这个里面的信息
while ( octave_gen.NextOctave( octave ) )// 遍历octave
{
std::vector<sift::Keypoint> keys;
// Find Keypoints 创建一个 SIFT 关键点检测器
sift::SIFT_KeypointExtractor keypointDetector(
params_.peak_threshold_ / octave_gen.NbSlice(),
params_.edge_threshold_);
//参数1：这里就是原理篇中讲的对比度阈值计算方法，其中peak_threshold_为经验值
//参数2：边缘阈值，详见上文
keypointDetector(octave, keys); // 检测关键点，详细过程见下文
// 创建一个 SIFT 描述符提取器
sift::Sift_DescriptorExtractor descriptorExtractor;
descriptorExtractor(octave, keys); // 计算当前octave中关键点的方向和描述符
// Concatenate the found keypoints
std::move(keys.begin(), keys.end(), std::back_inserter(keypoints));
}
for (const auto & k : keypoints)
{
// Feature masking 特征掩码处理
if (mask)
{
const image::Image<unsigned char> & maskIma = *mask;
if (maskIma(k.y, k.x) == 0) // 如果关键点位于掩码的非感兴趣区域，则跳过该关键点
continue;
}
// Create the SIFT region 创建 SIFT 区域
{
// 将关键点的描述符转换为无符号字符类型，并添加到 Regions 对象中
regions->Descriptors().emplace_back(k.descr.cast<unsigned char>());
// 将关键点的坐标、尺度和方向添加到 Regions 对象中
regions->Features().emplace_back(k.x, k.y, k.sigma, k.theta);
}
}
}
return regions;
};在构建高斯金字塔时，往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的数量。
virtual void SetImage(const image::Image<float> & img)
{
const double sigma_extra =
sqrt(Square(m_params.sigma_min) - Square(m_params.sigma_in)) / m_params.delta_min; // 计算额外的高斯模糊标准差，用于在图像上应用额外的高斯滤波，见上文分析
if (m_params.delta_min == 1.0f) // 如果最小采样距离为 1.0f，直接对输入图像应用高斯滤波
{
image::ImageGaussianFilter(img, sigma_extra, m_cur_base_octave_image);
}
else// 如果最小采样距离为 0.5f，需要对图像进行上采样
{
if (m_params.delta_min == 0.5f)
{
image::Image<float> tmp;
ImageUpsample(img, tmp); //上采样
image::ImageGaussianFilter(tmp, sigma_extra, m_cur_base_octave_image);
}
else
{
OPENMVG_LOG_ERROR
<< "Upsampling or downsampling with delta equal to: "
<< m_params.delta_min << " is not yet implemented";
}
}
//-- Limit the size of the last octave to be at least 32x32 pixels
// 保证最后一个octave至少为 32x32 像素，这个要和传入的最大octave数量做个对比，取小的那个作为octave数量的确定值
const int nbOctaveMax = std::ceil(std::log2( std::min(m_cur_base_octave_image.Width(), m_cur_base_octave_image.Height())/32));
m_nb_octave = std::min(m_nb_octave, nbOctaveMax);
}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迭代时，包含两个关键步骤，极值提取、极值修正。
默认参数
virtual bool NextOctave(Octave & octave)
{
if (m_cur_octave_id >= m_nb_octave)
{
return false;// 如果达到最大数量，返回 false 表示没有更多的octave可供计算
}
else
{
octave.octave_level = m_cur_octave_id;
if (m_cur_octave_id == 0) // 如果是第一个octave
{
octave.delta = m_params.delta_min; //设置采样率为最小采样率
}
else
{
octave.delta *= 2.0f; // 对于后续的octave，采样率翻倍（相当于图像缩小）
}
// init the "blur"/sigma scale spaces values
octave.slices.resize(m_nb_slice + m_params.supplementary_levels);
octave.sigmas.resize(m_nb_slice + m_params.supplementary_levels);
for (int s = 0; s < m_nb_slice+ m_params.supplementary_levels; ++s)
{
octave.sigmas =
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}，详见原理篇文章
}
// Build the octave iteratively
octave.slices = m_cur_base_octave_image;
for (int s = 1; s < octave.sigmas.size(); ++s)
{
// Iterative blurring the previous image
const image::Image<float> & im_prev = octave.slices;
image::Image<float> & im_next = octave.slices;
const double sig_prev = octave.sigmas;
const double sig_next = octave.sigmas;
const double sigma_extra = sqrt(Square(sig_next) - Square(sig_prev)) / octave.delta;// 下一层的sigma是由当前层再进行一次卷积得到的，这个额外的卷积核计算方法遵循勾股定理
image::ImageGaussianFilter(im_prev, sigma_extra, im_next); //高斯模糊生成该层图像，每一层都是由上一层进行卷积的到，而非由原始图像卷积得到
}
// Prepare for next octave computation -> Decimation
++m_cur_octave_id;
if (m_cur_octave_id < m_nb_octave)
{
// Decimate => sigma * 2 for the next iteration
const int index = (m_params.supplementary_levels == 0) ? 1 : m_params.supplementary_levels;
// 对选定的切片进行下采样，作为下一个octave的基础图像
ImageDecimate(octave.slices, m_cur_base_octave_image);
}
return true;
}
}核心流程
float peak_threshold, // 对比度阈值，详见上文
float edge_threshold, // 边缘去除阈值，详见上文
int nb_refinement_step = 5 // 极值修正是一个迭代过程，这里限制了最大迭代次数极值提取 Find_3d_discrete_extrema
在 DoG 空间中查找 3D 离散极值点
// 重载函数调用运算符，计算给定高斯octave的差分高斯（DoG），并在 DoG 空间中查找离散极值点
//对这些关键点的位置进行细化，以得到更精确的关键点位置。
void operator()( const Octave & octave , std::vector<Keypoint> & keypoints )
{
if (!ComputeDogs(octave)) // 差分高斯（DoG），就是两个图像相减，不再赘述
return;
Find_3d_discrete_extrema(keypoints, 0.8f); //找到离散极值点，这里的0.8是对对比度阈值的一次缩放，可以初步筛出更多的极值点；默认值为1，对原阈值不做修改
Keypoints_refine_position(keypoints); //对极值点进行修正
}极值修正 Keypoints_refine_position
细化关键点的位置（空间位置和尺度），丢弃无法细化的关键点
//遍历输入的关键点列表，对每个关键点进行位置细化。//通过二次函数插值来优化关键点的位置，直到满足收敛条件或达到最大迭代次数。//对细化后的关键点进行峰值阈值、边缘响应和边界检查，只有通过所有检查的关键点才会被保留。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// 这个就是kconst float sigma_ratio = m_Dogs.sigmas / m_Dogs.sigmas;// 边缘响应阈值, 计算原理详见上文const float edge_thres = Square(m_edge_threshold + 1) / m_edge_threshold;const Octave & octave = m_Dogs;const int w = octave.slices.Width();const int h = octave.slices.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 * pow(sigma_ratio, ofstS); // logarithmic scale kp.val = val; // Edge check 高斯差分对边缘的响应也很高，要把这些边缘点去除掉 if (Compute_edge_response(kp) >=0 && std::abs(kp.edgeResp)
页:
[1]