OpenCV+OpenGL 双目立体视觉三维重建

0.绪论

这篇文章主要为了研究双目立体视觉的最终目标——三维重建，系统的介绍了三维重建的整体步骤。双目立体视觉的整体流程包括：图像获取，摄像机标定，特征提取（稠密匹配中这一步可以省略），立体匹配，三维重建。我在做双目立体视觉问题时，主要关注的点是立体匹配，本文主要关注最后一个步骤三维重建中的：三角剖分和纹理贴图以及对应的OpenCV+OpenGL代码实现。

1.视差计算

1.1基于视差信息的三维重建

特征提取
由双目立体视觉进行三位重建的第一步是立体匹配，通过寻找两幅图像中的对应点获取视差。OpenCV 中的features2d库中包含了很多常用的算法，其中特征点定位的算法有FAST, SIFT, SURF ,MSER, HARRIS等，特征点描述算法有SURF, SIFT等，还有若干种特征点匹配算法。这三个步骤的算法可以任选其一，自由组合，非常方便。经过实验，选择了一种速度、特征点数量和精度都比较好的组合方案：FAST角点检测算法+SURF特征描述子+FLANN(Fast Library for Approximate Nearest Neighbors) 匹配算法。

在匹配过程中需要有一些措施来过滤误匹配。一种比较常用的方法是比较第一匹配结果和第二匹配结果的得分差距是否足够大，这种方法可以过滤掉一些由于相似造成的误匹配。还有一种方法是利用已经找到的匹配点，使用RANSAC算法求得两幅视图之间的单应矩阵，然后将左视图中的坐标P用单应矩阵映射到右视图的Q点，观察与匹配结果Q’的欧氏距离是否足够小。当然由于图像是具有深度的，Q与Q’必定会有差距，因此距离阈值可以设置的稍微宽松一些。我使用了这两种过滤方法。

另外，由于图像有些部分的纹理较多，有些地方则没有什么纹理，造成特征点疏密分布不均匀，影响最终重建的效果，因此我还采取了一个措施：限制特征点不能取的太密。如果新加入的特征点与已有的某一特征点距离太小，就舍弃之。最终匹配结果如下图所示，精度和均匀程度都较好。

代码：

// choose the corresponding points in the stereo images for 3d reconstruction
void GetPair( Mat &imgL, Mat &imgR, vector<Point2f> &ptsL, vector<Point2f> &ptsR ) 
{
    Mat descriptorsL, descriptorsR;
    double tt = (double)getTickCount();

   Ptr<FeatureDetector> detector = FeatureDetector::create( DETECTOR_TYPE ); // factory mode
    vector<KeyPoint> keypointsL, keypointsR; 
    detector->detect( imgL, keypointsL );
    detector->detect( imgR, keypointsR );

    Ptr<DescriptorExtractor> de = DescriptorExtractor::create( DESCRIPTOR_TYPE );
    //SurfDescriptorExtractor de(4,2,true);
    de->compute( imgL, keypointsL, descriptorsL );
    de->compute( imgR, keypointsR, descriptorsR );

    tt = ((double)getTickCount() - tt)/getTickFrequency(); // 620*555 pic, about 2s for SURF, 120s for SIFT

    Ptr<DescriptorMatcher> matcher = DescriptorMatcher::create( MATCHER_TYPE );
    vector<vector<DMatch>> matches;
    matcher->knnMatch( descriptorsL, descriptorsR, matches, 2 ); // L:query, R:train

    vector<DMatch> passedMatches; // save for drawing
    DMatch m1, m2;
    vector<Point2f> ptsRtemp, ptsLtemp;
    for( size_t i = 0; i < matches.size(); i++ )
    {
        m1 = matches[i][0];
        m2 = matches[i][1];
        if (m1.distance < MAXM_FILTER_TH * m2.distance)
        {
            ptsRtemp.push_back(keypointsR[m1.trainIdx].pt);
            ptsLtemp.push_back(keypointsL[i].pt);
            passedMatches.push_back(m1);
        }
    }

    Mat HLR;
    HLR = findHomography( Mat(ptsLtemp), Mat(ptsRtemp), CV_RANSAC, 3 );
    cout<<"Homography:"<<endl<<HLR<<endl;
    Mat ptsLt; 
    perspectiveTransform(Mat(ptsLtemp), ptsLt, HLR);

    vector<char> matchesMask( passedMatches.size(), 0 );
    int cnt = 0;
    for( size_t i1 = 0; i1 < ptsLtemp.size(); i1++ )
    {
        Point2f prjPtR = ptsLt.at<Point2f>((int)i1,0); // prjx = ptsLt.at<float>((int)i1,0), prjy = ptsLt.at<float>((int)i1,1);
         // inlier
        if( abs(ptsRtemp[i1].x - prjPtR.x) < HOMO_FILTER_TH &&
            abs(ptsRtemp[i1].y - prjPtR.y) < 2) // restriction on y is more strict
        {
            vector<Point2f>::iterator iter = ptsL.begin();
            for (;iter!=ptsL.end();iter++)
            {
                Point2f diff = *iter - ptsLtemp[i1];
                float dist = abs(diff.x)+abs(diff.y);
                if (dist < NEAR_FILTER_TH) break;
            }
            if (iter != ptsL.end()) continue;

            ptsL.push_back(ptsLtemp[i1]);
            ptsR.push_back(ptsRtemp[i1]);
            cnt++;
            if (cnt%1 == 0) matchesMask[i1] = 1; // don't want to draw to many matches
        }
    }

    Mat outImg;
    drawMatches(imgL, keypointsL, imgR, keypointsR, passedMatches, outImg, 
        Scalar::all(-1), Scalar::all(-1), matchesMask, DrawMatchesFlags::NOT_DRAW_SINGLE_POINTS);
    char title[50];
    sprintf_s(title, 50,