Mat srcImg1 = imread("basketball1.png", 1);
Mat srcImg2 = imread("basketball2.png", 1);
if (srcImg1.empty() || srcImg2.empty())
    return -1;

Ptr<cv::Feature2D> detector = xfeatures2d::SIFT::create(128);
std::vector<KeyPoint> keypoints1, keypoints2;
Mat descriptors1, descriptors2;
// 提取图片的特征点及特征点描述
detector->detect(srcImg1, keypoints1);
detector->compute(srcImg1, keypoints1, descriptors1);
rootSift(descriptors1);

detector->detect(srcImg2, keypoints2);
detector->compute(srcImg2, keypoints2, descriptors2);
rootSift(descriptors2);

void rootSift(Mat &descriptors, const float eps = 1e-7)
{
    // Compute sums for L1 Norm
    Mat sums_vec;
    descriptors = abs(descriptors); //otherwise we draw sqrt of negative vals
    reduce(descriptors, sums_vec, 1 /*sum over columns*/, CV_REDUCE_SUM, CV_32FC1);
    for (int row = 0; row < descriptors.rows; row++) {
        int offset = row * descriptors.cols;
        for (int col = 0; col < descriptors.cols; col++) {
            descriptors.at<float>(offset + col) = sqrt(descriptors.at<float>(offset + col) /
                (sums_vec.at<float>(row) + eps) /*L1-Normalize*/);
        }
        // L2 distance
        normalize(descriptors.row(row), descriptors.row(row), 1.0, 0.0, NORM_L2);

    }
    return;
}

size_t maxCount = descriptors1.rows;
int M = 32;
int efConstruction = 100;// TestFindSelf，140以下有不能查找到自己的情况
int vecdim = descriptors1.cols;// 特征维度
hnswlib::InnerProductSpace *L2Space = new hnswlib::InnerProductSpace(vecdim);
hnswlib::HierarchicalNSW<float> *hnsw = new hnswlib::HierarchicalNSW<float>(L2Space, maxCount, M, efConstruction);
int count = 0;
for (size_t i = 0; i < descriptors1.rows; i++)
{
    float *mass = new float[vecdim];
    for (size_t j = 0; j < vecdim; j++)
    {
        mass[j] = descriptors1.at<float>(i, j);
    }
    hnsw->addPoint((void*)mass, count);
    ++count;
    delete[] mass;
}

int nn = 3;
std::vector<DMatch> matchPoints;
Mat indices, dists;
for (size_t i = 0; i < descriptors2.rows; i++)
{
    float *mass = new float[vecdim];
    for (size_t j = 0; j < vecdim; j++)
    {
        mass[j] = descriptors2.at<float>(i, j);
    }
    std::priority_queue<std::pair<float, hnswlib::labeltype>> candidates = hnsw->searchKnn((void*)mass, nn);
    delete[] mass;

    for (int k = 0; k < nn; k++)
    {
        double dAngle = abs(keypoints2.at(i).angle - keypoints1.at(candidates.top().second).angle);
        if (abs(candidates.top().first) < 0.1 && dAngle < 10)
        {
            std::cout << i << " angle diff: " << dAngle << ", dist: " << candidates.top().first << " "<< candidates.top().second << std::endl;
            DMatch dmatches(candidates.top().second, i, candidates.top().first);
            matchPoints.push_back(dmatches);
        }
        candidates.pop();
    }
}