手写ORB图像特征（附代码）

ORB特征由key points和descriptor两部分组成，也就是找出代表性的点，而且给出一个向量来表示点周围像素的信息。

它的keypoint称为"Oriented FAST"，是一种改进的FAST角点（FAST不在这里介绍）
它的descriptor称为BRIFF，是一种二进制的描述方式，在后面的代码里有体现。

ORB特征提取有以下2步：

1. FAST角点提取，相比于原版FAST，它加入了特征点的主方向，也就是方向信息，为BRIFF描述子提供了旋转不变性。
2. BRIFF描述子，就是在前一步key point的周围提取有描述信息的向量。

第一步的FAST角点提取不在这里介绍，会直接调用opencv。

第二步简单介绍下BRIFF描述子：
它是一种二进制描述子，就是说特征向量里面只有0和1，怎么确定是0还是1？
这时取key point（计算FAST角点得到）附近的随机的两个像素点p和q，如果p<q就是1，反之是0。只取2个点肯定不足以表达特征吧，于是代码中对每个key point附近随机取256个点，就得到256维的0，1向量。
由于是二进制，所以用一个整型就能储存。

再说旋转不变性，这里在一个小的图像块中找到质心，图像块的几何中心O与质心C会形成一个向量，这个向量就代表了一个方向角theta。
还记得旋转变换么，
$\left(\begin{array}{c}{x}^{{}^{\prime }}\\ y^{{}^{\prime }}\end{array}\right)=\left(\begin{array}{cc}cos\theta & -sin\theta \\ sin\theta & cos\theta \end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)$

对每个随机的256个p,q做旋转，加到keypoint上（p,q是以keypoint为中心的(x,y), 可理解成keypoint的随机偏移向量），这样会得到旋转后的256个随机点形成的01特征向量，这就是旋转之后的"Steer BRIFF"特征，具有较好的旋转不变性。

下面在代码里体现以上步骤

//输入参数：图像，FAST角点，保存描述子的vector
void computeORB(const Mat &img, vector<KeyPoint> &keypoints, vector<DescType> &descriptors) {
    const int half_patch_size = 8;
    const int half_boundary = 16; //设的一圈边界
    int bad_points = 0;
    for(auto &kp : keypoints) { //auto &获取集合中数据的引用，不会copy数据，for中修改数据会反映到原数据
        //keypoint中的pt指(x,y)坐标
        if(kp.pt.x < half_boundary || kp.pt.y < half_boundary ||
        kp.pt.x >= img.cols - half_boundary || kp.pt.y >= img.rows - half_boundary) {
            //outside
            bad_points ++;
            descriptors.push_back({});
            continue;
        }
        float m01 = 0, m10 = 0;  //质心
        //以keypoint为中心
        for(int dx = -half_patch_size; dx < half_patch_size; dx++) {
            for(int dy = -half_patch_size; dy < half_patch_size; dy++) {
                uchar pixel = img.at<uchar>(kp.pt.y + dy, kp.pt.x + dx);
                m01 += dx * pixel;
                m10 += dy * pixel;
            }
        }

        //angle = arc tan(m01/m10);
        float m_sqrt = sqrt(m01 * m01 + m10 * m10);
        float sin_theta = m01 / m_sqrt;
        float cos_theta = m10 / m_sqrt;

        //描述子8x256
        DescType desc(8, 0);
        for(int i = 0; i < 8; i++) {
            uint32_t d = 0;
            for(int k = 0; k < 32; k++) {
                int idx_pq = i*32 + k;
                //随机选两个偏移量p, q, 其中p,q为点(x,y)格式，为可重复性，定义在pattern
                Point2f p(ORB_pattern[idx_pq * 4], ORB_pattern[idx_pq * 4 + 1]);
                Point2f q(ORB_pattern[idx_pq * 4 + 2], ORB_pattern[idx_pq * 4 + 3]);

                //对偏移量p,q做旋转变换，旋转角度为theta，加到keypoint上得到随机的两点pp, qq
                Point2f pp = Point2f(cos_theta * p.x - sin_theta * p.y, sin_theta * p.x + cos_theta * p.y) + kp.pt;
                Point2f qq = Point2f(cos_theta * q.x - sin_theta * q.y, sin_theta * q.x + cos_theta * q.y) + kp.pt;
                //比较pp, qq处的像素大小，决定0 or 1
                if(img.at<uchar>(pp.y, pp.x) < img.at<uchar>(qq.y, qq.x)) {
                    d |= 1 << k;  //第k位设为1
                }
            }
            desc[i] = d;  //每个i 对应32位的0，1特征
        }
        descriptors.push_back(desc);
    }
    cout << "bad/total: " << bad_points << "/" << keypoints.size() << endl;
}

描述子，也就是特征向量，会保存在descriptor里面，如果输入两幅图，会得到两个descriptor，下面对两个descriptor中的特征做matching

void BfMatch(
    const vector<DescType> &desc1, const vector<DescType> &desc2, vector<DMatch> &matches
        ) {
    const int d_max = 40;

    for(size_t i1 = 0; i1 < desc1.size(); i1++) {
        if(desc1[i1].empty()) continue;
        //每个i1对应一个最小距离的i2
        //DMatch包含queryIdx(特征1的idx),trainIdx(与特征1相匹配的特征2的idx),distance(特征1与特征2的距离）
        DMatch m{(int)i1, 0 ,256};
        for(size_t i2 = 0; i2 < desc2.size(); i2++) {
            if(desc2[i2].empty()) continue;
            int distance = 0;
            for(int k = 0; k < 8; k++) {    //特征为8x32
                //统计两个distance中不同bit的个数
                distance += _mm_popcnt_u32((unsigned int)(desc1[i1][k] ^ desc2[i2][k]));
            }
            //找出最小距离和对应的desc2的index
            if(distance < d_max && distance < m.distance) {
                m.distance = distance;
                m.trainIdx = i2;
            }
            if(m.distance < d_max) {
                matches.push_back(m);
            }

        }
    }
}

测试代码：

void ORB_test() {
    Mat first_image = imread("../imgs/1.png", 0); //灰度图
    Mat second_image = imread("../imgs/2.png", 0);
    assert(first_image.data != nullptr && second_image.data != nullptr);

    //detect Fast keypoints using thres=40
    vector<KeyPoint> keypoints1;
    FAST(first_image, keypoints1, 40);
    vector<DescType> descriptor1;
    //计算ORB特征
    computeORB(first_image, keypoints1, descriptor1);

    //same for the second
    vector<KeyPoint> keypoints2;
    vector<DescType> descriptor2;
    FAST(second_image, keypoints2, 40);
    computeORB(second_image, keypoints2, descriptor2);

    //find matches
    vector<DMatch> matches;
    BfMatch(descriptor1, descriptor2, matches);

    //plot matches
    Mat image_show;
    drawMatches(first_image, keypoints1, second_image, keypoints2, matches, image_show);
    imshow("matches", image_show);
    waitKey();
}

参考链接