一：构造函数ORBextractor()
说明:
特征点提取器的构造函数
该函数位于Tracking.cc的构造函数中Tracking::Tracking(...)，该构造函数主要从配置文件中加载相机参数,并计算ORB特征点有关的参数,并新建特征点提取器
mpORBextractorLeft = new ORBextractor(nFeatures,fScaleFactor,nLevels,fIniThFAST,fMinThFAST);
参数说明
int nfeatures | 整个图像金字塔中，要提取的特征点数目 |
|---|---|
double scaleFactor | 图像金字塔层与层之间的缩放因子 |
int nlevels | 图像金字塔的层数 |
int iniThFAST | 初始的FAST响应值阈值 |
int minThFAST | 最小的FAST响应值阈值 |
vector<int> mnFeaturesPerLevel | 分配到每层图像中，要提取的特征点数目 |
|---|---|
vector<int> umax | 1/4圆形的图像区域每行u轴的边界 |
vector<float> mvScaleFactor | 每层图像的缩放因子 |
vector<float> mvInvScaleFactor | 每层缩放因子的倒数- |
vector<float> mvLevelSigma2 | 每层图像的缩放因子的平方 |
vector<float> mvInvLevelSigma2 | 每层图像的缩放因子的平方的倒数 |
函数流程
构造函数列表初始化
初始化nfeatures，scaleFactor，nlevels，iniThFAST，minThFAST这几个参数
ORBextractor::ORBextractor(int _nfeatures,float _scaleFactor,int _nlevels,int _iniThFAST, int _minThFAST):
nfeatures(_nfeatures), scaleFactor(_scaleFactor), nlevels(_nlevels),iniThFAST(_iniThFAST), minThFAST(_minThFAST)
计算每层金字塔的缩放系数和它的几个好哥俩
mvScaleFactor[0] | 1.0 |
|---|---|
mvScaleFactor[1] | 1.2 |
mvScaleFactor[2] | 1.2 * 1.2 |
mvScaleFactor[3 ] | 1.2 * 1.2 * 1.2 |
| … | … |
for(int i=1; i<nlevels; i++)
{
mvScaleFactor[i]=mvScaleFactor[i-1]*scaleFactor;
mvLevelSigma2[i]=mvScaleFactor[i]*mvScaleFactor[i];
}
for(int i=0; i<nlevels; i++)
{
mvInvScaleFactor[i]=1.0f/mvScaleFactor[i];
mvInvLevelSigma2[i]=1.0f/mvLevelSigma2[i];
}
计算每层金字塔需要提取的特征点的数目
金字塔层数越高，图像面积越小，所分配的特征点数目就越少。
每层图像分配的特征点数目是成比例的，比例系数为scaleFactor(1.2)
设N为特征点总数,n为层数，则第0层分配的特征点数目为
N
0
=
N
(
1
−
s
)
(
1
−
s
n
)
N0 = \frac{N(1-s)}{(1-s^{n})}
N0=(1−sn)N(1−s)
N0也就是代码中的nDesiredFeaturesPerScale，然后依次按比例计算1-6层特征点数目，第7层数目是总数减去前6层数目。
mnFeaturesPerLevel.resize(nlevels);
float factor = 1.0f / scaleFactor;
//第0层图像应该分配的特征点数量
float nDesiredFeaturesPerScale = nfeatures*(1 - factor)/(1 - (float)pow((double)factor, (double)nlevels));
int sumFeatures = 0;
//先不计算顶层图像
for( int level = 0; level < nlevels-1; level++ )
{
mnFeaturesPerLevel[level] = cvRound(nDesiredFeaturesPerScale);
sumFeatures += mnFeaturesPerLevel[level];
nDesiredFeaturesPerScale *= factor;
}
//由于前面的特征点个数取整操作，可能会导致剩余一些特征点个数没有被分配，所以这里就将这个余出来的特征点分配到最高的图层中
mnFeaturesPerLevel[nlevels-1] = std::max(nfeatures - sumFeatures, 0);
计算1/4圆u坐标最大值
这是为了使用灰度质心法实现BRIEF描述子的旋转不变性的。
何为灰度质心？灰度质心就是一个图像区域内计算的，以图像快灰度值作为权重的中心（就像一块质量分布不均匀的平板一样计算它的重心一样）。
暂且先不关心灰度质心法和旋转不变性的含义，先来看上面所说的图像区域，实际上是一个圆，称之为灰度质心圆。为了方便后续灰度质心的计算，首先要求出umax
参数说明
const int HALF_PATCH_SIZE = 15 | 使用灰度质心法计算特征点的方向信息时，图像块的半径 |
|---|---|
vector<int> umax | 1/4灰度质心圆的每一行的u坐标边界 |
umax.resize(HALF_PATCH_SIZE + 1);
int v,v0,vmax = cvFloor(HALF_PATCH_SIZE * sqrt(2.f) / 2 + 1);
//这里的二分之根2就是对应那个45°圆心角
int vmin = cvCeil(HALF_PATCH_SIZE * sqrt(2.f) / 2);
const double hp2 = HALF_PATCH_SIZE*HALF_PATCH_SIZE;//半径的平方
//利用圆的方程计算每行像素的u坐标边界（max）
for (v = 0; v <= vmax; ++v)
umax[v] = cvRound(sqrt(hp2 - v * v));
//使用了对称的方式计算上八分之一的圆周上的umax
for (v = HALF_PATCH_SIZE, v0 = 0; v >= vmin; --v)
{
while (umax[v0] == umax[v0 + 1])
++v0;
umax[v] = v0;
++v0;
}
然而图像的像素坐标都是整数值，求解完的图像区域只能是一个近似的圆。所以为了尽量提高对称性，使选取的点的区域既关于两个坐标轴对称，也关于u=v对称。
第一个for循环计算时由于四舍五入有连续相同的u坐标的点，比如v从0增加到3，umax[0]=umax[3]=15保持不变；第二个for循环就是先判断u坐标连续相同的次数，保持v=15使u从0增到3，那么umax[15]=3，依次类推使围成的点域关于u=v对称
计算之后的效果大概是这样的
二：提取特征点函数
说明 :
提取图像的ORB特征点，提取的关键点存放在mvKeys，描述子存放在mDescriptors,对应于仿函数的_keypoints和_descriptors
void Frame::ExtractORB(int flag, const cv::Mat &im)
{
if(flag==0)
// 左图的话就套使用左图指定的特征点提取器，并将提取结果保存到对应的变量中
// 仿函数，重载括号运算符
(*mpORBextractorLeft)(im,cv::Mat(),mvKeys, mDescriptors);
else
// 右图的话就需要使用右图指定的特征点提取器，并将提取结果保存到对应的变量中
(*mpORBextractorRight)(im,cv::Mat(),mvKeysRight,mDescriptorsRight);
}
使用仿函数完成提取特征点的工作
输出
vector<KeyPoint>& _keypoints各个图像金字塔中特征点坐标都转化成第0层的wholeSize尺寸下特征点坐标，
_descriptors存储整个金字塔的特征点，一列（32维）对应一个特征点。
函数流程
void ORBextractor::operator()( InputArray _image, InputArray _mask, vector<KeyPoint>& _keypoints,OutputArray _descriptors){
/**
* @param[in] _image 输入原始图的图像
* @param[in] _mask 掩膜mask
* @param[in & out] _keypoints 输出：存储特征点关键点的向量
* @param[in & out] _descriptors 输出：存储特征点描述子的矩阵
*/
if(_image.empty())
return;
//获取图像的大小
Mat image = _image.getMat();
//判断图像的格式是否正确，要求是单通道灰度值
assert(image.type() == CV_8UC1 );
}
1.构建图像金字塔
概念：图像金字塔是图像多尺度表达的一种，是一种以多分辨率来解释图像的有效但概念简单的结构。一幅图像的图像金字塔是一系列以金字塔形状（自下而上）逐步降低，且来源于同一张原始图的图像分辨率集合。其通过梯次向下采样获得，直到达到某个终止条件才停止采样。我们将一层一层的图像比喻成金字塔，层级越高，则图像越小，分辨率越低。来自于百度百科
在相机移动过程中，相机距离观测的目标有时近有时远，人们直观的理解物体是远小近大的，但是相机的分辨率是固定的，导致远处的图像所提取的特征点在近处观测可能不再是特征点了（远看一个点，近看是个圆），但是明明观测的是同一个物体，却因为尺度的变化使两张图像无法匹配上，干着急。
为了解决这个问题，构建图像金字塔。其思路就是对原始图像依次按某个比例系数缩放比例进行降采样得到共计8张图片，分别对每张图片进行特征提取，并记录特征所在金字塔的第几层。
图像金字塔示例
ComputePyramid(image);
2.计算图像的特征点，并且将特征点进行均匀化。均匀的特征点可以提高位姿计算精度
// 存储所有的特征点，第一维存储金字塔层数，第二维存储该层金字塔里提取的所有特征点
vector < vector<KeyPoint> > allKeypoints;
//使用四叉树的方式计算每层图像的特征点并进行分配
ComputeKeyPointsOctTree(allKeypoints);
3.拷贝图像描述子到新的矩阵descriptors
4.执行下面的for循环，遍历每一层图像
//设置了Offset变量来保存“寻址”时的行偏移量，辅助进行在总描述子mat中的定位
int offset = 0;
for (int level = 0; level < nlevels; ++level)
{
//获取在allKeypoints中当前层特征点容器的句柄
vector<KeyPoint>& keypoints = allKeypoints[level];
//本层的特征点数
int nkeypointsLevel = (int)keypoints.size();
//如果特征点数目为0，跳出本次循环，继续下一层金字塔
if(nkeypointsLevel==0)
continue;
... ...对图像进行高斯模糊... ...
... ...计算高斯模糊后图像的描述子... ...
... ...对非第0层图像中的特征点的坐标恢复到第0层图像（原图像）的坐标系下... ...
}
5 对图像进行高斯模糊
为减少噪声干扰，先对图像进行高斯滤波，这样后续计算的描述子才更加准确。
注意：提取特征点的时候，使用的是清晰的原图像；这里计算描述子的时候，为了避免图像噪声的影响，使用了高斯模糊
GaussianBlur(workingMat, //源图像
workingMat, //输出图像
Size(7, 7), //高斯滤波器kernel大小，必须为正的奇数
2, //高斯滤波在x方向的标准差
2, //高斯滤波在y方向的标准差
BORDER_REFLECT_101);//边缘拓展点插值类型
6.计算高斯模糊后图像的描述子
// desc存储当前图层的描述子
Mat desc = descriptors.rowRange(offset, offset + nkeypointsLevel);
// Step
computeDescriptors(workingMat, //高斯模糊之后的图层图像
keypoints, //当前图层中的特征点集合
desc, //存储计算之后的描述子
pattern); //随机采样模板
// 更新偏移量的值
offset += nkeypointsLevel;
7.对非第0层图像中的特征点的坐标恢复到第0层图像坐标系下
// Step 6 对非第0层图像中的特征点的坐标恢复到第0层图像（原图像）的坐标系下
// 对于第0层的图像特征点，他们的坐标就不需要再进行恢复了
if (level != 0)
{
// 获取当前图层上的缩放系数
float scale = mvScaleFactor[level];
// 遍历本层所有的特征点
for (vector<KeyPoint>::iterator keypoint = keypoints.begin(),
keypointEnd = keypoints.end(); keypoint != keypointEnd; ++keypoint)
// 特征点本身直接乘缩放倍数就可以了
keypoint->pt *= scale;
}
// 将keypoints中内容插入到_keypoints 的末尾
// keypoint其实是对allkeypoints中每层图像中特征点的引用，这样allkeypoints中的所有特征点在这里被转存到输出的_keypoints
_keypoints.insert(_keypoints.end(), keypoints.begin(), keypoints.end());
相关函数
1 构建金字塔
对于第0层图像，直接将原图image拷贝到temp的中间，并对其四周进行对称式的边界扩展；对于其他层的图像，先将原图按照当前层的缩放系数sz缩小，再将缩小后的图像拷贝到当前层的扩展边界图像temp中，并对其四周进行对称式的边界扩展。
其中，定义的Size wholeSize是图像进行扩充之后的尺寸
在原函数基础上作一定修改
void ORBextractor::ComputePyramid(cv::Mat image)
{
//开始遍历所有的图层
for (int level = 0; level < nlevels; ++level){
//获取本层图像的缩放系数
float scale = mvInvScaleFactor[level];
//计算本层图像的像素尺寸大小
Size sz(cvRound((float)image.cols*scale), cvRound((float)image.rows*scale));
//全尺寸图像。包括无效图像区域的大小。将图像进行“补边”，EDGE_THRESHOLD区域外的图像不进行FAST角点检测
Size wholeSize(sz.width + EDGE_THRESHOLD*2, sz.height + EDGE_THRESHOLD*2);
// 定义了两个变量：temp是扩展了边界的图像，masktemp并未使用
Mat temp(wholeSize, image.type()), masktemp;
// mvImagePyramid 刚开始时是个空的vector<Mat>
// 我把下面这句代码注释掉了，因为写不写不影响最后的金字塔结果
//mvImagePyramid[level] = temp(Rect(EDGE_THRESHOLD, EDGE_THRESHOLD, sz.width, sz.height));
//计算第0层以上resize后的图像
if( level != 0 ){
resize(image, //输入图像，源代码这个参数是mvImagePyramid[level-1]，显然是不正确的
mvImagePyramid[level], sz, 0,0,cv::INTER_LINEAR);
//把源图像拷贝到目的图像的中央，四面填充指定的像素。图片如果已经拷贝到中间，只填充边界
copyMakeBorder(mvImagePyramid[level], //源图像
temp, //目标图像（此时其实就已经有大了一圈的尺寸了）
EDGE_THRESHOLD, EDGE_THRESHOLD, //top & bottom 需要扩展的border大小
EDGE_THRESHOLD, EDGE_THRESHOLD, //left & right 需要扩展的border大小
BORDER_REFLECT_101+BORDER_ISOLATED); //扩充方式
}
else{
copyMakeBorder(image,temp,EDGE_THRESHOLD,EDGE_THRESHOLD,
EDGE_THRESHOLD,EDGE_THRESHOLD,BORDER_REFLECT_101);
}
//此处加上这句，将扩展边后的temp赋给金字塔
mvImagePyramid[level] = temp;
}
}
//几种 copyMakeBorder函数图像扩充的方式
/*Various border types, image boundaries are denoted with '|'
* BORDER_REPLICATE: aaaaaa|abcdefgh|hhhhhhh
* BORDER_REFLECT: fedcba|abcdefgh|hgfedcb
* BORDER_REFLECT_101: gfedcb|abcdefgh|gfedcba
* BORDER_WRAP: cdefgh|abcdefgh|abcdefg
* BORDER_CONSTANT: iiiiii|abcdefgh|iiiiiii with some specified 'i'
*/
最后结果示意图，每层图像在原图基础上四边都扩充了EDGE_THRESHOLD
构
2 提取特征点
函数大概流程分析
将每层金字塔分成若干个cell，分别在每个cell中提取特征点，然后采用四叉树对特征点均匀化处理。
对于每一层图像，计算这层图像的有效边界，请看图
每层金字塔的边界是上图的最外围边框，计算金字塔函数中称之为wholeSize，生成这个外框的目的是进行图像金子塔的生成时，需要对图像进行高斯滤波处理，最里面区域是原图像，而中间的绿色边框则是在原图基础上四边各扩充了3得到的。给它起个名字叫“有效图像边界”，扩充3的目的，是因为在计算FAST特征点的时候，需要建立一个半径为3的圆，如图所示
首先计算有效图像边界的横纵坐标最小值以及最大值，分别是minBorderX、minBorderY、maxBorderX、maxBorderY，计算有效图像边界的宽和高，分别为width、height，进而计算实际每个cell所占像素的行数和列数nCols、nRows，注意行数和列数是大于等于30的，一般会是30（概率很小）或31或32。
下面开始对每个cell进行遍历，首先计算了iniY和maxY，以及iniX和maxX，它们是每个网格实际的边界，在width * height的基础上分别扩充3，以计算FAST特征点，而width * height是可以提取特征点的区域。此外，下一个网格的iniY比当前的maxY小6，也就是有一部分重叠，这样才可以使提取FAST特征点的区域是连续的。
函数的一个小bug：判断行初始坐标是否有效时时if(iniY>=maxBorderY-3)，判断
列初始坐标是否有效时 if(iniX>=maxBorderX-6)，为什么减的不一样呢？答案应该是相同的。那么应该是-3还是-6呢，我的理解是-6，原因如下：
第i+1个格子的iniY比第i个格子的maxY小6，那么如果第n个格子的iniY大于maxBorderY-6的话，那说明第n-1个格子的maxY已经大于maxBorderY了，已经到了有效图像边界的最右边了，“已经结束啦”。-6包括-3的情况，而且是合理的判断条件，因此认为iniY的判断条件是if(iniY>=maxBorderY-6)
示意图：
在网格中提取的特征点存放在 vToDistributeKeys中，其坐标是基于图像坐标边界的，调用DistributeOctTree将筛选后的特征点保存在keypoints中，vector<KeyPoint> & keypoints = allKeypoints[level]。最后 allKeypoints的坐标是基于所在层的wholeSize而言的。
//计算四叉树的特征点
void ORBextractor::ComputeKeyPointsOctTree(vector<vector<KeyPoint> >& allKeypoints)
{
allKeypoints.resize(nlevels);
//图像cell的尺寸，是个正方形
const float W = 30;
// 对每一层图像做处理
for (int level = 0; level < nlevels; ++level)
{
const int minBorderX = EDGE_THRESHOLD-3;
const int minBorderY = minBorderX;
const int maxBorderX = mvImagePyramid[level].cols-EDGE_THRESHOLD+3;
const int maxBorderY = mvImagePyramid[level].rows-EDGE_THRESHOLD+3;
//存储需要进行平均分配的特征点
vector<cv::KeyPoint> vToDistributeKeys;
vToDistributeKeys.reserve(nfeatures*10);
//计算进行特征点提取的图像区域尺寸
const float width = (maxBorderX-minBorderX);
const float height = (maxBorderY-minBorderY);
//计算网格在当前层的图像有的行数和列数
const int nCols = width/W;
const int nRows = height/W;
//计算每个图像网格所占的像素行数和列数
const int wCell = ceil(width/nCols);
const int hCell = ceil(height/nRows);
//开始遍历图像网格，还是以行开始遍历的
for(int i=0; i<nRows; i++){
const float iniY =minBorderY+i*hCell;//计算当前网格初始行坐标
float maxY = iniY+hCell+6; //最大行坐标
//如果初始的行坐标就已经超过了maxBorderY-3？
if(iniY>=maxBorderY-3)
continue;
//如果图像的大小导致不能够正好划分出来整齐的图像网格，那么就要委屈最后一行了
if(maxY>maxBorderY)
maxY = maxBorderY;
//开始列的遍历
for(int j=0; j<nCols; j++){
const float iniX =minBorderX+j*wCell;//计算初始的列坐标
float maxX = iniX+wCell+6;//这列网格的最大列坐标
//如果初始的列坐标就已经超过了maxBorderX-6？
if(iniX>=maxBorderX-6)
continue;
//如果最大坐标越界那么委屈一下
if(maxX>maxBorderX)
maxX = maxBorderX;
vector<cv::KeyPoint> vKeysCell;//这个向量存储这个cell中的特征点
FAST(mvImagePyramid[level].rowRange(iniY,maxY).colRange(iniX,maxX), //待检测的图像，这里就是当前遍历到的图像块
vKeysCell, //存储角点位置的容器
iniThFAST, //检测阈值
true); //使能非极大值抑制
//如果这个图像块中使用默认的FAST检测阈值没有能够检测到角点
if(vKeysCell.empty()){
//那么就使用更低的阈值来进行重新检测
FAST(mvImagePyramid[level].rowRange(iniY,maxY).colRange(iniX,maxX),
vKeysCell,minThFAST,true);
}
//当图像cell中检测到FAST角点的时候执行下面的语句
if(!vKeysCell.empty()){
//遍历其中的所有FAST角点
for(vector<cv::KeyPoint>::iterator vit=vKeysCell.begin(); vit!=vKeysCell.end();vit++){
//到目前为止，这些角点的坐标都是基于图像cell的，先将其恢复到图像坐标边界下的坐标
(*vit).pt.x+=j*wCell;
(*vit).pt.y+=i*hCell;
//然后将其加入到”等待被分配“的特征点容器中
vToDistributeKeys.push_back(*vit);
}
}
}
}
//声明一个对当前图层的特征点的容器的引用
vector<KeyPoint> & keypoints = allKeypoints[level];
keypoints.reserve(nfeatures);
//返回值是一个保存有特征点的vector容器，含有剔除后的保留下来的特征点
//得到的特征点的坐标，依旧是在当前图层下来讲的
keypoints = DistributeOctTree(vToDistributeKeys, minBorderX, maxBorderX，
minBorderY, maxBorderY,mnFeaturesPerLevel[level],level);
//PATCH_SIZE是对于底层的初始图像来说的，现在要根据当前图层的尺度缩放倍数进行缩放得到缩放后的PATCH大小 和特征点的方向计算有关
const int scaledPatchSize = PATCH_SIZE*mvScaleFactor[level];
//获取剔除过程后保留下来的特征点数目
const int nkps = keypoints.size();
//然后开始遍历这些特征点，恢复其在当前图层图像坐标系下的坐标
for(int i=0; i<nkps ; i++){
//对每一个保留下来的特征点，恢复到相对于当前图层“边缘扩充图像下”的坐标系的坐标
keypoints[i].pt.x+=minBorderX;
keypoints[i].pt.y+=minBorderY;
//记录特征点来源的图像金字塔图层
keypoints[i].octave=level;
//记录计算方向的patch，缩放后对应的大小， 又被称作为特征点半径
keypoints[i].size = scaledPatchSize;
}
}//对于每层图像执行上述操作
//然后计算这些特征点的方向信息，注意这里还是分层计算的
for (int level = 0; level < nlevels; ++level)
computeOrientation(mvImagePyramid[level], //对应的图层的图像
allKeypoints[level], //这个图层中提取并保留下来的特征点容器
umax); //以及PATCH的横坐标边界
}
3 函数中调用的DistributeOctTree
说明 使用四叉树法对一个图像金字塔图层中的特征点进行平均和分发
传入vToDistributeKeys，返回 keypoints（每层图像均匀化后的特征点，其坐标仍然是基于有效图像边界的）
原理：
- 如果图片的宽度比较宽，就先把分成左右w/h份。一般的640×480的图像开始的时候只有一个node。
- 如果node里面的点数>1，把每个node分成四个node，如果node里面的特征点为空，就不要了，删掉。
- 新分的node的点数>1，就再分裂成4个node。如此，一直分裂。
- 终止条件为：node的总数量> [公式] ，或者无法再进行分裂。
- 然后从每个node里面选择一个质量最好的FAST点。
参考知乎的一篇文章使用四叉树法对特征点进行平均和分发
4 函数调用的computeOrientation
说明用于计算特征点的方向
for (vector<KeyPoint>::iterator kpit = keypoints.begin();kpit !=keypoints.end(); ++keypoint)
{
// 调用IC_Angle 函数计算这个特征点的方向
keypoint->angle = IC_Angle(image, //特征点所在的图层的图像
keypoint->pt, //特征点在这张图像中的坐标
umax); //每个特征点所在图像区块的每行的边界 u_max 组成的vector
IC_Angle()函数
该函数计算了特征点的方向，返回是一个角度值。
计算特征点方向是为了使提取的特征点具有旋转不变性。
计算特征点方向的方式是灰度质心法（Intensity Centriod）
个人理解是为了使描述子具有旋转不变形，因为计算描述子的pattern是固定的，所有图像都用这个pattern,所有人都用。如果图像旋转了一个角度，pattern还在原地踏步，计算出来的描述子和旋转之前对应的特征点的描述子就很不一样了（你站着的时候我认识你，现在不能因为你趴地上了我就认不出来你了），应该说即使图像旋转了，对应的特征点的描述子也应该是相似的，于是通过计算关键点的角度，即使图像发生旋转，让pattern也跟着关键点一起旋转，最后计算的描述子的距离才比较接近，匹配更准。因此说，计算了关键点角度以此计算的描述子才叫做Steered BRIEF。
灰度质心法原理相信大家应该比较熟悉了，十四讲啥的上面也写的很清楚了，或者看这篇文章灰度质心法原理
这里之前计算的u_max就派上用场了，我们处理的这块图像区域近似是一个圆形，是一行一行计算的，首先计算v=0这行的图像矩，然后依次计算关于u=0对称的两行像素的图像矩，两行两行计算，最后累加即可。
因为计算过程是按行索引的，因此需要知道每行像素坐标范围，还好，之前算过了。就是[-umax[v],umax[v]]
这里的const uchar *center是一个指向图像某坐标点灰度值的指针，本代码中指向的点是特征点的坐标。center[u]取出了[pt.(x+u),pt.y]点处的灰度值（取的都是v=0这个坐标轴的灰度值）。
那么怎么取出[pt.(x+u),pt.(y+v)]的像素值呢？使用(int)image.step1()表示图像一行的字节总数，其像素值为center[u + v * step]（在center[u]移动了v行）。其他的问题就都很好理解了。
static float IC_Angle(const Mat &image, Point2f pt, const vector<int> &u_max)
{
//图像的矩，前者是按照图像块的y坐标加权，后者是按照图像块的x坐标加权
int m_01 = 0, m_10 = 0;
//获得这个特征点所在的图像块的中心点坐标灰度值的指针center
const uchar *center = &image.at<uchar>(cvRound(pt.y), cvRound(pt.x));
//这条v=0中心线的计算需要特殊对待
//后面是以中心行为对称轴，成对遍历行数，所以PATCH_SIZE必须是奇数
for (int u = -HALF_PATCH_SIZE; u <= HALF_PATCH_SIZE; ++u)
m_10 += u * center[u];
//这里的step1表示这个图像一行包含的字节总数。参考
int step = (int)image.step1();
//注意这里是以v=0中心线为对称轴，然后对称地每成对的两行之间进行遍历，这样处理加快了计算速度
for (int v = 1; v <= HALF_PATCH_SIZE; ++v)
{
int v_sum = 0;
int d = u_max[v];
// 假设每次处理的两个点坐标，中心线下方为(x,y),中心线上方为(x,-y)
// 对于某次待处理的两个点：m_10 = Σ x*I(x,y) = x*I(x,y) + x*I(x,-y) = x*(I(x,y) + I(x,-y))
// 对于某次待处理的两个点：m_01 = Σ y*I(x,y) = y*I(x,y) - y*I(x,-y) = y*(I(x,y) - I(x,-y))
for (int u = -d; u <= d; ++u)
{
int val_plus = center[u + v * step], val_minus = center[u - v * step];
//在v（y轴）上，2行所有像素灰度值之差
v_sum += (val_plus - val_minus);
//u轴（x轴）方向上用u坐标加权和，相当于同时计算两行
m_10 += u * (val_plus + val_minus);
}
//将这一行上的和按照y坐标加权
m_01 += v * v_sum;
}
//输出为[0,360)角度，精度为0.3°
return fastAtan2((float)m_01, (float)m_10);
}
最终的返回值（角度）赋给了 keypoint->angle
5 函数中的computeDescriptors
说明 全局静态函数，static修饰符限定其只能够被本文件中的函数调用,计算某层金字塔图像上特征点的描述子
BRIEF描述子定义_descriptors.create(nkeypoints, 32,CV_8U)，可见描述子一共有nkeypoints行，每行（也就是每个特征点）对应一个描述子，描述子是32维的向量，每个元素存储了8bit的信息(CV_8U)
计算描述子需要使用随机采样点集合pattern，它定义在ORBextractor中。将成员变量static int bit_pattern_31_[256 * 4]内容复制给vector<Point> pattern。
256*4表示的是每4个值组成一对点进行像素值的比较，得到1bit的数据（比出个0或者1），一共比出了256个bit，其中每8bit的信息就是BRIEF描述子的一个元素（前面说描述子元素类型是CV_8U），所以得到了32个BRIEF描述子的值，组成了_descriptors的一列，也就是一个特征点所对应的描述子。
而bit_pattern_31_这个数组都是int类型的，需要将其强制类型转换成cv::point类型，用来计算描述子。这就是下面这段代码所做的。
//成员变量pattern的长度，也就是点的个数，这里的512表示512个点（
const int npoints = 512;
//获取用于计算BRIEF描述子的随机采样点点集头指针
//注意到pattern0数据类型为Points*,bit_pattern_31_是int[]型，所以这里需要进行强制类型转换
const Point *pattern0 = (const Point *)bit_pattern_31_;
//将在全局变量区域的、int格式的随机采样点以cv::point格式复制到当前类对象中的成员变量中
std::copy(pattern0, pattern0 + npoints, std::back_inserter(pattern));
下面是计算描述子的函数，里面又调用了computeOrbDescriptor函数
static void computeDescriptors(const Mat &image, vector<KeyPoint> &keypoints,
Mat &descriptors,const vector<Point> &pattern){
for (size_t i = 0; i < keypoints.size(); i++)
computeOrbDescriptor(keypoints[i],image, &pattern[0],descriptors.ptr((int)i));
}
这里计算了描述子，先将pattern中的点对按照计算的角度进行旋转 x'= xcos(θ) - ysin(θ), y'= xsin(θ) + ycos(θ)，获取其像素值center[ y'* step + x']，进行两两比较，比较的结果。desc[i] = (uchar)val是一种8-bit无符号整形数据，范围为[0, 255].
结果类似于这样
const float factorPI = (float)(CV_PI / 180.f); //乘数因子，一度对应着多少弧度
static void computeOrbDescriptor(const KeyPoint &kpt, const Mat &img, const Point *pattern, uchar *desc)
{
float angle = (float)kpt.angle * factorPI; //得到特征点的角度，用弧度制表示。
float a = (float)cos(angle), b = (float)sin(angle);//计算这个角度的余弦值和正弦值
const uchar *center = &img.at<uchar>(cvRound(kpt.pt.y), cvRound(kpt.pt.x)); //获得图像中心指针
const int step = (int)img.step;//获得图像的每行的字节数
//获得采样点中某个idx所对应的点的灰度值,这里旋转前坐标为(x,y), 旋转后坐标(x',y')
// 下面表示 y'* step + x'
#define GET_VALUE(idx) center[cvRound(pattern[idx].x * b + pattern[idx].y * a) * step + cvRound(pattern[idx].x * a - pattern[idx].y * b)]
//其中每一位是来自于两个像素点灰度的直接比较，所以每比较出8bit结果，需要16个随机点，这也就是pattern需要+=16的原因
for (int i = 0; i < 32; ++i, pattern += 16)
{
int t0, //参与比较的第1个特征点的灰度值
t1, //参与比较的第2个特征点的灰度值
val; //描述子这个字节的比较结果，0或1
t0 = GET_VALUE(0);
t1 = GET_VALUE(1);
val = t0 < t1; //描述子本字节的bit0
t0 = GET_VALUE(2);
t1 = GET_VALUE(3);
val |= (t0 < t1) << 1; //描述子本字节的bit1
t0 = GET_VALUE(4);
t1 = GET_VALUE(5);
val |= (t0 < t1) << 2; //描述子本字节的bit2
t0 = GET_VALUE(6);
t1 = GET_VALUE(7);
val |= (t0 < t1) << 3; //描述子本字节的bit3
t0 = GET_VALUE(8);
t1 = GET_VALUE(9);
val |= (t0 < t1) << 4; //描述子本字节的bit4
t0 = GET_VALUE(10);
t1 = GET_VALUE(11);
val |= (t0 < t1) << 5; //描述子本字节的bit5
t0 = GET_VALUE(12);
t1 = GET_VALUE(13);
val |= (t0 < t1) << 6; //描述子本字节的bit6
t0 = GET_VALUE(14);
t1 = GET_VALUE(15);
val |= (t0 < t1) << 7; //描述子本字节的bit7
desc[i] = (uchar)val;
}
#undef GET_VALUE
}
432
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
