本次主要讲解的是ORB特征点的具体提取函数void ORBextractor::operator()
,这个函数完成了Fast角点提取,特征点均匀化,描述子计算等关键功能。这里的operator()
重载了()
运算符,使得这个类可以使用调用函数的方式进行调用,也称仿函数
,关于仿函数
的具体细节有兴趣的可以百度。
首先是函数的参数:
void ORBextractor::operator()( InputArray _image, // 输入图像
InputArray _mask, // 掩膜,这里没有用到
vector<KeyPoint>& _keypoints, // 特征点Vector
OutputArray _descriptors) // 保存描述子的矩阵
在这个仿函数中主要完成了以下功能:
-
检查输入图像的有效性
// Step 1 检查图像有效性。如果图像为空,那么就直接返回 if(_image.empty()) return; //获取图像 Mat image = _image.getMat(); //判断图像的格式是否正确,要求是单通道灰度值 assert(image.type() == CV_8UC1 );
-
构建图像金字塔
// Step 2 构建图像金字塔 ComputePyramid(image);
这个函数我们详细展开:
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指向temp的中间部分(这里为浅拷贝,内存相同) mvImagePyramid[level] = temp(Rect(EDGE_THRESHOLD, EDGE_THRESHOLD, sz.width, sz.height)); //计算第0层以上resize后的图像 if( level != 0 ) { //将上一层金字塔图像根据设定sz缩放到当前层级 resize(mvImagePyramid[level-1], //输入图像 mvImagePyramid[level], //输出图像 sz, //输出图像的尺寸 0, //水平方向上的缩放系数,留0表示自动计算 0, //垂直方向上的缩放系数,留0表示自动计算 cv::INTER_LINEAR); //图像缩放的差值算法类型,这里的是线性插值算法 // //! 原代码mvImagePyramid 并未扩充,上面resize应该改为如下 // resize(image, //输入图像 // mvImagePyramid[level], //输出图像 // sz, //输出图像的尺寸 // 0, //水平方向上的缩放系数,留0表示自动计算 // 0, //垂直方向上的缩放系数,留0表示自动计算 // cv::INTER_LINEAR); //图像缩放的差值算法类型,这里的是线性插值算法 //把源图像拷贝到目的图像的中央,四面填充指定的像素。图片如果已经拷贝到中间,只填充边界 //这样做是为了能够正确提取边界的FAST角点 //EDGE_THRESHOLD指的这个边界的宽度,由于这个边界之外的像素不是原图像素而是算法生成出来的,所以不能够在EDGE_THRESHOLD之外提取特征点 copyMakeBorder(mvImagePyramid[level], //源图像 temp, //目标图像(此时其实就已经有大了一圈的尺寸了) EDGE_THRESHOLD, EDGE_THRESHOLD, //top & bottom 需要扩展的border大小 EDGE_THRESHOLD, EDGE_THRESHOLD, //left & right 需要扩展的border大小 BORDER_REFLECT_101+BORDER_ISOLATED); //扩充方式 } else { //对于第0层未缩放图像,直接将图像深拷贝到temp的中间,并且对其周围进行边界扩展。此时temp就是对原图扩展后的图像 copyMakeBorder(image, //这里是原图像 temp, EDGE_THRESHOLD, EDGE_THRESHOLD, EDGE_THRESHOLD, EDGE_THRESHOLD, BORDER_REFLECT_101); } // //! 原代码mvImagePyramid 并未扩充,应该添加下面一行代码 // mvImagePyramid[level] = temp; } }
这里有几个知识点需要注意:
- 我们需要明白
EDGE_THRESHOLD
具体做了什么,如下图所示
其中深灰色的为原始图像,而浅灰色的为四个方向各扩充了EDGE_THRESHOLD
个像素的图像,绿色的为四个方向各扩充了3个像素的图像,这么做的主要原因是为了处理在图像边缘附近的特征点提取和描述子计算问题。 - 关于图像拓充方式,
opencv
给出了以下几种选项及对应的含义
- 这里作者写的有问题,如代码中注释,作者并没有将每层金字塔指向拓充后的图像,而我们后续的操作的基本前提都是金字塔已经指向拓充后的图像。
- 我们需要明白
-
计算图像特征点并进行均匀化
// Step 3 计算图像的特征点,并且将特征点进行均匀化。均匀的特征点可以提高位姿计算精度 // 存储所有的特征点,注意此处为二维的vector, // 第一维存储的是金字塔的层数,第二维存储的是那一层金字塔图像里提取的所有特征点 vector < vector<KeyPoint> > allKeypoints; //使用四叉树的方式计算每层图像的特征点并进行分配 ComputeKeyPointsOctTree(allKeypoints);
对这个函数我们也进行展开,里面针对每一层金字塔图像主要进行了以下三步:
-
按网格区域提取特征点,这里的几个图像边界参数可能让人迷惑,参考上面的那张图
//重新调整图像层数 allKeypoints.resize(nlevels); //图像cell的尺寸,是个正方形,可以理解为边长in像素坐标 const float W = 30; // 对每一层图像做处理 //遍历所有图像 for (int level = 0; level < nlevels; ++level) { //计算这层图像的坐标边界,也就是图中的绿色区域 const int minBorderX = EDGE_THRESHOLD-3; //这里的3是因为在计算FAST特征点的时候,需要建立一个半径为3的圆 const int minBorderY = minBorderX; //minY的计算就可以直接拷贝上面的计算结果了 const int maxBorderX = mvImagePyramid[level].cols-EDGE_THRESHOLD+3; const int maxBorderY = mvImagePyramid[level].rows-EDGE_THRESHOLD+3; //存储需要进行平均分配的特征点 vector<cv::KeyPoint> vToDistributeKeys; //一般地都是过量采集,所以这里预分配的空间大小是nfeatures*10 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; //计算当前网格最大的行坐标,这里的+6=+3+3,即考虑到了多出来3是为了cell边界像素进行FAST特征点提取用 //前面的EDGE_THRESHOLD指的应该是提取后的特征点所在的边界,所以minBorderY是考虑了计算半径时候的图像边界 //目测一个图像网格的大小是25*25啊 float maxY = iniY+hCell+6; //如果初始的行坐标就已经超过了有效的图像边界了,这里的“有效图像”是指原始的、可以提取FAST特征点的图像区域 if(iniY>=maxBorderY-3) //那么就跳过这一行 continue; //如果图像的大小导致不能够正好划分出来整齐的图像网格,那么就要委屈最后一行了 if(maxY>maxBorderY) maxY = maxBorderY; //开始列的遍历 for(int j=0; j<nCols; j++) { //计算初始的列坐标 const float iniX =minBorderX+j*wCell; //计算这列网格的最大列坐标,+6的含义和前面相同 float maxX = iniX+wCell+6; //判断坐标是否在图像中 //如果初始的列坐标就已经超过了有效的图像边界了,这里的“有效图像”是指原始的、可以提取FAST特征点的图像区域。 //并且应该同前面行坐标的边界对应,都为-3 //!BUG 正确应该是maxBorderX-3 if(iniX>=maxBorderX-6) continue; //如果最大坐标越界那么委屈一下 if(maxX>maxBorderX) maxX = maxBorderX; // FAST提取兴趣点, 自适应阈值 //这个向量存储这个cell中的特征点 vector<cv::KeyPoint> vKeysCell; //调用opencv的库函数来检测FAST角点 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++) { //NOTICE 到目前为止,这些角点的坐标都是基于图像cell的,现在我们要先将其恢复到当前的【坐标边界】下的坐标 //这样做是因为在下面使用八叉树法整理特征点的时候将会使用得到这个坐标 //在后面将会被继续转换成为在当前图层的扩充图像坐标系下的坐标 (*vit).pt.x+=j*wCell; (*vit).pt.y+=i*hCell; //然后将其加入到”等待被分配“的特征点容器中 vToDistributeKeys.push_back(*vit); }//遍历图像cell中的所有的提取出来的FAST角点,并且恢复其在整个金字塔当前层有效特征点拓充后图像下的坐标 }//当图像cell中检测到FAST角点的时候执行下面的语句 }//开始遍历图像cell的列 }//开始遍历图像cell的行 //声明一个对当前图层的特征点的容器的引用 vector<KeyPoint> & keypoints = allKeypoints[level]; //并且调整其大小为欲提取出来的特征点个数(当然这里也是扩大了的,因为不可能所有的特征点都是在这一个图层中提取出来的) keypoints.reserve(nfeatures);
-
使用四叉树对特征点进行均匀化
//声明一个对当前图层的特征点的容器的引用 vector<KeyPoint> & keypoints = allKeypoints[level]; //并且调整其大小为欲提取出来的特征点个数(当然这里也是扩大了的,因为不可能所有的特征点都是在这一个图层中提取出来的) keypoints.reserve(nfeatures); // 根据mnFeatuvector<KeyPoint> & keypoints = allKeypoints[level];resPerLevel,即该层的兴趣点数,对特征点进行剔除 //返回值是一个保存有特征点的vector容器,含有剔除后的保留下来的特征点 //得到的特征点的坐标,依旧是在当前图层下来讲的 keypoints = DistributeOctTree(vToDistributeKeys, //当前图层提取出来的特征点,也即是等待剔除的特征点 //NOTICE 注意此时特征点所使用的坐标都是在“半径扩充图像”下的 minBorderX, maxBorderX, //当前图层图像的边界,而这里的坐标却都是在“边缘扩充图像”下的 minBorderY, maxBorderY, mnFeaturesPerLevel[level], //希望保留下来的当前层图像的特征点个数 level); //当前层图像所在的图层
这里要注意半径扩充图像是前面图中的绿色区域,而边缘拓充图像是前面图中的浅灰色区域。
这个函数十分重要,我们先来看一下他的原理:
- 如果图片的宽度比较宽,就先把分成左右
w/h
份。一般的640×480的图像开始的时候只有一个初始node。然后将所有的特征点分配到初始node中。 - 如果node里面的点数>1,把每个node分成四个node,如果node里面的特征点为空,就直接删除掉这个节点。
- 新分的node的特征点数>1,就再分裂成4个node。如此,一直分裂。在其中某一步如果现存node的数目+可继续分的节点的数目*3大于需要的特征点数目,则进入优先度分裂模式,优先度的排序为node内包含的特征点数目。
- 终止条件为:node的总数量>需要的特征点数目,或者无法再进行分裂。
然后从每个node里面选择一个质量最好的FAST点。,具体的代码如下:
// Step 1 根据宽高比确定初始节点数目 //计算应该生成的初始节点个数,根节点的数量nIni是根据边界的宽高比值确定的,一般是1或者2 // ! bug: 如果宽高比小于0.5,nIni=0, 后面hx会报错 const int nIni = round(static_cast<float>(maxX-minX)/(maxY-minY)); //一个初始的节点的x方向有多少个像素 const float hX = static_cast<float>(maxX-minX)/nIni; //存储有提取器节点的链表 list<ExtractorNode> lNodes; //存储初始提取器节点指针的vector vector<ExtractorNode*> vpIniNodes; //重新设置其大小 vpIniNodes.resize(nIni); // Step 2 生成初始提取器节点 for(int i=0; i<nIni; i++) { //生成一个提取器节点 ExtractorNode ni; //设置提取器节点的图像边界 //注意这里和提取FAST角点区域相同,都是“半径扩充图像”,特征点坐标从0开始 ni.UL = cv::Point2i(hX*static_cast<float>(i),0); //UpLeft ni.UR = cv::Point2i(hX*static_cast<float>(i+1),0); //UpRight ni.BL = cv::Point2i(ni.UL.x,maxY-minY); //BottomLeft ni.BR = cv::Point2i(ni.UR.x,maxY-minY); //BottomRight //重设节点中储存的特征点的个数 ni.vKeys.reserve(vToDistributeKeys.size()); //将刚才生成的提取节点添加到链表中 //虽然这里的ni是局部变量,但是由于这里的push_back()是拷贝参数的内容到一个新的对象中然后再添加到列表中 //所以当本函数退出之后这里的内存不会成为“野指针” lNodes.push_back(ni); //存储这个初始的提取器节点句柄 vpIniNodes[i] = &lNodes.back(); } // Step 3 将特征点分配到子提取器节点中 for(size_t i=0;i<vToDistributeKeys.size();i++) { //获取这个特征点对象 const cv::KeyPoint &kp = vToDistributeKeys[i]; //按特征点的横轴位置,分配给属于那个图像区域的提取器节点(最初的提取器节点) vpIniNodes[kp.pt.x/hX]->vKeys.push_back(kp); } // Step 4 遍历此提取器节点列表,标记那些不可再分裂的节点,删除那些没有分配到特征点的节点 list<ExtractorNode>::iterator lit = lNodes.begin(); while(lit!=lNodes.end()) { //如果初始的提取器节点所分配到的特征点个数为1 if(lit->vKeys.size()==1) { //那么就标志位置位,表示此节点不可再分 lit->bNoMore=true; //更新迭代器 lit++; } ///如果一个提取器节点没有被分配到特征点,那么就从列表中直接删除它 else if(lit->vKeys.empty()) //注意,由于是直接删除了它,所以这里的迭代器没有必要更新;否则反而会造成跳过元素的情况 lit = lNodes.erase(lit); else //如果上面的这些情况和当前的特征点提取器节点无关,那么就只是更新迭代器 lit++; } //结束标志位清空 bool bFinish = false; //记录迭代次数,只是记录,并未起到作用 int iteration = 0; //声明一个vector用于存储节点的vSize和句柄对 //这个变量记录了在一次分裂循环中,那些可以再继续进行分裂的节点中包含的特征点数目和其句柄 vector<pair<int,ExtractorNode*> > vSizeAndPointerToNode; //调整大小,这里的意思是一个初始化节点将“分裂”成为四个 vSizeAndPointerToNode.reserve(lNodes.size()*4); // Step 5 利用四叉树方法对图像进行划分区域,均匀分配特征点 while(!bFinish) { //更新迭代次数计数器,只是记录,并未起到作用 iteration++; //保存当前节点个数,prev在这里理解为“保留”比较好 int prevSize = lNodes.size(); //重新定位迭代器指向列表头部 lit = lNodes.begin(); //需要展开的节点计数,这个一直保持累计,不清零 int nToExpand = 0; //因为是在循环中,前面的循环体中可能污染了这个变量,所以清空 //这个变量也只是统计了某一个循环中的点 //这个变量记录了在一次分裂循环中,那些可以再继续进行分裂的节点中包含的特征点数目和其句柄 vSizeAndPointerToNode.clear(); // 将目前的子区域进行划分 //开始遍历列表中所有的提取器节点,并进行分解或者保留 while(lit!=lNodes.end()) { //如果提取器节点只有一个特征点, if(lit->bNoMore) { // If node only contains one point do not subdivide and continue //那么就没有必要再进行细分了 lit++; //跳过当前节点,继续下一个 continue; } else { // If more than one point, subdivide //如果当前的提取器节点具有超过一个的特征点,那么就要进行继续分裂 ExtractorNode n1,n2,n3,n4; //再细分成四个子区域 lit->DivideNode(n1,n2,n3,n4); // Add childs if they contain points //如果这里分出来的子区域中有特征点,那么就将这个子区域的节点添加到提取器节点的列表中 //注意这里的条件是,有特征点即可 if(n1.vKeys.size()>0) { //注意这里也是添加到列表前面的 lNodes.push_front(n1); //再判断其中子提取器节点中的特征点数目是否大于1 if(n1.vKeys.size()>1) { //如果有超过一个的特征点,那么待展开的节点计数加1 nToExpand++; //保存这个特征点数目和节点指针的信息 vSizeAndPointerToNode.push_back(make_pair(n1.vKeys.size(),&lNodes.front())); //?这个访问用的句柄貌似并没有用到? // lNodes.front().lit 和前面的迭代的lit 不同,只是名字相同而已 // lNodes.front().lit是node结构体里的一个指针用来记录节点的位置 // 迭代的lit 是while循环里作者命名的遍历的指针名称 lNodes.front().lit = lNodes.begin(); } } //后面的操作都是相同的,这里不再赘述 if(n2.vKeys.size()>0) { lNodes.push_front(n2); if(n2.vKeys.size()>1) { nToExpand++; vSizeAndPointerToNode.push_back(make_pair(n2.vKeys.size(),&lNodes.front())); lNodes.front().lit = lNodes.begin(); } } if(n3.vKeys.size()>0) { lNodes.push_front(n3); if(n3.vKeys.size()>1) { nToExpand++; vSizeAndPointerToNode.push_back(make_pair(n3.vKeys.size(),&lNodes.front())); lNodes.front().lit = lNodes.begin(); } } if(n4.vKeys.size()>0) { lNodes.push_front(n4); if(n4.vKeys.size()>1) { nToExpand++; vSizeAndPointerToNode.push_back(make_pair(n4.vKeys.size(),&lNodes.front())); lNodes.front().lit = lNodes.begin(); } } //当这个母节点expand之后就从列表中删除它了,能够进行分裂操作说明至少有一个子节点的区域中特征点的数量是>1的 // 分裂方式是后加的节点先分裂,先加的后分裂 lit=lNodes.erase(lit); //继续下一次循环,其实这里加不加这句话的作用都是一样的 continue; }//判断当前遍历到的节点中是否有超过一个的特征点 }//遍历列表中的所有提取器节点 // Finish if there are more nodes than required features or all nodes contain just one point //停止这个过程的条件有两个,满足其中一个即可: //1、当前的节点数已经超过了要求的特征点数 //2、当前所有的节点中都只包含一个特征点 if((int)lNodes.size()>=N //判断是否超过了要求的特征点数 || (int)lNodes.size()==prevSize) //prevSize中保存的是分裂之前的节点个数,如果分裂之前和分裂之后的总节点个数一样,说明当前所有的 //节点区域中只有一个特征点,已经不能够再细分了 { //停止标志置位 bFinish = true; } // Step 6 当再划分之后所有的Node数大于要求数目时,就慢慢划分直到使其刚刚达到或者超过要求的特征点个数 //可以展开的子节点个数nToExpand x3,是因为一分四之后,会删除原来的主节点,所以乘以3 else if(((int)lNodes.size()+nToExpand*3)>N) { // //如果再分裂一次那么数目就要超了,这里想办法尽可能使其刚刚达到或者超过要求的特征点个数时就退出 //这里的nToExpand和vSizeAndPointerToNode不是一次循环对一次循环的关系,而是前者是累计计数,后者只保存某一个循环的 //一直循环,直到结束标志位被置位 while(!bFinish) { //获取当前的list中的节点个数 prevSize = lNodes.size(); //保留那些还可以分裂的节点的信息, 这里是深拷贝 vector<pair<int,ExtractorNode*> > vPrevSizeAndPointerToNode = vSizeAndPointerToNode; //清空 vSizeAndPointerToNode.clear(); // 对需要划分的节点进行排序,对pair对的第一个元素进行排序,默认是从小到大排序 // 优先分裂特征点多的节点,使得特征点密集的区域保留更少的特征点 //! 注意这里的排序规则非常重要!会导致每次最后产生的特征点都不一样。建议使用 stable_sort sort(vPrevSizeAndPointerToNode.begin(),vPrevSizeAndPointerToNode.end()); //遍历这个存储了pair对的vector,注意是从后往前遍历 for(int j=vPrevSizeAndPointerToNode.size()-1;j>=0;j--) { ExtractorNode n1,n2,n3,n4; //对每个需要进行分裂的节点进行分裂 vPrevSizeAndPointerToNode[j].second->DivideNode(n1,n2,n3,n4); // Add childs if they contain points //其实这里的节点可以说是二级子节点了,执行和前面一样的操作 if(n1.vKeys.size()>0) { lNodes.push_front(n1); if(n1.vKeys.size()>1) { //因为这里还有对于vSizeAndPointerToNode的操作,所以前面才会备份vSizeAndPointerToNode中的数据 //为可能的、后续的又一次for循环做准备 vSizeAndPointerToNode.push_back(make_pair(n1.vKeys.size(),&lNodes.front())); lNodes.front().lit = lNodes.begin(); } } if(n2.vKeys.size()>0) { lNodes.push_front(n2); if(n2.vKeys.size()>1) { vSizeAndPointerToNode.push_back(make_pair(n2.vKeys.size(),&lNodes.front())); lNodes.front().lit = lNodes.begin(); } } if(n3.vKeys.size()>0) { lNodes.push_front(n3); if(n3.vKeys.size()>1) { vSizeAndPointerToNode.push_back(make_pair(n3.vKeys.size(),&lNodes.front())); lNodes.front().lit = lNodes.begin(); } } if(n4.vKeys.size()>0) { lNodes.push_front(n4); if(n4.vKeys.size()>1) { vSizeAndPointerToNode.push_back(make_pair(n4.vKeys.size(),&lNodes.front())); lNodes.front().lit = lNodes.begin(); } } //删除母节点,在这里其实应该是一级子节点 lNodes.erase(vPrevSizeAndPointerToNode[j].second->lit); //判断是是否超过了需要的特征点数?是的话就退出,不是的话就继续这个分裂过程,直到刚刚达到或者超过要求的特征点个数 //作者的思想其实就是这样的,再分裂了一次之后判断下一次分裂是否会超过N,如果不是那么就放心大胆地全部进行分裂(因为少了一个判断因此 //其运算速度会稍微快一些),如果会那么就引导到这里进行最后一次分裂 if((int)lNodes.size()>=N) break; }//遍历vPrevSizeAndPointerToNode并对其中指定的node进行分裂,直到刚刚达到或者超过要求的特征点个数 //这里理想中应该是一个for循环就能够达成结束条件了,但是作者想的可能是,有些子节点所在的区域会没有特征点,因此很有可能一次for循环之后 //的数目还是不能够满足要求,所以还是需要判断结束条件并且再来一次 //判断是否达到了停止条件 if((int)lNodes.size()>=N || (int)lNodes.size()==prevSize) bFinish = true; }//一直进行nToExpand累加的节点分裂过程,直到分裂后的nodes数目刚刚达到或者超过要求的特征点数目 }//当本次分裂后达不到结束条件但是再进行一次完整的分裂之后就可以达到结束条件时 }// 根据兴趣点分布,利用4叉树方法对图像进行划分区域 // Step 7 保留每个区域响应值最大的一个兴趣点 //使用这个vector来存储我们感兴趣的特征点的过滤结果 vector<cv::KeyPoint> vResultKeys; //调整容器大小为要提取的特征点数目 vResultKeys.reserve(nfeatures); //遍历这个节点链表 for(list<ExtractorNode>::iterator lit=lNodes.begin(); lit!=lNodes.end(); lit++) { //得到这个节点区域中的特征点容器句柄 vector<cv::KeyPoint> &vNodeKeys = lit->vKeys; //得到指向第一个特征点的指针,后面作为最大响应值对应的关键点 cv::KeyPoint* pKP = &vNodeKeys[0]; //用第1个关键点响应值初始化最大响应值 float maxResponse = pKP->response; //开始遍历这个节点区域中的特征点容器中的特征点,注意是从1开始哟,0已经用过了 for(size_t k=1;k<vNodeKeys.size();k++) { //更新最大响应值 if(vNodeKeys[k].response>maxResponse) { //更新pKP指向具有最大响应值的keypoints pKP = &vNodeKeys[k]; maxResponse = vNodeKeys[k].response; } } //将这个节点区域中的响应值最大的特征点加入最终结果容器 vResultKeys.push_back(*pKP); } //返回最终结果容器,其中保存有分裂出来的区域中,我们最感兴趣、响应值最大的特征点 return vResultKeys;
注意这里的nToExpand的方式是有问题的,也应该每次清零,虽然作者后面的while循环解决了这个BUG可能带来的问题,但是这里逻辑是有问题的。
- 如果图片的宽度比较宽,就先把分成左右
-
对每一个特征点恢复坐标到边缘拓充图像坐标系上,也就是前面图中深灰色对应的坐标下,因为我们后续计算描述子需要在全尺寸的图像上进行
//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的横坐标边界
对这个函数进行展开,里面主要进行了一下步骤:
// 遍历所有的特征点 for (vector<KeyPoint>::iterator keypoint = keypoints.begin(), keypointEnd = keypoints.end(); keypoint != keypointEnd; ++keypoint) { // 调用IC_Angle 函数计算这个特征点的方向 keypoint->angle = IC_Angle(image, //特征点所在的图层的图像 keypoint->pt, //特征点在这张图像中的坐标 umax); //每个特征点所在图像区块的每行的边界 u_max 组成的vector }
其中的
IC_Angle
就是计算灰度重心和角度的关键函数static float IC_Angle(const Mat& image, // 要进行操作的某层金字塔图像 Point2f pt, // 当前特征点的坐标 const vector<int> & u_max) // 图像块的每一行的坐标边界 u_max { //图像的矩,前者是按照图像块的y坐标加权,后者是按照图像块的x坐标加权 int m_01 = 0, m_10 = 0; //获得这个特征点所在的图像块的中心点坐标灰度值的指针center const uchar* center = &image.at<uchar> (cvRound(pt.y), cvRound(pt.x)); // Treat the center line differently, v=0 //这条v=0中心线的计算需要特殊对待 //后面是以中心行为对称轴,成对遍历行数,所以PATCH_SIZE必须是奇数 for (int u = -HALF_PATCH_SIZE; u <= HALF_PATCH_SIZE; ++u) //注意这里的center下标u可以是负的!中心水平线上的像素按x坐标(也就是u坐标)加权 m_10 += u * center[u]; // Go line by line in the circular patch //这里的step1表示这个图像一行包含的字节总数。参考[https://blog.csdn.net/qianqing13579/article/details/45318279] int step = (int)image.step1(); //注意这里是以v=0中心线为对称轴,然后对称地每成对的两行之间进行遍历,这样处理加快了计算速度 for (int v = 1; v <= HALF_PATCH_SIZE; ++v) { // Proceed over the two lines //本来m_01应该是一列一列地计算的,但是由于对称以及坐标x,y正负的原因,可以一次计算两行 int v_sum = 0; // 获取某行像素横坐标的最大范围,注意这里的图像块是圆形的! int d = u_max[v]; //在坐标范围内挨个像素遍历,实际是一次遍历2个 // 假设每次处理的两个点坐标,中心线下方为(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) { //得到需要进行加运算和减运算的像素灰度值 //val_plus:在中心线下方x=u时的的像素灰度值 //val_minus:在中心线上方x=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坐标加权和(u坐标也有正负符号),相当于同时计算两行 m_10 += u * (val_plus + val_minus); } //将这一行上的和按照y坐标加权 m_01 += v * v_sum; } //为了加快速度还使用了fastAtan2()函数,输出为[0,360)角度,精度为0.3° return fastAtan2((float)m_01, (float)m_10); }
这里作者采用了一个技巧,利用
opencv
Mat类的step
来一次遍历关于u轴对称的两个点
-
为计算描述子分配空间
// Step 4 拷贝图像描述子到新的矩阵descriptors Mat descriptors; //统计整个图像金字塔中的特征点 int nkeypoints = 0; //开始遍历每层图像金字塔,并且累加每层的特征点个数 for (int level = 0; level < nlevels; ++level) nkeypoints += (int)allKeypoints[level].size(); //如果本图像金字塔中没有任何的特征点 if( nkeypoints == 0 ) //通过调用cv::mat类的.realse方法,强制清空矩阵的引用计数,这样就可以强制释放矩阵的数据了 //参考[https://blog.csdn.net/giantchen547792075/article/details/9107877] _descriptors.release(); else { //如果图像金字塔中有特征点,那么就创建这个存储描述子的矩阵,注意这个矩阵是存储整个图像金字塔中特征点的描述子的 _descriptors.create(nkeypoints, //矩阵的行数,对应为特征点的总个数 32, //矩阵的列数,对应为使用32*8=256位描述子 CV_8U); //矩阵元素的格式 //获取这个描述子的矩阵信息 // ?为什么不是直接在参数_descriptors上对矩阵内容进行修改,而是重新新建了一个变量,复制矩阵后,在这个新建变量的基础上进行修改? descriptors = _descriptors.getMat(); } //清空用作返回特征点提取结果的vector容器 _keypoints.clear(); //并预分配正确大小的空间 _keypoints.reserve(nkeypoints); ```
-
计算描述子
//因为遍历是一层一层进行的,但是描述子那个矩阵是存储整个图像金字塔中特征点的描述子,所以在这里设置了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; // preprocess the resized image // Step 5 对图像进行高斯模糊 // 深拷贝当前金字塔所在层级的图像 Mat workingMat = mvImagePyramid[level].clone(); // 注意:提取特征点的时候,使用的是清晰的原图像;这里计算描述子的时候,为了避免图像噪声的影响,使用了高斯模糊 GaussianBlur(workingMat, //源图像 workingMat, //输出图像 Size(7, 7), //高斯滤波器kernel大小,必须为正的奇数 2, //高斯滤波在x方向的标准差 2, //高斯滤波在y方向的标准差 BORDER_REFLECT_101);//边缘拓展点插值类型 // Compute the descriptors 计算描述子 // desc存储当前图层的描述子 Mat desc = descriptors.rowRange(offset, offset + nkeypointsLevel); // Step 6 计算高斯模糊后图像的描述子 computeDescriptors(workingMat, //高斯模糊之后的图层图像 keypoints, //当前图层中的特征点集合 desc, //存储计算之后的描述子 pattern); //随机采样模板 // 更新偏移量的值 offset += nkeypointsLevel; // Scale keypoint coordinates // Step 6 对非第0层图像中的特征点的坐标恢复到第0层图像(原图像)的坐标系下 // ? 得到所有层特征点在第0层里的坐标放到_keypoints里面 // 对于第0层的图像特征点,他们的坐标就不需要再进行恢复了 if (level != 0) { // 获取当前图层上的缩放系数 float scale = mvScaleFactor[level]; // 遍历本层所有的特征点 for (vector<KeyPoint>::iterator keypoint = keypoints.begin(), keypointEnd = keypoints.end(); keypoint != keypointEnd; ++keypoint) // 特征点本身直接乘缩放倍数就可以了 keypoint->pt *= scale; } // And add the keypoints to the output // 将keypoints中内容插入到_keypoints 的末尾 // keypoint其实是对allkeypoints中每层图像中特征点的引用,这样allkeypoints中的所有特征点在这里被转存到输出的_keypoints _keypoints.insert(_keypoints.end(), keypoints.begin(), keypoints.end()); }
在这里需要注意两个知识点:
-
作者在计算描述子时候对图像进行了高斯模糊,这样可以去除一些噪点的影响,高斯模糊函数的具体参数如下
其中各个参数的意义如下:- src - 输入图像;图像可以具有任意数量的通道,这些通道可以独立处理,但深度应为CV_8U,CV_16U,CV_16S,CV_32F或CV_64F。
- dst - 和输入图像类型和大小相同的输出图像
- ksize高斯内核大小。 ksize.width和ksize.height可以不同,但它们都必须为正数和奇数,也可以为零,然后根据sigma计算得出。
- sigmaX X方向上的高斯核标准偏差
- sigmaY Y方向上的高斯核标准差;如果sigmaY为零,则将其设置为等于sigmaX;如果两个sigma都为零,则分别从ksize.width和ksize.height计算得出
高斯模糊的效果如下:
-
计算描述子时候采用了位操作来加速,并且将描述子分成了32*8bit的格式,具体的针对每一个描述子的计算函数如下
/** * @brief 计算ORB特征点的描述子。注意这个是全局的静态函数,只能是在本文件内被调用 * @param[in] kpt 特征点对象 * @param[in] img 提取特征点的图像 * @param[in] pattern 预定义好的采样模板 * @param[out] desc 用作输出变量,保存计算好的描述子,维度为32*8 = 256 bit */ static void computeOrbDescriptor(const KeyPoint& kpt, const Mat& img, const Point* pattern, uchar* desc) { //得到特征点的角度,用弧度制表示。其中kpt.angle是角度制,范围为[0,360)度 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; //原始的BRIEF描述子没有方向不变性,通过加入关键点的方向来计算描述子,称之为Steer BRIEF,具有较好旋转不变特性 //具体地,在计算的时候需要将这里选取的采样模板中点的x轴方向旋转到特征点的方向。 //获得采样点中某个idx所对应的点的灰度值,这里旋转前坐标为(x,y), 旋转后坐标(x',y'),他们的变换关系: // x'= xcos(θ) - ysin(θ), y'= xsin(θ) + ycos(θ) // 下面表示 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)] //brief描述子由32*8位组成 //其中每一位是来自于两个像素点灰度的直接比较,所以每比较出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 }
在计算描述子之前需要先将
pattern
的点根据我们前面计算出来的灰度质心角度进行旋转从而保证特征点的旋转不变性,注意这里是主动旋转,和二维的位姿变换对应的被动旋转不同,二者之间差了一个转置。
-
-