以下内容均来自对计算机视觉life公众号中资料的整理,作为自己的笔记使用

ORB-SLAM2中特征点的处理

1 ORB特征点简介

ORB特征由关键点和描述子组成

FAST关键点:
1.选取像素 p p p,假设其亮度为 I p I_p Ip​;
2.设置一个阈值 T T T(比如 I p I_p Ip​的20%);
3.以像素 p p p为中心，选取半径为3的圆上的16个像素点；
4.假如选取的圆上，有连续的N个点的亮度大于 I p + T I_p+T Ip​+T或 小 于 I p − T 小于I_p-T 小于Ip​−T，那么像素p可以被认为是特征点；
5.循环以上4步，对每一个像素执行相同操作。

Brief描述子:

BRIEF是一种二进制描述子,描述量由许多个0和1组成.这些0和1编码了关键点附近两个随机像素(如p和q)的大小关系,如果p比q大则取1,反之取0.

BRIEF算法的核心思想是在关键点P的周围以一定模式选取N个点对，把这N个点对的比较结果组合起来作为描述子。为了保持踩点固定，工程上采用特殊设计的固定的pattern来做.

2.灰度质心法

原始的FAST关键点没有方向信息,如果图像发生旋转,brief描述子也会发生变化,使得特征点对旋转不鲁棒,所以通过灰度质心法计算出特征点的方向.

灰度质心:可以理解为,图像区域内,以一像素灰度值作为权重的中心,形心指向质心的向量即为关键点的主方向.
或者说的不规范一点,灰度质心可以理解为:这一区域灰度值分布的中心点

当图像发生旋转的时候,质心一定也会跟着转,这样我们构建的坐标系也随着质心一起转,从而保证了旋转不变性.

计算灰度质心的步骤:

• 1.我们定义区域图像的矩为:

m p q = ∑ x , y x p y q I ( x , y ) , p , q = { 0 , 1 } m_{pq}=\sum_{x,y}x^py^qI(x,y),\quad \quad \quad p,q=\{0,1\} mpq​=x,y∑​xpyqI(x,y),p,q={0,1}

其中,p,q取值为0或1, I ( x , y ) I(x,y) I(x,y)表示像素在坐标 ( x , y ) (x,y) (x,y)处图像的灰度值; m p q m_{pq} mpq​表示图像的矩.
在半径为R的图形图像区域,沿两个坐标轴x,y方向的图像的矩为:
m 10 = ∑ x = − R R ∑ y = − R R x I ( x , y ) m 01 = ∑ x = − R R ∑ y = − R R y I ( x , y ) m_{10}=\sum_{x=-R}^R\sum_{y=-R}^RxI(x,y)\\ m_{01}=\sum_{x=-R}^R\sum_{y=-R}^RyI(x,y) m10​=x=−R∑R​y=−R∑R​xI(x,y)m01​=x=−R∑R​y=−R∑R​yI(x,y)
上两个式子也可以理解为区域内图像的灰度值在分别 x x x方向和 y y y方向的加和.然后下式 m 00 m_{00} m00​是圆内所有像素的灰度值的加和:
m 00 = ∑ x = − R R ∑ y = − R R I ( x , y ) m_{00}=\sum_{x=-R}^R\sum_{y=-R}^RI(x,y) m00​=x=−R∑R​y=−R∑R​I(x,y)

• 2.然后分别将x,y方向的加与总和做比值,得到图像的质心:

C = ( c x , c y ) = ( m 10 m 00 , m 01 m 00 ) C=(c_x,c_y)=\left( \frac{m_{10}}{m_{00}},\frac{m_{01}}{m_{00}} \right) C=(cx​,cy​)=(m00​m10​​,m00​m01​​)

• 3.然后关键点的"主方向"可以表示为从圆形图像形心O指向质心C的方向向量 O C → \overrightarrow{OC} OC ,于是关键点的旋转角度记为:
  θ = a r c t a n 2 ( c y , c x ) = a r c t a n 2 ( m 01 , m 10 ) \theta=arctan2(c_y,c_x)=arctan2(m_{01},m_{10}) θ=arctan2(cy​,cx​)=arctan2(m01​,m10​)

在ORB-SLAM2中的实现为,先计算四分之一圆的下半部分的横坐标(u轴)边界,然后根据圆对称性计算出圆的上半部分的边界

具体代码如下:

首先确定圆的边界:

ORBextractor::ORBextractor( ...  ):	nfeatures(_nfeatures), ... ,  minThFAST(_minThFAST)
{
	...
    ...
        
	//预先计算圆形patch中行的结束位置,+1中的1表示那个圆的中间行,HALF_PATCH_SIZE 即半径,取值为15
    umax.resize(HALF_PATCH_SIZE + 1); 
    int v, v0,		//循环辅助变量,辅助变量
		vmax = cvFloor(HALF_PATCH_SIZE * sqrt(2.f) / 2 + 1);	//计算圆的最大行号，+1应该是把中间行也给考虑进去了
    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));		//结果都是大于0的结果，表示x坐标在这一行的边界

    // Make sure we are symmetric
	//这里其实是使用了对称的方式计算上四分之一的圆周上的umax，目的也是为了保持严格的对称（如果按照常规的想法做，由于cvRound就会很容易出现不对称的情况，
	//同时这些随机采样的特征点集也不能够满足旋转之后的采样不变性了）
	for (v = HALF_PATCH_SIZE, v0 = 0; v >= vmin; --v)
    {
        while (umax[v0] == umax[v0 + 1])
            ++v0;
        umax[v] = v0;
        ++v0;
    }
}

然后计算方向:

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
    }
}

IC_Angle()函数如下:

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)
		//注意这里的center下标u可以是负的！中心水平线上的像素按x坐标（也就是u坐标）加权
        m_10 += u * center[u];

	//这里的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)
    {
		//本来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);
}

3.特征提取与分配

这部分是也是ORB-SLAM2非常核心也是非常关键的地方

4.特征点的分配(每层分配策略)

ORB-SLAM2中特征点分配的策略为:(比如一张图像一共要提取1000个特征点)
根据输入图像构建图像金字塔 -> 计算出整个图像金字塔的总面积S -> 计算出各层金字塔图像面积占总面积的百分比 -> 根据各层占比分配特征点提取数量

具体计算方法如下:

假设需要提取的特征点数目为 N N N,金字塔共 m m m层,第0层图像的宽为 W W W,高为 H H H,对应的面积为 C = H ⋅ W C=H·W C=H⋅W,图像金字塔的缩放因子为 s , 0 < S < 1 s,0<S<1 s,0<S<1,在ORB-SLAM2中, m = 8 , s = 1 1.2 m=8,s=\frac{1}{1.2} m=8,s=1.21​.

那么整个图像金字塔的总面积是:
S = H ⋅ W ⋅ ( s 2 ) 0 + H ⋅ W ⋅ ( s 2 ) 1 + ⋅ ⋅ ⋅ + H ⋅ W ⋅ ( s 2 ) m − 1 = H W 1 − ( s 2 ) m 1 − s 2 = C 1 − ( s 2 ) m 1 − s 2 \begin{aligned} S&=H·W·(s^2)^0+H·W·(s^2)^1+···+H·W·(s^2)^{m-1}\\ &=HW\frac{1-(s^2)^m}{1-s^2}\\ &=C\frac{1-(s^2)^m}{1-s^2}\\ \end{aligned} S​=H⋅W⋅(s2)0+H⋅W⋅(s2)1+⋅⋅⋅+H⋅W⋅(s2)m−1=HW1−s21−(s2)m​=C1−s21−(s2)m​​
单位面积应该分配的特征点数量为:
N a v g = N S = N C 1 − ( s 2 ) m 1 − s 2 = N ( 1 − s 2 ) C ( 1 − ( s 2 ) m ) \begin{aligned} N_{avg}&=\frac{N}{S}\\ &=\frac{N}{C\frac{1-(s^2)^m}{1-s^2}}\\ &=\frac{N(1-s^2)}{C(1-(s^2)^m)} \end{aligned} Navg​​=SN​=C1−s21−(s2)m​N​=C(1−(s2)m)N(1−s2)​​
第0层应该分配的特征点数量为:
N 0 = N a v g ⋅ C = N ( 1 − s 2 ) 1 − ( s 2 ) m \begin{aligned} N_0&=N_{avg}·C\\ &=\frac{N(1-s^2)}{1-{(s^2)^m}} \end{aligned} N0​​=Navg​⋅C=1−(s2)mN(1−s2)​​
第i层应该分配的特征点数量为:
N i = N ( 1 − s 2 ) C ( 1 − ( s 2 ) m ) C ( s 2 ) i = N ( 1 − s 2 ) 1 − ( s 2 ) m ( s 2 ) i \begin{aligned} N_i&=\frac{N(1-s^2)}{C(1-(s^2)^m)}C(s^2)i\\ &=\frac{N(1-s^2)}{1-(s^2)^m}(s^2)^i \end{aligned} Ni​​=C(1−(s2)m)N(1−s2)​C(s2)i=1−(s2)mN(1−s2)​(s2)i​
在ORB-SLAM2 的代码里，不是按照面积均摊的，而是按照面积的开方来均摊特征点的，也就是将上述公式中的 s 2 s^2 s2换成 s s s即可。

ORBextractor::ORBextractor( ... ) :		...
{
	//存储每层图像缩放系数的vector调整为符合图层数目的大小
    mvScaleFactor.resize(nlevels);  
	//存储这个+子的平方
    mvLevelSigma2.resize(nlevels);
	//对于初始图像，这两个参数都是1
    mvScaleFactor[0]=1.0f;
    mvLevelSigma2[0]=1.0f;
	//然后逐层计算图像金字塔中图像相当于初始图像的缩放系数 
    for(int i=1; i<nlevels; i++)  
    {
		//其实就是这样累乘计算得出来的
        mvScaleFactor[i]=mvScaleFactor[i-1]*scaleFactor;
		//原来这里的sigma^2就是每层图像相对于初始图像缩放因子的平方
        mvLevelSigma2[i]=mvScaleFactor[i]*mvScaleFactor[i];
    }

    //接下来的两个向量保存上面的参数的倒数
    mvInvScaleFactor.resize(nlevels);
    mvInvLevelSigma2.resize(nlevels);
    for(int i=0; i<nlevels; i++)
    {
        mvInvScaleFactor[i]=1.0f/mvScaleFactor[i];
        mvInvLevelSigma2[i]=1.0f/mvLevelSigma2[i];
    }

    //调整图像金字塔vector以使得其符合设定的图像层数
    mvImagePyramid.resize(nlevels);

	//每层需要提取出来的特征点个数，这个向量也要根据图像金字塔设定的层数进行调整
    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++ )
    {
		//分配 cvRound : 返回个参数最接近的整数值
        mnFeaturesPerLevel[level] = cvRound(nDesiredFeaturesPerScale);
		//累计
        sumFeatures += mnFeaturesPerLevel[level];
		//乘系数
        nDesiredFeaturesPerScale *= factor;
    }
    //由于前面的特征点个数取整操作，可能会导致剩余一些特征点个数没有被分配，所以这里就将这个余出来的特征点分配到最高的图层中
    mnFeaturesPerLevel[nlevels-1] = std::max(nfeatures - sumFeatures, 0);
}

未完待续…