ORB特征学习

原理

资料1
ORB一种改进的fast角点。
分为两部分：
（1）：FAST角点提取:找到图像的角点，计算特征点的主方向，为后续的BRIEF特征添加了旋转不变特性。
（2）：BRIEF描述子：对前一步特征提取出特征点的周围图像区域进行描述。ORB对BRIEF进行了改进，在BRIEF使用了先前计算的方向信息。

Fast 角点

主要检测局部像素灰度变化明显的地方，如果像素与邻域的像素差别较大（亮或暗），那么他更可能是角点，Fast角点相对于其他角点检测速度快。
步骤如下：

1. 选取像素$p$，假设他的灰度值为$I_p$
2. 设置一个阈值T(比如，$I_p$的20%)
3. 以像素p为中心，选取半径为3的圆上16个像素点。
4. 判断元上16个点，连续大于$I_p+T$或小于$I_p-T$的点的个数，当N=12,认为式特征点FAST-12，(当N取值为9或者11为FAST-9和FAST-11)
5. 循环上述4步循环整个像素
6. 极大值抑制保留相应极大值的角点，避免角点重复。

优缺点：
（1）：检测速度快
（2）：重复性不强，分布不均匀，不具有尺度和方向性。

图像金字塔

资料1
解决Fast角点不具有尺度不变性
对原始图像进行固定倍率的缩放，较小的看成远处看过来的场景。
高斯金字塔（Gaussian pyramid）

1. 图像的向下取样操作，即缩小图像
   得到的图像即为G_i+1的图像，显而易见，结果图像只有原图的四分之一。通过对输入图像G_i(原始图像)不停迭代以上步骤就会得到整个金字塔。同时我们也可以看到，向下取样会逐渐丢失图像的信息。以上就是对图像的向下取样操作，即缩小图像。
   （1）：对图像G_i进行高斯内核卷积，进行高斯模糊。
   （2）：将所有偶数行和列去除。

2. 对图像的向上取样，即放大图像
   （1）：将图像在每个方向扩大为原来的两倍，新增的行和列以0填充
   （2）：使用先前同样的内核(乘以4)与放大后的图像卷积，获得 “新增像素”的近似值

灰度质心法

矩的概念
以图像块灰度值作为权重的中心，解决图像FAST不有旋转不变性的缺点。
步骤如下：

1. 在图像块B中，定义图像块的据为： m p q = ∑ x , y ∈ B x p y q I ( x , y ) , p , q = { 0 , 1 } m_{pq}=\sum_{x,y\in B}x^py^qI(x,y),p,q=\{0,1\} mpq​=∑x,y∈B​xpyqI(x,y),p,q={0,1}
2. 通过矩找到图像块的中心：
   $$C=(\frac{m_{10}}{m_{00}},\frac{m_{01}}{M_{00}})$$
3. 连接几何中心O与灰度质心C，得到一个新的方向$\overrightarrow{OC}$,得到特征点方向为：
   $$\theta =arctan(\frac{m_{01}}{m_{10}})$$

BRIEF描述子

资料1
资料2
这是论文的较为核心部分。首先需要利用其它特征点提取算法，如FAST等生成特征点。随后，在特征点的领域一范围内$S\times S(S一般取31)$内随机选取点对，提供5种方式，如下：

1. x,y方向平均分布采样
2. x,y均服从$Gauss(0,S^2/25)$各向同性采样
3. x服从$Gauss(0,S^2/25)$，y服从$Gauss(0,S^2/100)$采样
4. x,y从网格中随机获取
5. x一直在(0,0)，y从网格中随机选取

生成描述子:
当随机生成点对后，通过点对之间灰度值比较来确定该点对对应的数值是0还是1
ROB算法：
使用随机选点比较，并利用FAST特征点阶段提取的关键点方向，计算旋转之后的"streer BRIEF"特征时ORB矩阵具有旋转不变性,即：
$D:为特征点在S\times S（S一般为31）邻域内点的集合$
$D_{\theta }:旋转后的BRIEF选择区域$
$\theta :灰度质心法计算的夹角$
$$D_{\theta }=R_{\theta }D$$

rBRIEF-改进特征点描述子的相关性

使用steeredBRIEF方法得到的特征描述子具有旋转不变性，但是却在另外一个性质上不如原始的BRIEF算法，即描述符的可区分性（相关性）。为了解决描述子的可区分性和相关性的问题，ORB论文中没有使用原始BRIEF算法中选取点对时的5种方法中的任意一种，而是使用统计学习的方法来重新选择点对集合。

1. 对每个特征点选取31x31领域，每个领域选择5x5的平均灰度值代替原来单个像素值进行比对，因此可以得到N=(31-5+1)x(31-5+1) = 729个可以比对的子窗口（patch），可以使用积分图像加快求取5x5邻域灰度平均值的速度。一共有M = 1+2+3+…+N = 265356种点对组合，也就是一个长度为M的01字符串。显然M远大于256，我们得筛选。
2. 重组所有点以及对应的初始二值串得到矩阵O，行数为提取得到的点数，每行是每个点对应的初始二值描述子
3. 对重组后的矩阵​O，按照每列均值与0.5的绝对差从小到大排序，得到矩阵T
4. 贪心选择：把T中第一列放进矩阵R（一开始为空）中，并从T中移除依次选择T的每列，与R中所有的列进行比较，如果相似度超过一定阈值，忽略，进行下一列，否则放进R中，并从T中移除重复以上过程直到选择​256个列，这样每个特征点就有256个0,1组成的描述子。如果不足256个，则降低阈值直到满足256就可，R即为最终特征描述矩阵。

特征点匹配

FLNN（近似最近邻算法）
资料2
创建KD树;KD最临近查找：

#include <iostream>
#include <Eigen/Core>
#include <vector>
#include <Eigen/Dense>
#include <opencv/cv.hpp>
using namespace std;
using namespace Eigen;
using namespace cv;
typedef struct KD_node//定义KD树节点
{
    
    cv::Point2d Node_piont;//每个点的数据为2维
    double Range[4];//节点代表的空间范围，0-1：x范围;2-3：y范围
    int split_number;//垂直于分割超平面的方向轴序号
    KD_node *left;//由位于该节点分割超平面左子空间内所有数据点所构成的k-d树
    KD_node *right;//由位于该节点分割超平面右子空间内所有数据点所构成的k-d树
    KD_node *parent;//父节点
}KDnode;
bool cmpx(Point2d &p1,Point2d &p2)//sort的回调函数
{
    return p1.x<p2.x;//从小到大排序
};
bool cmpy(Point2d &p1,Point2d &p2)//sort的回调函数
{
    return p1.y<p2.y;//从小到大排序
};
KDnode* Create_KD_trees(std::vector<Point2d> points,double Range[4], KD_node *parent)//初始化KD树
{
    const int init_tree=0;//初始化树
    const int sub_tree=1;
    const int compute_big_value=2;
    const int x_main_state=3;
    const int y_main_state=4;
    const int new_points=5;

    int now_state=init_tree;
    KDnode *now=(KD_node *)malloc(sizeof(KD_node));//当前节点结构体
    double left_range[4],right_range[4];//节点左边/右边点的集合
            double x_max,x_min,y_max,y_min;//初始化参数

    std::vector<double> x,y;//用作排序
    std::vector<Point2d> left_pionts,right_pionts;//左边点，右边点的集合
    switch (now_state)
    {
    case init_tree:{
        now->parent=parent;
         x_max=points[0].x;
         x_min=points[0].x;
         y_max=points[0].y;
         y_min=points[0].y;
        memcpy(left_range,Range,4*sizeof(double));
        memcpy(right_range,Range,4*sizeof(double));
        memcpy(now->Range,Range,4*sizeof(double));
        now_state=compute_big_value;
        
        }
    case compute_big_value:{
        int size=points.size();
        for(int i=0;i<points.size();i++){   //确定展开轴     
            if(x_max<points[i].x){x_max=points[i].x;};
            if(x_min>points[i].x){x_min=points[i].x;};
            if(y_max<points[i].y){y_max=points[i].y;};
            if(y_min>points[i].y){y_min=points[i].y;};
            };
        if(abs(x_max-x_min)>=(y_max-y_min)){
            sort(points.begin(),points.end(),cmpx);
            now_state=x_main_state; //x范围方差更大
            now->split_number=0;
            now->Node_piont.x=points[int(size/2)].x;
            now->Node_piont.y=points[int(size/2)].y;
            left_range[1]=points[int(size/2)].x;
            right_range[0]=points[int(size/2)].x;
            
            // cout<<"now->split_number"<<now->split_number<<endl;
            }   
        else{
            sort(points.begin(),points.end(),cmpy);
            now_state=y_main_state; //y范围方差更大
            now->Node_piont.x=points[int(size/2)].x;
            now->Node_piont.y=points[int(size/2)].y;
            now->split_number=1;
            left_range[3]=points[int(size/2)].y;
            right_range[2]=points[int(size/2)].y;
            // cout<<"now->split_number"<<now->split_number<<endl;
            }
            now_state=new_points;
        }
    
    case new_points:{
            int size=points.size();
            // cout<<"size:"<<size<<endl;
            if(points.size()>1){
                left_pionts.assign(points.begin(),points.begin()+int(size/2));
                right_pionts.assign(points.begin()+1+int(size/2),points.end());
                // cout<<"left_pionts:"<<left_pionts.size()<<endl;
                // cout<<"right_pionts:"<<right_pionts.size()<<endl;
                if(left_pionts.size()>0){now->left=Create_KD_trees(left_pionts,left_range,now);}
                if(right_pionts.size()>0){now->right=Create_KD_trees(right_pionts,right_range,now);}
            }
            else
            {

                now->left=NULL;now->right=NULL;
            }
            return now;
        }
    }



    
};

KDnode* KD_tree_search(KDnode* KD,Point2d target,std::vector<KDnode*> &line){//返回目标值所在子数
   Point2d kd_node=KD->Node_piont;
//    cout<<"Node_piont:"<<kd_node.x<<","<<kd_node.y<<endl;
//    cout<<"KD->split_number"<<KD->split_number<<endl;
   line.push_back(KD);
   if((kd_node.x==target.x)&&(kd_node.y==target.y)){
       return KD;
   }
   switch (KD->split_number)
   {
   case 0:{//x列分割
        if((kd_node.x)>=(target.x)){
            if(KD->left==NULL){return KD;}
            KD_tree_search(KD->left,target,line);
        }
        else{
            if(KD->right==NULL){return KD;}
            KD_tree_search(KD->right,target,line);     
        }
           
        break;
   }
    case 1:{//y列分割
        if((kd_node.y)>=(target.y)){
            if(KD->left==NULL){return KD;}
            KD_tree_search(KD->left,target,line);
        }
        else{
            if(KD->right==NULL){return KD;}
            KD_tree_search(KD->right,target,line); 
        }      
        break;
   }
   default:
       break;
   }
};
bool is_boundary(KDnode *KD,double distance){//判断是否有边界交点
    double left_diff=KD->Node_piont.x-distance;//左边界
    double right_diff=KD->Node_piont.x+distance;
    double up_diff=KD->Node_piont.y+distance;
    double down_diff=KD->Node_piont.y-distance;
    if((left_diff>KD->Range[0])&&(right_diff<KD->Range[1])&&(down_diff>KD->Range[2])&&(up_diff<KD->Range[3]))
    {
        return 1;//满足条件直接返回
    }
    return 0;
}
KDnode* traback(Point2d target,double distance,std::vector<KDnode*> &line){
   
    
    const int ordinary=1;//新的距离大于最近点
    const int update_distance=2;//新的距离小于最近点
    const int update_routes=3;
    const int next_point=4;
    static int now_state=ordinary;
    double R;//当前点与目标点的欧氏距离
    KDnode *KD;
    KDnode *nearest;
    while (line.size()!=0)
    {
        cout<<"line.size():"<<line.size()<<endl;
        cout<<"distance"<<distance<<endl;
        switch (now_state)
        {
        case ordinary:{
            cout<<"1:"<<endl;
            KD=line.back();  //得到末尾元素
            line.pop_back();       //删除末尾元素
            R=sqrtf((KD->Node_piont.x-target.x)*(KD->Node_piont.x-target.x)+(KD->Node_piont.y-target.y)*(KD->Node_piont.y-target.y));
            cout<<"R:"<<R<<endl;
            if(R==0){return nearest;}
            if(KD->parent==NULL){return nearest;}
            if(R<=distance){
                    now_state=update_distance;
                    nearest=KD;
                }
            break;
            }
        case update_distance:{
            cout<<"2:"<<endl;
            distance=R;//更新距离值
            if(is_boundary(KD,distance)){//以当前点为原点，原点目标距离为半径划圆，查看是否超出边界
                return nearest;
            }
            now_state=update_routes;
            break;
        }
        case update_routes:{//重新计算路线
            cout<<"3:"<<endl;
                switch (KD->split_number)//根据划分的主维重新筛选
            {
                case 0:{//x列分割
                    cout<<"4:"<<endl;
                    if((KD->Node_piont.x)<(target.x)){       //这里与原始规律相反，防止路线还是原来值
                        if(KD->left!=NULL)
                        KD_tree_search(KD->left,target,line);
                    }
                    else{
                        if(KD->right!=NULL)
                        KD_tree_search(KD->right,target,line);     
                    }            
                    break;
                }
                case 1:{//y列分割
                    cout<<"5:"<<endl;
                    if((KD->Node_piont.y)<(target.y)){             //这里与原始规律相反，防止路线还是原来值        
                        if(KD->left!=NULL)
                        KD_tree_search(KD->left,target,line);
                    }
                    else{
                        if(KD->right!=NULL)
                        KD_tree_search(KD->right,target,line); 
                    }      
                    break;
                }
                default:
                    break;
            }
            now_state=ordinary;
            break;
        }
        default:
            break;
        }
    }
    
}
int main(int argc, char const *argv[])
{
    std::vector<Point2d> points;
    Point2d target_point=Point2d(2.2,4.5);
    points.push_back(Point2d(2,3));
    points.push_back(Point2d(5,4));
    points.push_back(Point2d(9,6));
    points.push_back(Point2d(4,7));
    points.push_back(Point2d(8,1));
    points.push_back(Point2d(7,2));
    double Range[4]={0,10,0,10};
    KDnode* KD=Create_KD_trees(points,Range,NULL);
    std::vector<KDnode*> line;
    KDnode* kd_sereach=KD_tree_search(KD,target_point,line);
    Point2d search_point=line.back()->Node_piont;
    double distance=sqrtf((search_point.x-target_point.x)*(search_point.x-target_point.x)+(search_point.y-target_point.y)*(search_point.y-target_point.y));
    KDnode* nearest=traback(target_point,distance,line);
    cout<<line.size()<<endl;
    cout<<"nearest:"<<nearest->Node_piont.x<<","<<nearest->Node_piont.y<<endl;

    
}

ORB实现

opencv keypoint类
Feature2D类

1. 进行特征匹配的一般步骤
   实例化特征点检测器，进行特征点检测
   实例化描述子提取器，对计算得到特征点提取描述子
   实例化匹配器，根据描述子进行匹配
   筛选优秀匹配结果并绘图

2. 调库实现

#include <iostream>
#include <opencv/cv.hpp>
#include <chrono>

using namespace std;
using namespace cv;

bool cmp( DMatch &m1,  DMatch &m2){
    return (m1.distance<m2.distance);
};
int main(int argc, const char** argv) {
    //读取图像
    Mat img1=imread("/home/n1/notes/ORB/ORB/1.png",CV_LOAD_IMAGE_COLOR);
    Mat img2=imread("/home/n1/notes/ORB/ORB/2.png",CV_LOAD_IMAGE_COLOR);
    assert(img1.data!=nullptr&&img2.data!=nullptr);

    //初始化
    vector<KeyPoint> key_point_1,key_point_2;//关键点的信息
    //opencv keypoint类KeyPoint::KeyPoint(pt(x,y),, float _size, float _angle=-1, float _response=0, int _octave=0, int _class_id=-1)
    //point2f(x,y) 坐标
    //_size关键点的邻域直径大小
    //_angle 关键点方向大小
    //_response 为角点的可能性
    //_octave 金字塔第几层获得的角点
    //_class_id 对角点进行分类，未设定时为-1
    Mat description_1,description_2;
    Ptr<FeatureDetector> detetor=ORB::create();//特征提取器
    //Ptr<> opencv提供的智能类模板指针
    //ORB::create（nfeatures = 500, 最大关键点数
            //  scaleFactor = 1.2, 金字塔抽取率，scaleFactor = 2表示经典金字塔
            //  nlevels = 8,        金字塔层数
            //  edgeThreshold = 31, 未检测到要素的边框大小 应与patchSize大小一致
            //  firstLevel = 0,  定义金字塔第一层的序号
            //  WTA_K = 2,          生成的BRIEF描述符时，选择同时几个元素进行比较。可以设置为2,3和4，其中2为默认值。
            //  scoreType = HARRIS_SCORE,     对角点的好坏进行排序，默认的HARRIS_SCORE，还可以设置为FAST_SCORE，速度快精度差
            //  patchSize = 31,     滑动窗口区
            //  fastThreshold = 20）
    Ptr<DescriptorExtractor> descriptor=ORB::create();//特征描述子
    Ptr<DescriptorMatcher> matcher=DescriptorMatcher::create("BruteForce-Hamming");//特征匹配"BruteForce-Hamming"是暴力匹配方法

    //第一步 Oriented Fast角点检测
    chrono::steady_clock::time_point t1=chrono::steady_clock::now();//获取当前时间
    detetor->detect(img1,key_point_1);//提取特征点
    detetor->detect(img2,key_point_2);

    //第二步 计算BRIEF描述子
    descriptor->compute(img1,key_point_1,description_1);
    descriptor->compute(img2,key_point_2,description_2);
    chrono::steady_clock::time_point t2=chrono::steady_clock::now();//获取当前时间
    chrono::duration<double> time_used = chrono::duration_cast<chrono::duration<double>>(t2 - t1);
    cout << "extract ORB cost = " << time_used.count() << " seconds. " << endl;
    Mat output1;
    drawKeypoints(img1,key_point_1,output1,Scalar::all(-1),DrawMatchesFlags::DEFAULT);
    //Scalar::all(-1):随机像素值
    //cv2.DRAW_MATCHES_FLAGS_DEFAULT：创建输出图像矩阵，使用现存的输出图像绘制匹配对和特征点，对每一个关键点只绘制中间点
    //cv2.DRAW_MATCHES_FLAGS_DRAW_OVER_OUTIMG：不创建输出图像矩阵，而是在输出图像上绘制匹配对
    //cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS：对每一个特征点绘制带大小和方向的关键点图形
    //cv2.DRAW_MATCHES_FLAGS_NOT_DRAW_SINGLE_POINTS：单点的特征点不被绘制
    imshow("ORB features",output1);

    //第三步 BRIEF 描述子进行匹配
    vector<DMatch> matches;
    //CV_WRAP DMatch( int _queryIdx, int _trainIdx, int _imgIdx,float _distance )
    //DMatch主要用来储存匹配信息的结构体，
    //_queryIdx是要匹配的描述子id，
    //_trainIdx是被匹配的描述子id，
    //_imgIdx匹配图像的id
    //_distance距离
    t1 = chrono::steady_clock::now();
    matcher->match(description_1,description_2,matches);
    //description_1:query
    //description_2:trainI
    //输出的信息
    time_used = chrono::duration_cast<chrono::duration<double>>(t2 - t1);
    cout << "match ORB cost = " << time_used.count() << " seconds. " << endl;
    
    //第4步：匹配点筛选
    auto min_max=minmax_element(matches.begin(),matches.end(),cmp);
    //返回值：make_pair类型
    double min=min_max.first->distance;
    double max=min_max.second->distance;
    cout<<"-- min dist :"<<min<<endl;
    cout<<"-- max dist :"<<max<<endl;
    //简单的数据筛选
    //当描述子之间的距离大于两倍的最小距离时,即认为匹配有误.但有时候最小距离会非常小,设置一个经验值30作为下限.
    std::vector<DMatch> good_matches;
    for (int i = 0; i < description_1.rows; i++) {
        if (matches[i].distance<=std::max(2 * min, 30.0)) {
        good_matches.push_back(matches[i]);
        }
    }

    //第5步：显示效果图
    Mat img_match;
    Mat img_goodmatch;
    drawMatches(img1, key_point_1, img2, key_point_2, matches, img_match);
    //void drawMatches（img1，key_point_1,img2,key_point_2,matches,img_match,matchColor,singlePointColor,matchesMask,flags)
    // img1:源图像1
    // keypoints1: 源图像1的特征点.
    // img2: 源图像2.
    // keypoints2: 源图像2的特征点
    // matches1to2: 源图像1的特征点匹配源图像2的特征点[matches[i]] .
    // outImg: 输出图像具体由flags决定.
    // matchColor: 匹配的颜色（特征点和连线),若matchColor==Scalar::all(-1)，颜色随机.
    // singlePointColor: 单个点的颜色，即未配对的特征点，若matchColor==Scalar::all(-1)，颜色随机.
    // matchesMask: Mask决定哪些点将被画出，若为空，则画出所有匹配点.
    // flags:0:默认值创建输出图像;输出图像可以重用;将绘制两个源图像、匹配项和单个关键点;对于每个关键点，仅绘制中心点(带关键点尺寸和方向的圆圈)
    //       1:不创建输出图像矩阵;在输出图像上输出图像
    //       2：不绘制单个关键点
    //       4：对于每个关键点，围绕关键点的圆圈以及方向
    drawMatches(img1, key_point_1, img2, key_point_2, good_matches, img_goodmatch);
    imshow("all matches", img_match);
    imshow("good matches", img_goodmatch);
    waitKey(0);
    return 0;
}

3. 手写实现

• List item