之前结合不同人的资料理解了sift的原理，这里通过opencv中的代码来加深对sift的实现的理解。
使得能够从原理性理解到源码级的理解。不过该博文还是大量基于《赵春江, opencv2.4.9 源码分析，SIFT http://blog.youkuaiyun.com/zhaocj》的。
在opencv3.0中，已经看不到sift.cpp源代码了，在2.4.10中还是有的：opencv\sources\modules\nonfree\src下面。可以看出这是一个非免费代码。使用需要付费的，毕竟sift是哥伦比亚大学的专利。
对于sift（1）原理的对照来解释代码。
ps：在2.3与2.4版本中该sift的代码是两个不同的版本。
ps：2.4.10相比较于2.4.4中该源码的546处增加了两行，修复了一个小bug。
#define CV_EXPORTS __declspec(dllexport)
#define CV_EXPORTS_W CV_EXPORTS
#define CV_WRAP
这两个关键字在这里有人回答了。
2、sift类的定义：
在文件“opencv\sources\modules\features2d\include\opencv2\features2d.hpp”中：
class CV_EXPORTS_W_SIMPLE KeyPoint{}；
在文件“opencv\sources\modules\nonfree\include\opencv2\nonfree\features2d.hpp”中：
#ifndef __OPENCV_NONFREE_FEATURES_2D_HPP__
#define __OPENCV_NONFREE_FEATURES_2D_HPP__
#include "opencv2/features2d/features2d.hpp"
//该头文件中包含了类Feature2D的定义。而看过里面这个Feature2D又是基于两个基类派生出来的,这里就不接着往下挖了。
//class CV_EXPORTS_W Feature2D : public FeatureDetector,
// public DescriptorExtracto{}
#ifdef __cplusplus
namespace cv
{
class CV_EXPORTS_W SIFT : public Feature2D
{
public:
CV_WRAP explicit SIFT( int nfeatures=0, int nOctaveLayers=3,
double contrastThreshold=0.04, double edgeThreshold=10,
double sigma=1.6);
//!返回描述子的大小(128)
CV_WRAP int descriptorSize() const;
//! 返回描述子类型
CV_WRAP int descriptorType() const;
//! 使用SIFT算法找到关键点
void operator()(InputArray img, InputArray mask,
vector<KeyPoint>& keypoints) const;
//! 使用SIFT算法找到关键点然后计算描述子。
//! （可选的）使用该操作计算用户提供的关键字的描述子
void operator()(InputArray img, InputArray mask,
vector<KeyPoint>& keypoints,
OutputArray descriptors,
bool useProvidedKeypoints=false) const;
AlgorithmInfo* info() const;
//!建立高斯金字塔
void buildGaussianPyramid( const Mat& base, vector<Mat>& pyr, int nOctaves ) const;
//！建立DoG金字塔
void buildDoGPyramid( const vector<Mat>& pyr, vector<Mat>& dogpyr ) const;
//！查找尺度空间的极值
void findScaleSpaceExtrema( const vector<Mat>& gauss_pyr, const vector<Mat>& dog_pyr,
vector<KeyPoint>& keypoints ) const;
protected:
void detectImpl( const Mat& image, vector<KeyPoint>& keypoints, const Mat& mask=Mat() ) const;
void computeImpl( const Mat& image, vector<KeyPoint>& keypoints, Mat& descriptors ) const;
CV_PROP_RW int nfeatures;//特征
CV_PROP_RW int nOctaveLayers;//组数
CV_PROP_RW double contrastThreshold;
CV_PROP_RW double edgeThreshold;//边缘阈值
CV_PROP_RW double sigma;//尺度
};
typedef SIFT SiftFeatureDetector;
typedef SIFT SiftDescriptorExtractor;
} /* namespace cv */
#endif /* __cplusplus */
#endif
/* End of file. */
上面就是SIFT类的定义，从中可以看出该类主要实现的是高斯金字塔的建立、DoG金字塔的建立、尺度空间中极值点的搜寻等等。
0、SIFT构造函数，查询用户是否提供关键点，然后调用1、2、3、4、5；
1、首先我们要做的就是建立初始图像createInitialImage();
2、建立高斯金字塔 SIFT::buildGaussianPyramid();
3、建立DoG金字塔 SIFT::buildDoGPyramid();
4、在DoG中查找尺度空间极值点SIFT::findScaleSpaceExtrema();其中会调用4.1和4.2
4.1、对找到的极值点进行曲线插值拟合，并过滤 adjustLocalExtrema();
4.2、计算方向直方图 calcOrientationHist();
5、计算描述子 calcDescriptors();
5.1、计算sift描述子calcSIFTDescriptor();
1、
//按照用户指定( bool doubleImageSize)是否需要用线性插值的方法扩展图像；然后进行高斯模糊。
//公式a
static Mat createInitialImage( const Mat& img, bool doubleImageSize, float sigma )
{
Mat gray, gray_fpt;
//如果输入图像是彩色图像，则需要转换成灰度图像
if( img.channels() == 3 || img.channels() == 4 )
cvtColor(img, gray, COLOR_BGR2GRAY);
else
img.copyTo(gray);
//调整图像的像素数据类型
//typedef float sift_wt;
//static const int SIFT_FIXPT_SCALE = 1
gray.convertTo(gray_fpt, DataType<sift_wt>::type, SIFT_FIXPT_SCALE, 0);
float sig_diff;
//如果需要扩大图像的长宽尺寸
if( doubleImageSize )
{
// SIFT_INIT_SIGMA 为0.5，即输入图像的尺度，SIFT_INIT_SIGMA×2=1.0，即图像扩大2 倍以后的尺度，sig_diff 为公式6 中的高斯函数所需要的方差
sig_diff = sqrtf( std::max(sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA * 4, 0.01f) );
Mat dbl;
//利用双线性插值法把图像的长宽都扩大2 倍
resize(gray_fpt, dbl, Size(gray.cols*2, gray.rows*2), 0, 0, INTER_LINEAR);
//利用上面提到的公式a对图像进行高斯平滑处理
GaussianBlur(dbl, dbl, Size(), sig_diff, sig_diff);
//输出图像矩阵
return dbl;
}
else
{
//不需要扩大图像的尺寸，sig_diff 为公式6 中的高斯函数所需要的方差
sig_diff = sqrtf( std::max(sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA, 0.01f) );
//利用上面提到的公式a对图像进行高斯平滑处理
GaussianBlur(gray_fpt, gray_fpt, Size(), sig_diff, sig_diff);
//输出图像矩阵
return gray_fpt;
}
}2、构建高斯金字塔
//输入一张初始图像
涉及到公式b：
void SIFT::buildGaussianPyramid( const Mat& base, vector<Mat>& pyr, int nOctaves ) const
{
//向量数组sig 表示每组中计算各层图像所需的方差，nOctaveLayers + 3 即为S = s+3;
vector<double> sig(nOctaveLayers + 3);
//定义高斯金字塔的总层数，nOctaves*(nOctaveLayers + 3)即组数×层数
pyr.resize(nOctaves*(nOctaveLayers + 3));
//提前计算好各层图像所需的方差:
// \sigma_{total}^2 = \sigma_{i}^2 + \sigma_{i-1}^2
//第一层图像的尺度为基准层尺度σ0
sig[0] = sigma;
//每一层的尺度倍数k
double k = pow( 2., 1. / nOctaveLayers );
//从第1层开始计算所有的尺度
for( int i = 1; i < nOctaveLayers + 3; i++ )
{
//由上面的公式计算前一层的尺度
double sig_prev = pow(k, (double)(i-1))*sigma;
//计算当前尺度
double sig_total = sig_prev*k;
//计算公式a 中高斯函数所需的方差，并存入sig 数组内
sig[i] = std::sqrt(sig_total*sig_total - sig_prev*sig_prev);
}
//遍历高斯金字塔的所有层，构建高斯金字塔
for( int o = 0; o < nOctaves; o++ )
{
for( int i = 0; i < nOctaveLayers + 3; i++ )
{
//dst 为当前层图像矩阵
Mat& dst = pyr[o*(nOctaveLayers + 3) + i];
//如果当前层为高斯金字塔的第0 组第0 层，则直接赋值
if( o == 0 && i == 0 )
//把由createInitialImage 函数得到的基层图像矩阵赋予该层
dst = base;
// 新组的第0层需要从之前组中进行降采样
//如果当前层是除了第0 组以外的其他组中的第0 层，则要进行降采样处理
else if( i == 0 )
{
//提取出当前层所在组的前一组中的倒数第3 层图像
const Mat& src = pyr[(o-1)*(nOctaveLayers + 3) + nOctaveLayers];
//隔点降采样处理
resize(src, dst, Size(src.cols/2, src.rows/2),
0, 0, INTER_NEAREST);
}
//除了以上两种情况以外的其他情况的处理，即正常的组内生成不同层
else
{
//提取出当前层的前一层图像
const Mat& src = pyr[o*(nOctaveLayers + 3) + i-1];
//根据公式3，由前一层尺度图像得到当前层的尺度图像
GaussianBlur(src, dst, Size(), sig[i], sig[i]);
}
}
}
3、构建DoG 金字塔：
根据第2步得到的高斯金字塔，这里得到DoG金字塔
void SIFT::buildDoGPyramid( const vector<Mat>& gpyr, vector<Mat>& dogpyr ) const
{
//计算金字塔的组的数量
int nOctaves = (int)gpyr.size()/(nOctaveLayers + 3);
//定义DoG 金字塔的总层数，DoG 金字塔比高斯金字塔每组少一层
dogpyr.resize( nOctaves*(nOctaveLayers + 2) );
//遍历DoG 的所有层，构建DoG 金字塔
for( int o = 0; o < nOctaves; o++ )
{
for( int i = 0; i < nOctaveLayers + 2; i++ )
{
//提取出高斯金字塔的当前层图像
const Mat& src1 = gpyr[o*(nOctaveLayers + 3) + i];
//提取出高斯金字塔的上层图像
const Mat& src2 = gpyr[o*(nOctaveLayers + 3) + i + 1];
//提取出DoG 金字塔的当前层图像
Mat& dst = dogpyr[o*(nOctaveLayers + 2) + i];
//DoG 金字塔的当前层图像等于高斯金字塔的当前层图像减去高斯金字塔的上层图像
subtract(src2, src1, dst, noArray(), DataType<sift_wt>::type);
}
}
4、在DoG 尺度空间内找到极值点
在DoG尺度空间中找极值点，并按照主曲率的对比度删除不好的极值点
void SIFT::findScaleSpaceExtrema( const vector<Mat>& gauss_pyr, const vector<Mat>& dog_pyr,
vector<KeyPoint>& keypoints ) const
{
//得到金字塔的组数
int nOctaves = (int)gauss_pyr.size()/(nOctaveLayers + 3);
//SIFT_FIXPT_SCALE = 1，设定一个阈值用于判断在DoG 尺度图像中像素的大小
int threshold = cvFloor(0.5 * contrastThreshold / nOctaveLayers * 255 * SIFT_FIXPT_SCALE);
//SIFT_ORI_HIST_BINS = 36，定义梯度方向直方图的柱的数量
const int n = SIFT_ORI_HIST_BINS;
//定义梯度方向直方图变量
float hist[n];
KeyPoint kpt;
keypoints.clear();//清空
for( int o = 0; o < nOctaves; o++ )
for( int i = 1; i <= nOctaveLayers; i++ )
{
//DoG 金字塔的当前层索引
int idx = o*(nOctaveLayers+2)+i;
//DoG 金字塔当前层的尺度图像
const Mat& img = dog_pyr[idx];
//DoG 金字塔下层的尺度图像
const Mat& prev = dog_pyr[idx-1];
//DoG 金字塔上层的尺度图像
const Mat& next = dog_pyr[idx+1];
int step = (int)img.step1();
//图像的长和宽
int rows = img.rows, cols = img.cols;
//SIFT_IMG_BORDER = 5，该变量的作用是保留一部分图像的四周边界
for( int r = SIFT_IMG_BORDER; r < rows-SIFT_IMG_BORDER; r++)
{
//DoG 金字塔当前层图像的当前行指针
const sift_wt* currptr = img.ptr<sift_wt>(r);
//DoG 金字塔下层图像的当前行指针
const sift_wt* prevptr = prev.ptr<sift_wt>(r);
//DoG 金字塔上层图像的当前行指针
const sift_wt* nextptr = next.ptr<sift_wt>(r);
for( int c = SIFT_IMG_BORDER; c < cols-SIFT_IMG_BORDER; c++)
{
//DoG 金字塔当前层尺度图像的像素值
sift_wt val = currptr[c];
// find local extrema with pixel accuracy
//精确定位局部极值点
//如果满足if 条件，则找到了极值点，即候选特征点
if( std::abs(val) > threshold &&
//像素值要大于一定的阈值才稳定，即要具有较强的对比度
//下面的逻辑判断被“与”分为两个部分，前一个部分要满足像素值大于0，在3×3×3 的立方体内与周围26 个邻近像素比较找极大值，
//后一个部分要满足像素值小于0，找极小值
((val > 0 && val >= currptr[c-1] && val >= currptr[c+1] &&
val >= currptr[c-step-1] && val >= currptr[c-step] && val >= currptr[c-step+1] &&
val >= currptr[c+step-1] && val >= currptr[c+step] && val >= currptr[c+step+1] &&
val >= nextptr[c] && val >= nextptr[c-1] && val >= nextptr[c+1] &&
val >= nextptr[c-step-1] && val >= nextptr[c-step] && val >= nextptr[c-step+1] &&
val >= nextptr[c+step-1] && val >= nextptr[c+step] && val >= nextptr[c+step+1] &&
val >= prevptr[c] && val >= prevptr[c-1] && val >= prevptr[c+1] &&
val >= prevptr[c-step-1] && val >= prevptr[c-step] && val >= prevptr[c-step+1] &&
val >= prevptr[c+step-1] && val >= prevptr[c+step] && val >= prevptr[c+step+1]) ||
(val < 0 && val <= currptr[c-1] && val <= currptr[c+1] &&
val <= currptr[c-step-1] && val <= currptr[c-step] && val <= currptr[c-step+1] &&
val <= currptr[c+step-1] && val <= currptr[c+step] && val <= currptr[c+step+1] &&
val <= nextptr[c] && val <= nextptr[c-1] && val <= nextptr[c+1] &&
val <= nextptr[c-step-1] && val <= nextptr[c-step] && val <= nextptr[c-step+1] &&
val <= nextptr[c+step-1] && val <= nextptr[c+step] && val <= nextptr[c+step+1] &&
val <= prevptr[c] && val <= prevptr[c-1] && val <= prevptr[c+1] &&
val <= prevptr[c-step-1] && val <= prevptr[c-step] && val <= prevptr[c-step+1] &&
val <= prevptr[c+step-1] && val <= prevptr[c+step] && val <= prevptr[c+step+1])))
{
//三维坐标，长、宽、层（层与尺度相对应）
int r1 = r, c1 = c, layer = i;
// adjustLocalExtrema 函数的作用是调整局部极值的位置，即找到亚像素级精度的特征点，该函数在后面给出详细的分析
//如果满足if 条件，说明该极值点不是特征点，继续上面的for循环
if( !adjustLocalExtrema(dog_pyr, kpt, o, layer, r1, c1,
nOctaveLayers, (float)contrastThreshold,
(float)edgeThreshold, (float)sigma) )
continue;
//计算特征点相对于它所在组的基准层的尺度，即公式5 所得尺度
float scl_octv = kpt.size*0.5f/(1 << o);
// calcOrientationHist 函数为计算特征点的方向角度，后面给出该
函数的详细分析
//SIFT_ORI_SIG_FCTR = 1.5f ， SIFT_ORI_RADIUS = 3 *
SIFT_ORI_SIG_FCTR
float omax = calcOrientationHist(gauss_pyr[o*(nOctaveLayers+3) + layer],
Point(c1, r1),
cvRound(SIFT_ORI_RADIUS * scl_octv),
SIFT_ORI_SIG_FCTR * scl_octv,
hist, n);
//SIFT_ORI_PEAK_RATIO = 0.8f
//计算直方图辅方向的阈值，即主方向的80%
float mag_thr = (float)(omax * SIFT_ORI_PEAK_RATIO);
//计算特征点的方向
for( int j = 
最低0.47元/天 解锁文章
扫一扫
14万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
