ORB特征点

Oriented FAST and Rotated BRIEF

1 Harris角点

因为在ORB中要用到Harris的一个结论，那就顺便复习一下Harris吧！

1.1 基本原理

人眼对角点的识别通常是在一个局部的小区域或小窗口完成的。如果在各个方向上移动这个特征的小窗口，窗口内区域的灰度发生了较大的变化，那么就认为在窗口内遇到了角点。如果这个特定的窗口在图像各个方向上移动时，窗口内图像的灰度没有发生变化，那么窗口内就不存在角点；如果窗口在某一个方向移动时，窗口内图像的灰度发生了较大的变化，而在另一些方向上没有发生变化，那么，窗口内的图像可能就是一条直线的线段。

1.2 推导

对于图像$I(x,y)$，当在点$(x,y)$处平移$(\Delta x,\Delta y)$后的自相似性，可以通过自相关函数给出：

$$c(x,y;\Delta x,\Delta y) = \sum_{(u,v)\in W(x,y)}w(u,v)(I(u,v) I(u+\Delta x,v+\Delta y))^2$$
其中， $W(x,y)$是以点 $(x,y)$为中心的窗口， $w(u,v)$为加权函数，它既可是常数，也可以是高斯加权函数。
根据泰勒展开，对图像 $I(x,y)$在平移 $(\Delta x,\Delta y)$后进行一阶近似：
$$I(u+\Delta x,v+\Delta y) = I(u,v)+I_x(u,v)\Delta x+I_y(u,v)\Delta y+O(\Delta x^2,\Delta y^2) \approx I(u,v)+I_x(u,v)\Delta x+I_y(u,v)\Delta y$$
其中， $I_x,I_y$是图像 $I(x,y)$的偏导数，则自相关函数可以近似为：
$$c(x,y;\Delta x,\Delta y)\approx \sum_w (I_x(u,v)\Delta x+I_y(u,v)\Delta y)^2=[\Delta x,\Delta y]M(x,y)\begin{bmatrix}\Delta x \\ \Delta y\end{bmatrix}$$
其中：
$$M(x,y)=\sum_w \begin{bmatrix}I_x(x,y)^2&I_x(x,y)I_y(x,y) \\ I_x(x,y)I_y(x,y)&I_y(x,y)^2\end{bmatrix} = \begin{bmatrix}\sum_w I_x(x,y)^2&\sum_w I_x(x,y)I_y(x,y) \\\sum_w I_x(x,y)I_y(x,y)&\sum_w I_y(x,y)^2\end{bmatrix}=\begin{bmatrix}A&C\\C&B\end{bmatrix}$$
也就是说图像 $I(x,y)$在点 $(x,y)$处平移 $(\Delta x,\Delta y)$后的自相关函数可以近似为二项函数：
$$c(x,y;\Delta x,\Delta y)\approx A\Delta x^2+2C\Delta x\Delta y+B\Delta y^2$$
其中：
$$A=\sum_w I_x^2, B=\sum_w I_y^2,C=\sum_w I_x I_y$$
二次项函数本质上就是一个椭圆函数。椭圆的扁率和尺寸是由 $M(x,y)$的特征值 $\lambda_1、\lambda_2$决定的，椭圆的方向是由 $M(x,y)$的特征矢量决定的，如下图所示，椭圆方程为：
$$[\Delta x,\Delta y]M(x,y)\begin{bmatrix}\Delta x \\ \Delta y\end{bmatrix} = 1$$

椭圆函数特征值与图像中的角点、直线（边缘）和平面之间的关系如下图所示。共可分为三种情况：

• 图像中的直线。一个特征值大，另一个特征值小，$\lambda_1\gg \lambda_2$或$\lambda_2\gg \lambda_1$。自相关函数值在某一方向上大，在其他方向上小。
• 图像中的平面。两个特征值都小，且近似相等；自相关函数数值在各个方向上都小。
• 图像中的角点。两个特征值都大，且近似相等，自相关函数在所有方向都增大。

根据二次项函数特征值的计算公式，我们可以求$M(x,y)$矩阵的特征值。但是Harris给出的角点差别方法并不需要计算具体的特征值，而是计算一个角点响应值$R$来判断角点。$R$的计算公式为：

$$R=det \boldsymbol{M} - \alpha(trace\boldsymbol{M})^2$$
式中， $det\boldsymbol{M}$为矩阵 $\boldsymbol{M}=\begin{bmatrix}A&B\\B&C\end{bmatrix}$的行列式； $trace\boldsymbol{M}$为矩阵 $\boldsymbol{M}$的直迹； $\alpha$为经常常数，取值范围为0.04~0.06。事实上，特征是隐含在 $det\boldsymbol{M}$和 $trace\boldsymbol{M}$中，因为， $$det\boldsymbol{M} = \lambda_1\lambda_2=AC-B^2$$

$$trace\boldsymbol{M}=\lambda_2+\lambda_2=A+C$$

2 ORB特征点

ORB特征是将FAST特征点的检测方法与BRIEF特征描述子结合起来，并在它们原来的基础上做了改进与优化。首先，它利用FAST特征点检测的方法来检测特征点，然后利用Harris角点的度量方法，从FAST特征点从挑选出Harris角点响应值最大的$N$个特征点。其中Harris角点的响应函数定义为：

$$R=det \boldsymbol{M} - \alpha(trace\boldsymbol{M})^2$$

2.1 旋转不变性

我们知道FAST特征点是没有尺度不变性的，所以我们可以通过构建高斯金字塔，然后在每一层金字塔图像上检测角点，来实现尺度不变性。那么，对于局部不变性，我们还差一个问题没有解决，就是FAST特征点不具有方向，ORB的论文中提出了一种利用灰度质心法来解决这个问题，灰度质心法假设角点的灰度与质心之间存在一个偏移，这个向量可以用于表示一个方向。对于任意一个特征点$p$来说，我们定义$p$的邻域像素的矩为：

$$m_{pq} = \sum_{x,y \in B}x^py^qI(x,y)$$

其中$I(x,y)$为点$(x,y)$处的灰度值，$B$ 表示一个图像块。那么我们可以得到图像的质心为：

$$C = \left(\frac{m_{10}}{m_{00}},\frac{m_{01}}{m_{00}}\right)$$

那么特征点与质心的夹角定义为FAST特征点的方向：

$$\theta = arctan(m_{01},m_{10}) \tag{2.1}$$

2.2 描述子

ORB选择了BRIEF作为特征描述方法，但是我们知道BRIEF是没有旋转不变性的，所以我们需要给BRIEF加上旋转不变性，把这种方法称为“Steer BREIF”。对于任何一个特征点来说，它的BRIEF描述子是一个长度为$n$的二值码串，这个二值串是由特征点周围$n$个点对（$2n$个点）生成的，现在我们将这$2n$个点$(x_i,y_i),i = 1,2,\cdots,2n$组成一个矩阵$S$

$$S=\begin{pmatrix}x_1&x_2&\cdots&x_{2n}\\y_1&y_2&\cdots&y_{2n}\end{pmatrix}$$

ORB中采用了一个更有效的方法：使用邻域方向$\theta$和对应的旋转矩阵$R_{\theta}$，构建$S$的一个校正版本$S_{\theta}$

$$S_{\theta} = R_{\theta}S$$

其中

$$R_{\theta} = \begin{bmatrix}cos\theta & sin\theta \\ –sin\theta &cos\theta\end{bmatrix}$$
而 $\theta$即（2.1式）中求得的主方向。
最后将 $S_\theta$中相邻的两个点的灰度值比较，如果前者大于后者则为1，反之为0. 这样就可以得到一个长度为 $n$的二值码串，可以使用汉明距离来度量这些二值码串之间的相似性。

2.3 OpenCV3中的ORB

OpenCV3中使用create函数创建ORB对象

static Ptr<ORB> create(int nfeatures=500, float scaleFactor=1.2f, int nlevels=8, int edgeThreshold=31,int firstLevel=0, int WTA_K=2, int scoreType=ORB::HARRIS_SCORE, int patchSize=31, int fastThreshold=20);

• nfeatures - 最多提取的特征点的数量；
• scaleFactor - 金字塔图像之间的尺度参数，类似于SIFT中的k；
• nlevels – 高斯金字塔的层数；
• edgeThreshold – 边缘阈值，这个值主要是根据后面的patchSize来定的，靠近边缘edgeThreshold以内的像素是不检测特征点的。
• firstLevel - 看过SIFT都知道，我们可以指定第一层的索引值，这里默认为0。
• WTA_K - 用于产生BIREF描述子的 点对的个数，一般为2个，也可以设置为3个或4个，那么这时候描述子之间的距离计算就不能用汉明距离了，而是应该用一个变种。OpenCV中，如果设置WET_K = 2，则选用点对就只有2个点，匹配的时候距离参数选择NORM_HAMMING，如果WET_K设置为3或4，则BIREF描述子会选择3个或4个点，那么后面匹配的时候应该选择的距离参数为NORM_HAMMING2。
• scoreType - 用于对特征点进行排序的算法，你可以选择HARRIS_SCORE，也可以选择FAST_SCORE，但是它也只是比前者快一点点而已。
• patchSize – 用于计算BIREF描述子的特征点邻域大小。
• fastThreshold - FAST角点检测中阈值，若记此阈值为T，假如选取的圆上连续有N个点的亮度大于$I_p+T$或者小于 $I_p+T$，那么像素p可以被认为是特征点。

  #include <opencv2/core/core.hpp> 
  #include <opencv2/highgui/highgui.hpp> 
  #include <opencv2/imgproc/imgproc.hpp> 
  #include <opencv2/features2d/features2d.hpp>
  #include <iostream>
  #include <string>
  
  using namespace cv;
  using namespace std;
  
  int main(int argc, char** argv) 
  {
      if(argc != 3)
      {
          cout << "Usage: orbdetect image1 image2" << endl;
          return 1;
      }
  
      string img1(argv[1]);
      string img2(argv[2]);
  
      Mat img_1 = imread(img1); 
      Mat img_2 = imread(img2);
  
      // -- Step 1: Detect the keypoints using STAR Detector 
      std::vector<KeyPoint> keypoints_1,keypoints_2; 
      Ptr<ORB> orb = ORB::create(300);
      orb->detect(img_1, keypoints_1); 
      orb->detect(img_2, keypoints_2);
  
      // -- Stpe 2: Calculate descriptors (feature vectors) 
      Mat descriptors_1, descriptors_2; 
      orb->compute(img_1, keypoints_1, descriptors_1); 
      orb->compute(img_2, keypoints_2, descriptors_2);
  
      //-- Step 3: Matching descriptor vectors with a brute force matcher 
      BFMatcher matcher(NORM_HAMMING); 
      std::vector<DMatch> mathces; 
      matcher.match(descriptors_1, descriptors_2, mathces); 
      // -- dwaw matches 
      Mat img_mathes; 
      drawMatches(img_1, keypoints_1, img_2, keypoints_2, mathces, img_mathes); 
      // -- show 
      imshow("Mathces", img_mathes);
  
      waitKey(0); 
      return 0; 
  }

  参考资料

  1. ORB特征点
  2. Harris角点
  3. 高翔. 视觉SLAM十四讲从理论到实践