ORB2单目读代码笔记5--利用灰度质心计算ORB特征点方向，实现旋转不变性

5.利用灰度质心计算ORB特征点方向，实现旋转不变性

ComputeKeyPointOctTree 跳转 computeOrientation
computeOrientation 跳转 IC_Angle
补 : 灰度质心法
IC_Angle 结束返回 computeOrientation
computeOrientation 结束返回ComputeKeyPointOctTree
ComputeKeyPointOctTree 结束返回operator()

均匀化之后

scaledPatchSize 表示当前金字塔层缩放倍数
nkps 表示保留的特征点数目

将保留的特征点统一成当前金字塔层的坐标,即加上minBorderX,Y

        //PATCH_SIZE是对于底层的初始图像来说的，现在要根据当前图层的尺度缩放倍数进行缩放得到缩放后的PATCH大小 和特征点的方向计算有关const int scaledPatchSize = PATCH_SIZE*mvScaleFactor[level];// Add border to coordinates and scale information//获取剔除过程后保留下来的特征点数目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;}}

使用computeOrientation计算所有层的特征点的方向

    // compute orientations//然后计算这些特征点的方向信息，注意这里还是分层计算的for (int level = 0; level < nlevels; ++level)computeOrientation(mvImagePyramid[level],   //对应的图层的图像allKeypoints[level],     //这个图层中提取并保留下来的特征点容器umax);                   //以及PATCH的横坐标边界
}

ComputeKeyPointOctTree 跳转 computeOrientation

这个函数实际只起到了遍历所有特征点的作用,在每个特征点都调用IC_Angle函数来计算特征点方向

static void computeOrientation(const Mat& image, vector<KeyPoint>& keypoints, const vector<int>& umax)
{// 遍历所有的特征点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}
}

computeOrientation 跳转 IC_Angle

基本原理就是灰度质心法:每个特征点以几何中心P为圆心,画圆,计算出灰度质心Q,则P -> Q就是该特征点的方向.后面比较特征点时,先把方向统一,然后比较特征点描述子,这样特征点就不会因为旋转而不能匹配了

补 : 灰度质心法

目的:使特征点具有旋转不变性

灰度质心:像素集中点

原理:即使图像旋转,图像几何中心与灰度质心的连线也是随之改变的,可以通过这个方向恢复原始图像

输入参数是:

image输入图像
pt当前特征点坐标
u_max是指定一行后的u轴边界坐标

static float IC_Angle(const Mat& image, Point2f pt,  const vector<int> & u_max)
{

定义的变量center 表示特征点坐标的指针

    //图像的矩，前者是按照图像块的y坐标加权，后者是按照图像块的x坐标加权int m_01 = 0, m_10 = 0;//获得这个特征点所在的图像块的中心点坐标灰度值的指针centerconst uchar* center = &image.at<uchar> (cvRound(pt.y), cvRound(pt.x));

遍历中心线(-R到R),求m10和m01,要乘加权center表示该点灰度值

不需要求m01,因为当前遍历的是u轴,v坐标为0

step 是一行的字节总数,后面遍历的时候增量乘step表示像素增量

   // 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从1开始,到半径的长度

    //注意这里是以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];

在每个v下,遍历u,前面定义过u_max[v],表示给定一个v坐标,求出最大u坐标,这里赋值给了d

遍历d相当于遍历了给定v之后的当前行

由两个公式m_10 m_01 需要当前u或v值乘以当前像素灰度值

val_plus相当于I(x,y),val_minus相当于I(x,-y)

有了这两个变量,下面的程序就很好懂了

fastAtan2为计算特征点方向的函数

        //在坐标范围内挨个像素遍历，实际是一次遍历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);
}

IC_Angle 结束返回 computeOrientation

computeOrientation 结束返回ComputeKeyPointOctTree

ComputeKeyPointOctTree 结束返回operator()