精华内容
下载资源
问答
  • 双目视觉三维重建

    2018-08-14 12:45:37
    一些双目三维重建代码,有matlab和c++的,效果不错。
  • opencv写的双目视觉摄像机标定和三维重建代码opencv写的双目视觉摄像机标定和三维重建代码opencv写的双目视觉摄像机标定和三维重建代码
  • vs2015+opencv2.4.10实现的双目立体视觉三维重建c++代码。SGBM立体匹配
  • OpenCV+OpenGL 双目立体视觉三维重建

    万次阅读 多人点赞 2016-08-08 00:02:47
    我在做双目立体视觉问题时,主要关注的点是立体匹配,本文主要关注最后一个步骤三维重建中的:三角剖分和纹理贴图以及对应的OpenCV+OpenGL代码实现。1.视差计算1.1基于视差信息的三维重建特征提

    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, "%.3f s, %d matches, %d passed", tt, matches.size(), cnt);
        imshow(title, outImg);
        waitKey();
    }

    p.s. 源代码中的基于特征点的视差计算有点问题,还在调试中,希望有经验的大牛共同解决一下。


    最新回复,特别鸣谢大神:G3fire(update 20180718)

    代码在opencv2.4.9版本下运行的,由于SIFT和SURF的专利约束需要nofree的引用.
    
    在Reconstuction3d.cpp中添加initModule_nonfree();
    同样在head.h中添加
    #pragma comment(lib,"opencv_nonfree249d.lib"),
    把Algorithm g_algo 改成= FEATURE_PT
    就可以运行基于特征点的视差计算了 
    
    楼24的修改特征点检测的创建方法没有运行通

    特别感谢24楼的回复

    博主的特征点匹配这边运行时会崩溃,我用VS2013+opencv2.4.10版本,然后修改特征点检测的创建方法就可以用了。

    例如:

    SurfFeatureDetector detector;
    detector.detect(imgL1, keypointsL); 
    detector.detect(imgR1, keypointsR);

    1.2基于块匹配的视差计算

    上面提取特征点的过程中实际上忽略了一个辅助信息:对应点应当是取在对应极线上的一个区间内的。利用这个信息可以大幅简化对应点的匹配,事实上只要用L1距离对一个像素周围的block计算匹配距离就可以了,也就是OpenCV中实现的块匹配算法的基本思路。比起特征点匹配,这是一种“稠密”的匹配算法,精度也可以接受。下图中浅色表示视差较大,对应深度较浅。左侧有一块区域是左右视图不相交的部分,因此无法计算视差。
    这里写图片描述
    可以发现视差计算结果中有很多噪声。事实上在纹理平滑的区域,还有左右视图中不同遮挡的区域,是很难计算视差的。因此我利用最近邻插值和数学形态学平滑的方法对视差图进行了修复(见cvFuncs2.cpp中的FixDisparity函数):
    这里写图片描述

    // roughly smooth the glitches on the disparity map
    void FixDisparity( Mat_<float> & disp, int numberOfDisparities ) 
    {
        Mat_<float> disp1;
        float lastPixel = 10;
        float minDisparity = 23;// algorithm parameters that can be modified
        for (int i = 0; i < disp.rows; i++)
        {
            for (int j = numberOfDisparities; j < disp.cols; j++)
            {
                if (disp(i,j) <= minDisparity) disp(i,j) = lastPixel;
                else lastPixel = disp(i,j);
            }
        }
         int an = 4;    // algorithm parameters that can be modified
        copyMakeBorder(disp, disp1, an,an,an,an, BORDER_REPLICATE);
        Mat element = getStructuringElement(MORPH_ELLIPSE, Size(an*2+1, an*2+1));
        morphologyEx(disp1, disp1, CV_MOP_OPEN, element);
        morphologyEx(disp1, disp1, CV_MOP_CLOSE, element);
        disp = disp1(Range(an, disp.rows-an), Range(an, disp.cols-an)).clone();
    }

    对应点的选取
    上面提到,为了获得较好的重构效果,特征点最好取在深度变化较大的区域。基于这种猜想,我首先对上面的视差图求梯度,然后找到梯度最大的点,观察梯度的方向,如果是偏x方向,就在该点左右若干像素各取一个点;否则就在上下若干像素各取一个点。然后根据这两个点的视差值就可以计算出另外一个视图中的对应点的坐标。特征点还不能分布过密,因此我取完一对特征点后,将其周围一圈像素的梯度置零,然后在寻找下一个梯度最大值,这样一直下去,直到取够特征点数。
    特征点也不能全取在深度变化剧烈的区域,在平坦的区域也可以取一些。最终我取的特征点如下图:
    这里写图片描述
    其中紫色的点是在较平坦的区域取到的,其他颜色是在边界区域取到的。这些算法实现在ChooseKeyPointsBM函数中。

    2.计算世界坐标

    一般双目立体视觉中使用的实验图像都是经过外极线矫正的,计算3D坐标也比较方便,其实利用外极线约束(以及其他的约束条件)可以极大的降低立体匹配的计算量。见下图:
    这里写图片描述
    如果(x1,y1),(x2,y2)用各自图像上的像素坐标表示,L和(X,Y,Z)用毫米表示,f用像素表示的话,用相似三角形的知识就可以推出:
    这里写图片描述
    其中W和H是图像的宽高(像素数),y是y1和y2的均值,Z加负号是为了保持右手坐标系,而Y加负号是由于图像成像过程中上下发生了倒转。三维世界原点取为左摄像机的焦点。计算的代码见cvFunc.cpp中的StereoTo3D函数。

     // calculate 3d coordinates.
    // for rectified stereos: pointLeft.y == pointRight.y
    // the origin for both image is the top-left corner of the left image.
    // the x-axis points to the right and the y-axis points downward on the image.
    // the origin for the 3d real world is the optical center of the left camera
    // object -> optical center -> image, the z value decreases.
    
    void StereoTo3D( vector<Point2f> ptsL, vector<Point2f> ptsR, vector<Point3f> &pts3D,
                    float focalLenInPixel, float baselineInMM, Mat img,
                    Point3f &center3D, Vec3f &size3D) // output variable, the center coordinate and the size of the object described by pts3D
    {
        vector<Point2f>::iterator iterL = ptsL.begin(),
            iterR = ptsR.begin();
    
        float xl, xr, ylr;
        float imgH = float(img.rows), imgW = float(img.cols);
        Point3f pt3D;
        float minX = 1e9, maxX = -1e9;
        float minY = 1e9, maxY = -1e9;
        float minZ = 1e9, maxZ = -1e9;
    
        Mat imgShow = img.clone();
        char str[100];
        int ptCnt = ptsL.size(), showPtNum = 30, cnt = 0;
        int showIntv = max(ptCnt/showPtNum, 1);
        for ( ; iterL != ptsL.end(); iterL++, iterR++)
        {
            xl = iterL->x;
            xr = iterR->x; // need not add baseline
            ylr = (iterL->y + iterR->y)/2;
    
            //if (yl-yr>5 || yr-yl>5) // may be wrong correspondence, discard. But vector can't be changed during iteration
            //{}
    
            pt3D.z = -focalLenInPixel * baselineInMM / (xl-xr); // xl should be larger than xr, if xl is shot by the left camera
            pt3D.y = -(-ylr + imgH/2) * pt3D.z / focalLenInPixel;
            pt3D.x = (imgW/2 - xl) * pt3D.z / focalLenInPixel;
    
            minX = min(minX, pt3D.x); maxX = max(maxX, pt3D.x);
            minY = min(minY, pt3D.y); maxY = max(maxY, pt3D.y);
            minZ = min(minZ, pt3D.z); maxZ = max(maxZ, pt3D.z);
            pts3D.push_back(pt3D);
    
            if ((cnt++)%showIntv == 0)
            {
                Scalar color = CV_RGB(rand()&64,rand()&64,rand()&64);
                sprintf_s(str, 100, "%.0f,%.0f,%.0f", pt3D.x, pt3D.y, pt3D.z);
                putText(imgShow, str, Point(xl-13,ylr-3), FONT_HERSHEY_SIMPLEX, .3, color);
                circle(imgShow, *iterL, 2, color, 3);
            }
    
        }
    
        imshow("back project", imgShow);
        waitKey();
    
        center3D.x = (minX+maxX)/2;
        center3D.y = (minY+maxY)/2;
        center3D.z = (minZ+maxZ)/2;
        size3D[0] = maxX-minX;
        size3D[1] = maxY-minY;
        size3D[2] = maxZ-minZ;
    }

    3.三角剖分

    3.1 三角剖分简介

    三角剖分是为了之后的纹理贴图,我用了OpenCV中的Delaunay三角剖分函数,这种剖分算法的可以使所形成的三角形的最小角最大。剖分的示例如下:

    这里写图片描述
    OpenCV使用Delaunay算法将平面分割成小的三角形区域(该三角形确保包括所有的分割点)开始不断迭代完成。在这种情况下,对偶划分就是输入的二维点集的Voronoi图表。这种划分可以用于对一个平面进行三维分段变换、形态变换、平面点的快速 定位以及建立特定的图结构(如NNG,RNG)。

    这里写图片描述
    同时由表可以看出,三角网生成法的时间效率最低,分治算法的时间效率最高,逐点插入法效率居中。

    3.2 Bowyer-Watson算法

    目前采用逐点插入方式生成的Delaunay三角网的算法主要基于Bowyer-Watson算法,Bowyer-Watson算法的主要步骤如下:

    1)建立初始三角网格:针对给定的点集V,找到一个包含该点集的矩形R,我们称R为辅助窗口。连接R的任意一条对角线,形成两个三角形,作为初始Delaunay三角网格。

    2)逐点插入:假设目前已经有一个Delaunay三角网格T,现在在它里面再插入一个点P,需要找到该点P所在的三角形。从P所在的三角形开始,搜索该三角形的邻近三角形,并进行空外接圆检测。找到外接圆包含点P的所有的三角形并删除这些三角形,形成一个包含P的多边形空腔,我们称之为Delaunay空腔。然后连接P与Delaunay腔的每一个顶点,形成新的Delaunay三角网格。

    3)删除辅助窗口R:重复步骤2),当点集V中所有点都已经插入到三角形网格中后,将顶点包含辅助窗口R的三角形全部删除。

    在这些步骤中,快速定位点所在的三角形、确定点的影响并构建Delaunay腔的过程是每插入一个点都会进行的。随着点数的增加,三角形数目增加很快,因此缩短这两个过程的计算时间,是提高算法效率的关键。
    算法执行图示如下:
    这里写图片描述

    3.3 三角剖分代码分析

    三角剖分的代码见cvFuncs.cpp中的TriSubDiv函数,我将特征点存储到一个vector变量中,剖分结果存储到一个vector变量中,Vec3i中存储的是3个表示顶点编号的整数。

    我们需要存储Delaunay的内存空间和一个外接矩形(该矩形盒子用来确定虚拟三角形)

    // STORAGE AND STRUCTURE FOR DELAUNAY SUBDIVISION //存储和结构 for三角剖分  
    //  
    CvRect rect = { 0, 0, 600, 600 };  //Our outer bounding box //我们的外接边界盒子  
    CvMemStorage* storage;    //Storage for the Delaunay subdivsion //用来存储三角剖分  
    storage = cvCreateMemStorage(0);    //Initialize the storage //初始化存储器  
    CvSubdiv2D* subdiv; //The subdivision itself // 细分  
    subdiv = init_delaunay( storage, rect);   //See this function below //函数返回CvSubdiv类型指针  
    init_delaunay函数如下,它是一个OpenCV函数,是一个包含一些OpenCV函数的函数包。
    //INITIALIZATION CONVENIENCE FUNCTION FOR DELAUNAY SUBDIVISION //为三角剖分初始化便利函数  
    //  
    CvSubdiv2D* init_delaunay(CvMemStorage* storage,CvRect rect) {  
    CvSubdiv2D* subdiv;  
    subdiv = cvCreateSubdiv2D(CV_SEQ_KIND_SUBDIV2D,sizeof(*subdiv),sizeof(CvSubdiv2DPoint),sizeof(CvQuadEdge2D),storage);//为数据申请空间  
    cvInitSubdivDelaunay2D( subdiv, rect ); //rect sets the bounds  
    return subdiv;//返回申请空间的指针  
    }  

    我们知道三角剖分是对散点集进行处理的,我们知道了散点集就可以获得点集的三角剖分。如何传入(插入)散点集呢?
    这些点必须是32位浮点型,并通过下面的方式插入点:

    CvPoint2D32f fp; //This is our point holder//这是我们点的持有者(容器)  
    for( i = 0; i < as_many_points_as_you_want; i++ ) {  
    // However you want to set points //如果我们的点集不是32位的,在这里我们将其转为CvPoint2D32f,如下两种方法。  
    //  
    fp = your_32f_point_list[i];  
    cvSubdivDelaunay2DInsert( subdiv, fp );  
    }  

    转换为CvPoint2D32f的两种方法:
    1)通过宏cvPoint2D32f(double x,double y)
    2)通过cxtype.h下的cvPointTo32f(CvPoint point)函数将整形点方便的转换为32位浮点型。
    当可以通过输入点(散点集)得到Delaunay三角剖分后,接下来,我们用一下两个函数设置和清除相关的Voronoi划分:

    cvCalcSubdivVoronoi2D( subdiv ); // Fill out Voronoi data in subdiv //在subdiv中填充Vornoi的数据  
    cvClearSubdivVoronoi2D( subdiv ); // Clear the Voronoi from subdiv//从subdiv中清除Voronoi的数据  

    CvSubdiv2D结构如下:

    #define CV_SUBDIV2D_FIELDS() \  
    CV_GRAPH_FIELDS() \  
    int quad_edges; \  
    int is_geometry_valid; \  
    CvSubdiv2DEdge recent_edge; \  
    CvPoint2D32f topleft; \  
    CvPoint2D32f bottomright;  
    typedef struct CvSubdiv2D  
    {  
    CV_SUBDIV2D_FIELDS()  
    }  
    CvSubdiv2D;  
    #define CV_GRAPH_FIELDS()               \  
    CV_SET_FIELDS() /* set of vertices */   \  
    CvSet *edges;  /* set of edges    */  
    #define CV_SET_FIELDS()                                            \  
    CV_SEQUENCE_FIELDS()             /*inherits from [#CvSeq CvSeq] */ \  
    struct CvSetElem* free_elems;   /*list of free nodes           */  

    整体代码如下:

    void TriSubDiv( vector<Point2f> &pts, Mat &img, vector<Vec3i> &tri ) 
    {
        CvSubdiv2D* subdiv;//The subdivision itself // 细分 
        CvMemStorage* storage = cvCreateMemStorage(0); ;//Storage for the Delaunay subdivsion //用来存储三角剖分 
        Rect rc = Rect(0,0, img.cols, img.rows); //Our outer bounding box //我们的外接边界盒子 
    
        subdiv = cvCreateSubdiv2D( CV_SEQ_KIND_SUBDIV2D, sizeof(*subdiv),
            sizeof(CvSubdiv2DPoint),
            sizeof(CvQuadEdge2D),
            storage );//为数据申请空间  
    
        cvInitSubdivDelaunay2D( subdiv, rc );//rect sets the bounds 
    
        //如果我们的点集不是32位的,在这里我们将其转为CvPoint2D32f,如下两种方法。
        for (size_t i = 0; i < pts.size(); i++)
        {
            CvSubdiv2DPoint *pt = cvSubdivDelaunay2DInsert( subdiv, pts[i] );
            pt->id = i;
        }
    
        CvSeqReader reader;
        int total = subdiv->edges->total;
        int elem_size = subdiv->edges->elem_size;
    
        cvStartReadSeq( (CvSeq*)(subdiv->edges), &reader, 0 );
        Point buf[3];
        const Point *pBuf = buf;
        Vec3i verticesIdx;
        Mat imgShow = img.clone();
    
        srand( (unsigned)time( NULL ) );   
        for( int i = 0; i < total; i++ ) 
        {   
            CvQuadEdge2D* edge = (CvQuadEdge2D*)(reader.ptr);   
    
            if( CV_IS_SET_ELEM( edge )) 
            {
                CvSubdiv2DEdge t = (CvSubdiv2DEdge)edge; 
                int iPointNum = 3;
                Scalar color = CV_RGB(rand()&255,rand()&255,rand()&255);
    
                //bool isNeg = false;
                int j;
                for(j = 0; j < iPointNum; j++ )
                {
                    CvSubdiv2DPoint* pt = cvSubdiv2DEdgeOrg( t );
                    if( !pt ) break;
                    buf[j] = pt->pt;
                    //if (pt->id == -1) isNeg = true;
                    verticesIdx[j] = pt->id;
                    t = cvSubdiv2DGetEdge( t, CV_NEXT_AROUND_LEFT );
                }
                if (j != iPointNum) continue;
                if (isGoodTri(verticesIdx, tri))
                {
                    //tri.push_back(verticesIdx);
                    polylines( imgShow, &pBuf, &iPointNum, 
                        1, true, color,
                        1, CV_AA, 0);
                    //printf("(%d, %d)-(%d, %d)-(%d, %d)\n", buf[0].x, buf[0].y, buf[1].x, buf[1].y, buf[2].x, buf[2].y);
                    //printf("%d\t%d\t%d\n", verticesIdx[0], verticesIdx[1], verticesIdx[2]);
                    //imshow("Delaunay", imgShow);
                    //waitKey();
                }
    
                t = (CvSubdiv2DEdge)edge+2;
    
                for(j = 0; j < iPointNum; j++ )
                {
                    CvSubdiv2DPoint* pt = cvSubdiv2DEdgeOrg( t );
                    if( !pt ) break;
                    buf[j] = pt->pt;
                    verticesIdx[j] = pt->id;
                    t = cvSubdiv2DGetEdge( t, CV_NEXT_AROUND_LEFT );
                }   
                if (j != iPointNum) continue;
                if (isGoodTri(verticesIdx, tri))
                {
                    //tri.push_back(verticesIdx);
                    polylines( imgShow, &pBuf, &iPointNum, 
                        1, true, color,
                        1, CV_AA, 0);
                    //printf("(%d, %d)-(%d, %d)-(%d, %d)\n", buf[0].x, buf[0].y, buf[1].x, buf[1].y, buf[2].x, buf[2].y);
                    //printf("%d\t%d\t%d\n", verticesIdx[0], verticesIdx[1], verticesIdx[2]);
                    //imshow("Delaunay", imgShow);
                    //waitKey();
                }
            }
    
            CV_NEXT_SEQ_ELEM( elem_size, reader );
    
        }
    
        //RemoveDuplicate(tri);
        char title[100];
        sprintf_s(title, 100, "Delaunay: %d Triangles", tri.size());
        imshow(title, imgShow);
        waitKey();
    }

    平面划分是将一个平面分割为一组不重叠的、能够覆盖整个平面的区域。结构CvSubdiv2D描述了建立在二维点集上的划分结构,其中点集互相连接且构成平面图形,该图形通过结合一些无线连接外部划分点(称为凸形点)的边缘,将一个平面用按照其边缘划分成很多小区域。

    对于每一个划分操作,都有一个对偶划分与之对应。对偶的意思是小区域与点(划分的顶点)变换角色,即在对偶划分中,小区域被当做一个顶点(以下称为虚拟点)而原始的划分顶点被当做小区域。如下图所示,原始的划分用实线表示,而对偶划分用虚线表示。

    4.三维重构

    为了保证三维重建的效果,一般地要对深度图像进行后续处理。要从深度图像中恢复高质量的视差图,对深度图像的要求有:
    ①深度图像中,物体的边界必需与图像中物体的边界对齐;
    ②在场景图中,深度图像要尽可能均勻和平滑,即对图像进行平滑处理。

    三维重构的思路很简单,用OpenGL中纹理贴图功能,将平面图像中的三角形逐个贴到计算出的三维坐标上去就可以了。为了便于观察3D效果,我还设计了交互功能:用方向键可以上下左右旋转重构的模型,用鼠标滚轮可以放大或缩小。用gluLookAt函数可以实现视点旋转的功能。三维重构的代码实现在glFuncs.cpp中。

    纹理贴图:

    GLuint Create3DTexture( Mat &img, vector<Vec3i> &tri, 
                           vector<Point2f> pts2DTex, vector<Point3f> &pts3D, 
                            Point3f center3D, Vec3f size3D ) 
    {
        GLuint tex = glGenLists(1);
        int error = glGetError();
        if (error != GL_NO_ERROR) 
            cout << "An OpenGL error has occured: " << gluErrorString(error) << endl;
        if (tex == 0) return 0;
    
        Mat texImg;
        cvtColor(img, img, CV_BGR2RGB);
        resize(img, texImg, Size(512,512)); // seems no need to do this
    
        glNewList(tex, GL_COMPILE);
    
        vector<Vec3i>::iterator iterTri = tri.begin();
        //vector<Point3f>::iterator iterPts3D = pts3D.begin();
        Point2f pt2D[3];
        Point3f pt3D[3];
    
        glDisable(GL_BLEND);
        glEnable(GL_TEXTURE_2D);
        for ( ; iterTri != tri.end(); iterTri++)
        {
            Vec3i &vertices = *iterTri;
            int ptIdx;
            for (int i = 0; i < 3; i++)
            {
                ptIdx = vertices[i];
                if (ptIdx == -1) break;
                //else cout<<ptIdx<<"\t";
                pt2D[i].x = pts2DTex[ptIdx].x / img.cols;
                pt2D[i].y = pts2DTex[ptIdx].y / img.rows;
                pt3D[i] = (pts3D[ptIdx] - center3D) * (1.f / max(size3D[0],size3D[1]));
                //pt3D[i].z -= offset;
            }
    
            if (ptIdx != -1)
            {
                MapTexTri(texImg, pt2D, pt3D);
                //cout<<endl;
            }
        }
        glDisable(GL_TEXTURE_2D);
    
        glEndList();
        return tex;
    
    }

    效果展示及不足
    Cloth图像是重构效果比较好的一组:

    这里写图片描述

    可以比较明显的看出3D效果,也比较符合直觉。然而其他图像效果就差强人意了:

    这里写图片描述

    仔细分析造成这种效果的原因,一方面,特征点的匹配可能有些误差,造成3D坐标的计算不太精确,但大部分坐标还是准确的。另一方面,左右视图可能会有不同的遮挡、偏移等情况,因此匹配得到的特征点可能实际上并不是3维世界中的同一点,这种误差是无法消除的。但造成效果变差的最重要的原因,还是图像中深度变化较大,而特征点选取的比较稀疏,因此正面看还比较正常,一旦旋转纹理就显得扭曲变形了。为了解决这个问题,应当试图把特征点取到深度变化较剧烈的地方,一般是图像中的边界处。然而特征点检测一般都检测出的是角点和纹理密集的区域,因此可以考虑更换对应点匹配的方法。

    如果要进一步改进效果,可以先对视差图像进行分割,将图像分成视差比较连续的几块区域分别贴图,视差变化剧烈的区域就不必把扭曲的纹理贴上去了。我尝试了以下分割的效果,如下图所示,应该可以达到更好的效果,不过由于时间所限,就没有进一步实现下去了。

    关于上面实现的两种求取视差的算法,在main函数的前面设置了一个变量g_algo,可以用来切换不同的算法。

    参考文献:

    立体匹配原理:

    http://blog.csdn.net/wangyaninglm/article/details/51533549
    http://blog.csdn.net/wangyaninglm/article/details/51531333
    三维重建原理:

    http://blog.csdn.net/wangyaninglm/article/details/51558656
    http://blog.csdn.net/wangyaninglm/article/details/51558310
    三角剖分原理:

    http://blog.csdn.net/newthinker_wei/article/details/45598769
    http://www.learnopencv.com/delaunay-triangulation-and-voronoi-diagram-using-opencv-c-python/
    这篇文章其实主要是针对早期下到的一个代码和文档的总结,和一些个人资料的总结,由于时间比较早,找不到出处了,如果原作者看到了觉的不妥,那我就把它改成转载啦,嘿嘿嘿。

    代码下载

    CSDN: http://download.csdn.net/detail/wangyaninglm/9597622
    github:https://github.com/wynshiter/OpenCV-OpenGL–Reconstuction3d

    展开全文
  • 双目视觉三维重建

    2015-04-08 16:44:42
    这是双目视觉三维重建代码,希望对大家有所帮助
  • 基于双目视觉的深度计算和三维重建双目视觉opencv opengl三维重建,简单的三维重建系统,代码可正常运行
  • 双目视觉opencv opengl三维重建双目视觉opencv opengl三维重建双目视觉opencv opengl三维重建,请修改代码中opencv对应版本号
  • Evision双目视觉关于双目视觉的一些总结相机模型标定视差算法:立体匹配重投影:测量,三维重建,重投影约束三维重建示例程序 关于双目视觉的一些总结 笔者2013年进入吉林大学软件学院,2014年开始写自己的第一个完整的...

    关于双目视觉的一些总结

    说明

    如果读者对于本文或者Evision程序有任何问题,请尽量通过github的issue向我反馈,如果我回复不及时,请发邮件给我,博客的回复往往不及时,请谅解.

    前言

    笔者2013年进入吉林大学软件学院,2014年开始写自己的第一个完整的程序,期间受到过无数前辈的帮助,正是这个程序的完成给了我极大的信心,也让我喜欢上编程.这个程序是"基于OpenCV的双目测距",他的主要代码来自于邹宇华老师的OpenCV例程,我只不过进行了一些小小的修改,然后做了一个界面,在收获了许多经验的同时,也发现了双目视觉,乃至机器视觉的中文社区环境的一个问题,那就是有效资源非常少,很多人的博客只不过是互相转载,很多人提供的代码要么语焉不详,要么就是其他示例程序甚至是官方实例的修修补补,这对于初学者非常的不友好,一知半解的教程会让新手走很多弯路.在这期间,我的双目视觉程序仍然维持着更新,我决定把相关的资料全部公开,如果您发现错误,请及时留言评论,同时欢迎转载.文末附有程序代码链接,如果觉得有帮助,请加star,谢谢.

    双目视觉程序: https://github.com/jiafeng5513/Evision
    新版程序演示视频:https://www.bilibili.com/video/av46024738
    旧版程序演示视频:https://www.bilibili.com/video/av8862669

    相机模型

    通常关于相机模型的文章只会提到小孔成像模型,这里要提醒的是,小孔成像模型并不能代表所有种类的相机,为了彻底搞清楚计算机视觉的基础,我们必须把相机模型弄明白.
    1. 相机的投影模型
    图1.相机投影

    • 球面投影模型可以概括所有的相机模型,或者说所有的相机模型都是球面投影模型的一种特殊情况.
    • 我们可以在光心和成像平面之间假想一个球面,物体先在球面上成像,再从球面上投影到成像平面上.
    • 透镜的光学性质都是中心对称的,所以球面假设很合理.
    • 穿过光心的光线是直线传播的,不同模型的角度差异的来源在于球面向平面投影的不同方法.
    • 从球面向平面的投影,最常见的用处是绘制世界地图(可以百度:墨卡托投影,圆锥投影,等积投影等),
      从球面向平面的投影有很多种方法,相机在制造时,镜头由很多片透镜组成,这些透镜的左右就是逼近某一种投影方法,使得透过镜头的光线直接透射在传感器平面上的情况等效于先投影在成像球面,再向成像平面投影,物体与光轴的夹角θ变化时,像点与图像中心的距离R也随之变化,按照变化关系的不同,经典相机模型被分为5类,其中透视投影模型又称为小孔模型.
    计算公式 英文名 中文名
    R=f*tan(θ) Perspective 透射投影
    R=f*θ Equidistant 等距投影
    R=2fsin(θ/2) Equisolid angle 等立体角投影(等积投影)
    R=f*sin(θ) Orthographic 正交投影
    R=2ftan(θ/2) Stereographic 球极投影(体视投影,保角投影)

    图2 五种投影模型
    关于鱼眼镜头的投影特点,推荐大家看这篇知乎专栏,里面有一个视频直观的描述了这些投影的成像特点.

    2. 不同类型的镜头

    • 由于镜头制造上的限制,无法制造出满足理想情况的镜头,所以设计镜头时采用一个多项式逼近理想情况,其中K称为径向畸变系数.
    • 由于实际镜头在拟合理想投影模型略有不同,这种性质是相机镜头的一种属性,用K描述.
    • 径向畸变是最主要的畸变来源.此外,相机还有其他畸变,例如切向畸变,薄棱镜畸变,传感器倾斜畸变.
      在这里插入图片描述
      通过观察R和θ的曲线,我们能获知这些投影方式的一些性质:
      图4 R和θ的关系
      根据图像我们能得到一些结论:
    • 标准镜头(使用透射投影)的视角无法达到180°
    • 鱼眼镜头能够实现180°的视角.
    • 虽然理论上四种投影类型的鱼眼镜头都能实现180°视角,但是其最大视角对应的r是不同的,在传感器尺寸固定的情况下,不同的鱼眼相机实际能实现的最大视角不相同.
    • 标准镜头在被摄物体到光心与光轴之间的夹角接近90°的过程中,单位角度变化所对应的成像点位移迅速增大,可以预见,图片边缘的情况一定和图片中心有所区别.
      关于标准模型,这里提醒大家:
      在一些涉及到相机模型的场景,要关注其中使用的相机模型是否符合预期,小孔模型使用的比较多,但是在使用广角相机和鱼眼相机时,小孔模型未必是有效的,此外在有些场景下,小孔模型过多的畸变系数可能会引起麻烦(例如训练深度神经网络时,过多的参数可能会引入复杂的可训练性的证明).
      3. 小孔成像模型
      在这里插入图片描述
      在这里插入图片描述
    • (X,Y,Z) 是被摄物体的世界坐标
    • (u,v) 是投影点在像素坐标系下的坐标
    • A是相机矩阵,也叫内参矩阵,描述从相机坐标系到图像坐标系的变换
    • (cx,cy)是主点坐标,通常位于图像中心
    • fx,fy是以像素单位表示的焦距
    • S是缩放系数
    • [R|t] 也叫外参矩阵,描述由世界坐标系变换到相机坐标系

    外参数把世界坐标变换到相机坐标,内参数把相机坐标变换到图像坐标
    在这里插入图片描述
    考虑一些畸变:
    在这里插入图片描述
    其中k是径向畸变系数,P是切向畸变系数,S是薄棱镜畸变系数,在大多数情况下,只考虑径向和切向畸变的低阶系数即可,在实验精度有限的前提下,使用高阶系数可能会引入很大的误差.此外,薄棱镜畸变和倾斜传感器畸变在大多数相机上都非常小
    在这里插入图片描述
    径向畸变K1的正负会影响畸变的视觉效果,负数时一般称为桶形畸变,正数时称为枕形畸变

    4. 极线约束,本质矩阵和基本矩阵
    在这里插入图片描述

    • 相机从两个位姿O1和O2拍摄同一个点P,分别成像p1和p2,这里p1和p2是相机坐标系下的坐标.PO1O2共面,两个位姿的成像面为
    • I1I_1I2I_2,分别交平面PO1O2PO_1O_2于交线l1:p1e1l_1:p_1e_1l2:p2e2l_2:p_2e_2
    • 把直线参数表示成: l=[abc]Tl=\begin{bmatrix}a & b & c\end{bmatrix}^{T}
    • 存在矩阵E,使得Ep1=l2Ep_1=l_2
    • p1Tl1=0,p2Tl2=0p_{1}^{T}l_1=0,p_{2}^{T}l_2=0
    • p2TEp1=0p_{2}^{T}Ep_1=0
    • p1p_1p2p_2在图像坐标系下的坐标分别为p1p'_1p2p'_2,K为相机矩阵
    • Kp1=p1,Kp2=p2Kp_1= p'_1, Kp_2= p'_2.也就是p1=K1p1,p2=K1p2p_1=K^{-1}p'_1, p_2=K^{-1}p'_2
    • p2TEp1=0=(K1p2)TEK1p1=(p2)T(K1)TEK1p1p_2^{T} Ep_1=0=(K^{-1}p'_2)^{T}EK^{-1}p'_1=(p'_2)^{T}(K^{-1})^{T}EK^{-1}p'_1
    • (K1)TEK1=F(K^{-1})^{T}EK^{-1}=F,则有:(p2)TF(p1)=0(p'_2)^{T}F(p'_1)=0
    • 其中E为Essential Matrix(本质矩阵)
    • F为Fundamental Matrix(基础矩阵)
    • 本质矩阵和基础矩阵对于一对图片是固定值
    • 使用本质矩阵或者基础矩阵,可以从一个视角的一个像点确定另一个视角的一条线,所以这个性质被称为极限约束

    在这里插入图片描述

    • 从另一角度理解:Pac=O1PP_{a}^{c}= \vec{O_1P}
    • 相机由O1O_1运行到O2O_2,相当于坐标系变换:T=O1O2Pbc=O2P=RPac+tT=\stackrel{\longrightarrow}{O_1O_2}\cdot P_{b}^{c}= \stackrel{\longrightarrow}{O_2P}=RP_{a}^{c}+t
    • 由于空间中O1O2PO_1O_2P三点共面,所以O2P(O1O2×O1P)=0\stackrel{\longrightarrow}{O_2P} \cdot(\stackrel{\longrightarrow}{O_1O_2}\times \stackrel{\longrightarrow}{O_1P})=0
    • O2O_2的坐标系下可以表示为:Pbc(T×(RPac))=0P_{b}^{c} \cdot(T\times(RP_{a}^{c}))=0

    5. 单应矩阵(Homography matrix)

    • Homography matrix表达两个平面之间的投影关系,例如从三维空间中的实际平面到像平面,或者一个像平面到另一个像平面.
    • 两张图片上的对应点的齐次坐标分别是(x’,y’,1)和(x,y,1),则有以下等式,其中H为单应矩阵:

    • 齐次化之后H中有8个自由量,至少需要8个方程也就是四对点来求解
    • 实际平面到像平面的单应矩阵中包含了相机的内外参数和畸变

    标定

    相机的标定,目的是获取相机的参数和畸变系数,大致分为以下三类:

    标定方法 优点 缺点 常用方法
    传统相机标定法 可使用于任意的相机模型、 精度高 需要标定物 Tsai两步法、张氏标定法
    主动视觉相机标定法 不需要标定物、算法简单、鲁棒性高 成本高、设备昂贵 主动系统控制相机做特定运动
    相机自标定法 灵活性强、可在线标定 精度低、鲁棒性差 分层逐步标定、基于Kruppa方程

    最常用的标定方法是张氏标定,下面简要介绍张氏标定的原理:
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

    视差算法:立体匹配

    立体匹配算法也称为视差算法,目的是求解两幅图片之间的视差值.首先我们来定义视差值这个概念:

    定义1:三维世界中有某一定点P,相机在两个位置对P成像,分别获取像点p1p_1p2p_2,则称p1p_1p2p_2同名点
    定义2:如果两张图片经过了畸变矫正,并且同名点的Y坐标在各自的像素坐标系下相等,我们称这两张图片为已校正图片对
    定义3:在双目立体视觉中,两个相机的视角方向朝大致相同的方向,此时观察者站在相机的后方,视线方向与两台相机一致,他左手边的相机我们称之为 左相机,拍摄的图片称为 左视图
    定义4:已校正图片对上,同名点的是X坐标之差,为视差值,逐像素计算式差值后组成的图片,称为 视差图 .特别地,在双目立体视觉中,以左视图为基准,通过同名点在右图上的坐标减去其在左图上的坐标得到的,称为左视差(图),相对应的称为右视差图

    通过定义我们发现,立体匹配问题实际上是求同名点问题,因为只要找到准确的同名点,视差值不过是一个坐标差.可以说绝大多数视差算法的核心都是寻找同名点(一些基于深度学习的视差方法并不显式的求同名点).更具体点来说,当要求左视差图的时候,就是从左视图中依次取出每个像素,然后到右视图的某个范围内找 这个像素的同名点,这样当我们找完左视图上所有像素的同名点时,就得到了一张和左视图一样大的视差图了.

    上面是通俗化的解释,接下来我们进行形式化的说明.视差算法是一个机器视觉的传统领域,经过多年的发展已经相当成熟,根据Schrstein和Szeliski的总结,双目立体匹配可划分为四个步骤:匹配代价计算、代价聚合、视差计算和视差优化。

    1. 匹配代价计算
      匹配代价计算的目的是衡量待匹配像素与候选像素之间的相关性。两个像素无论是否为同名点,都可以通过匹配代价函数计算匹配代价,代价越小则说明相关性越大,是同名点的概率也越大。

      每个像素在搜索同名点之前,往往会指定一个视差搜索范围D(Dmin ~ Dmax),视差搜索时将范围限定在D内,用一个大小为W×H×D(W为影像宽度,H为影像高度)的三维矩阵C来存储每个像素在视差范围内每个视差下的匹配代价值。矩阵C通常称为DSI(Disparity Space Image)。

      匹配代价计算的方法有很多,传统的摄影测量中,使用灰度绝对值差(AD,Absolute Differences)[1]、灰度绝对值差之和(SAD,Sum of Absolute Differences)、归一化相关系数(NCC,Normalized Cross-correlation)等方法来计算两个像素的匹配代价;计算机视觉中,多使用互信息(MI,Mutual Information)法[2,3]、Census变换(CT,Census Transform)法[4,5]、Rank变换(RT, Rank Transform)法[6,7]、BT(Birchfield and Tomasi)法[8]等作为匹配代价的计算方法。不同的代价计算算法都有各自的特点,对各类数据的表现也不尽相同,选择合适的匹配代价计算函数是立体匹配中不可忽视的关键步骤。

      无论是什么样的代价计算方法,都可以近似的看作这样的一个过程:首先在左图中确定一个像素或者一块像素区域,然后在右图中用某种形状的模板或者策略选定一块匹配区域,这个区域的大小和滑动步长都是超参数控制的,然后计算代价,移动模板窗口,再计算代价值,直到耗尽搜索范围,然后左图中选定第二块或者第二个像素,再重复这个过程。

    2. 代价聚合
      在这里插入图片描述
      代价聚合的根本目的是让代价值能够准确的反映像素之间的相关性。上一步匹配代价的计算往往只会考虑局部信息,通过两个像素邻域内一定大小的窗口内的像素信息来计算代价值,这很容易受到影像噪声的影响,而且当影像处于弱纹理或重复纹理区域,这个代价值极有可能无法准确的反映像素之间的相关性,直接表现就是真实同名点的代价值非最小。

      而代价聚合则是建立邻接像素之间的联系,以一定的准则,如相邻像素应该具有连续的视差值,来对代价矩阵进行优化,这种优化往往是全局的,每个像素在某个视差下的新代价值都会根据其相邻像素在同一视差值或者附近视差值下的代价值来重新计算,得到新的DSI,用矩阵S来表示。

      实际上代价聚合类似于一种视差传播步骤,信噪比高的区域匹配效果好,初始代价能够很好的反映相关性,可以更准确的得到最优视差值,通过代价聚合传播至信噪比低、匹配效果不好的区域,最终使所有影像的代价值都能够准确反映真实相关性。常用的代价聚合方法有扫描线法、动态规划法、SGM算法中的路径聚合法等。

    3. 视差计算
      视差计算即通过代价聚合之后的代价矩阵S来确定每个像素的最优视差值,例如使用赢家通吃算法(WTA,Winner-Takes-All)来计算,如图2所示,即某个像素的所有视差下的代价值中,选择最小代价值所对应的视差作为最优视差。这一步非常简单,这意味着聚合代价矩阵S的值必须能够准确的反映像素之间的相关性,也表明上一步的代价聚合步骤是立体匹配中极为关键的步骤,直接决定了算法的准确性。
      在这里插入图片描述

    4. 视差优化
      视差优化的目的是对上一步得到的视差图进行进一步优化,是一些后处理方法的组合。改善视差图的质量,包括剔除错误视差、适当平滑以及子像素精度优化等步骤,一般采用左右一致性检查(Left-Right Check)算法剔除因为遮挡和噪声而导致的错误视差;采用剔除小连通区域算法来剔除孤立异常点;采用中值滤波(Median Filter)、双边滤波(Bilateral Filter)等平滑算法对视差图进行平滑;另外还有一些有效提高视差图质量的方法如鲁棒平面拟合(Robust Plane Fitting)、亮度一致性约束(Intensity Consistent)、局部一致性约束(Locally Consistent)等也常被使用。

      由于WTA算法所得到的视差值是整像素精度,为了获得更高的子像素精度,需要对视差值进行进一步的子像素细化,常用的子像素细化方法是一元二次曲线拟合法,通过最优视差下的代价值以及左右两个视差下的代价值拟合一条一元二次曲线,取二次曲线的极小值点所代表的视差值为子像素视差值。如图3所示。

      在编程实践中这里实际上最容易出现问题。OpenCV实现的BM和SGBM并没有亚像素增强的步骤,所以它输出的视差数据是CV8U也就是unsingned char类型的,他的视差等级非常有限,比如说设置的最大视差是50,那他的视差就是0-49这50个整数,精度就很差了;如果进行亚像素增强,输出的数据格式就是float或者double的,这是不能直接在屏幕上显示的,必须要进行归一化,所以如果视差数据只是整个程序的中间结果,那么就不可以保存成图片再读取,这样就会损失精度。

    5. 代价函数
      代价函数的想法是,确定了同名点实际上视差就确定了,那么我们输入左点和视差,在给定的视差范围内计算代价值,代价最低的那个情况就最可能是同名点,此时的视差也最可信.下面列举几种常见的代价函数:
      在这里插入图片描述
      在这里插入图片描述
      在这里插入图片描述
      在这里插入图片描述
      在这里插入图片描述

    6. OpenCV-BM
      OpenCV-BM算法使用SAD作为代价函数,在源码中我们看到该方法被称为“基于SAD的快速立体匹配算法”.根据对源码的解析,我把OpenCV-BM的流程分为四个部分:校验、预处理、视差计算和后处理.

      校验:校验算法的输入和超参数的选取是否合理.左右视图大小相等,格式均为CV_8UC1,接收视差图的Mat为CV_16SC1或者CV_32FC1, preFilterSize是[5,255]之间的奇数, preFilterCap在[1,63], SADWindowSize是[5,255]之间的奇数且不能大于图片的宽或高, numDisparities必须能被16整除, textureThreshold必须是非负数, uniquenessRatio必须是非负数.

    超参数 说明
    1 preFilterType 前处理可以是NORMALIZED或XSOBEL
    2 preFilterSize 归一化窗口大小,对于XSobel没有意义
    3 preFilterCap 截断值.对于XSobel(X方向上的Sobel),如果滤波后结果在[-prefilterCap,preFilterCap]之间,对应取[0, 2*preFilterCap]
    4 SADWindowSize SAD窗口大小
    5 minDisparity 视差搜索起点
    6 numDisparities 视差窗口,一个推荐值是((width / 8) + 15) & (~0xfl)
    7 textureThreshold 低纹理区域的判断阈值. 如果当前 SAD 窗口内所有邻居像素点的 x-导数绝对值之和小于指定阈值,则意味着当前窗口内的图像变化率太低也就是纹理很差,则该窗口对应的像素点的视差值设为0.
    8 uniquenessRatio 视差唯一性百分比. 视差窗口范围内最低代价是次低代价的 (1 + uniquenessRatio / 100) 倍时,最低代价对应的视差值才是该像素点的视差, 否则该像素点的视差为 0
    9 speckleRange 散斑窗口内允许的最大波动值
    10 speckleWindowSize 散斑滤波窗口大小,<=0时不进行散斑滤波
    11 disp12MaxDiff 左右视差图的最大容许差异,超标清零,左右一致性检查的目的是找到遮挡点,OpenCV默认的方法直接去除掉了遮挡点的视差,小于0时跳过
    12 dispType 视差图数据类型

    上表是OpenCV-BM的超参数,可以看到OpenCV实现的BM算法有很多地方还是很暴力的,它的前处理和后处理都很简单,再后来的很多论文中着墨甚多的后处理步骤例如遮挡点和miss点的处理在OpenCV中都极其简陋.然而OpenCV的源码水平非常高,对于内存的使用和优化,变量类型的使用等问题绝不含糊,非常值得学习.有人对源码进行了注释[26].
    预处理:X方向的sobel边缘检测和归一化
    视差计算:SAD带价函数
    后处理:左右视差一致性检查,散斑滤波

    1. OpenCV-SGBM
      也称为SGM、半全局匹配算法等[9]。
      观察OpenCV-BM的实现可知,在计算某一点的视差(代价)时,只考虑一个窗口内的像素,窗口外面的图像还有很多像素,但是并不会参与到本点的视差计算当中.事实上有一些视差算法会考虑所有像素的影响,这种算法就叫全局算法,BM是局部算法,而SGM是半全局的.需要注意的是,OpenCV-SGBM(以下简称SGBM)对于原有的SGM(以下简称SGM)做出了一些调整.

      首先看SGM算法.SGM使用分层互信息代价函数,基于动态规划的代价聚合,与OpenCV-BM类似的前处理和后处理.互信息的计算如下:
      在这里插入图片描述
      其中, MII1,I2MI_{I_1,I_2}是互信息(Mutual Information), HIH_I是图II的熵, PiP_i代表某个点ii的概率分布,也就是灰度直方图为ii的点出现的概率;对应地,PI1,I2P_{I_1,I_2}就是两个图对应点i1i_1i2i_2的联合概率分布,也就是:
      在这里插入图片描述
      Kim等人[27]将HI1,I2H_{I_1,I_2}的计算方法做了一个改进,利用泰勒展开转化为求和问题:
      在这里插入图片描述
      其中⨂表示卷积运算,g(i,k)g(i,k)为高斯卷积核。相应地,边缘熵以及边缘概率的计算如下:
      在这里插入图片描述
      这样,互信息的定义更改为:
      在这里插入图片描述
      代价函数CMIC_{MI}定义为:
      在这里插入图片描述
      其中qq是点pp在视差为dd的情况下的对应校正点。原作者使用分层互信息HMI(HMI)进行计算,每一层尺寸减少一半。单次计算的时间复杂度是O(WHD)O(WHD),即width×height×DisparityRangewidth×height×DisparityRange,所以上次迭代将会是当前迭代速度的18\frac{1}{8}
      在这里插入图片描述
      这里1163\frac{1}{16^{3}} 要乘3的原因是小尺寸的随机视差图不靠谱,需要迭代3次。我们可以看到,相比于后文的BT方法仅仅慢了14%.

      以上是使用互信息熵代价函数的SGM,下面介绍SGBM中使用的BT方法,该方法由Birchfield和Tomasi[28]提供.

      对于一个匹配序列MM,其代价函数γ(M)γ(M)表示匹配结果不准确的程度,其值越小越好。
      在这里插入图片描述
      其中,kocck_{occ}表示未匹配的惩罚项(constant occlusion penalty),krk_r表示匹配的奖励项,NoccN_{occ}NrN_r分别表示未匹配和匹配的点数。

      代价聚合:使用动态规划计算[30].设置能量函数E(D):
      在这里插入图片描述
      其中P1和P2分别表示视差差值为1和视差差值大于1的惩罚系数,一般P1<P2。添加两个正则化项一是为了保持视差图平滑,二是为了保持边缘。我们要做的是找到D使得能量函数E(D)最小,但是不幸的是,在二维图像的这个问题是一个NP-完全问题。为了解决这个问题,原文选择沿着一圈8个或者16个方向进行优化。在OpenCV的实现中,默认情况下使用单通道图片,只考虑5个方向,但是当mode = StereoSGBM :: MODE_HH时,将会使用完整的算法,同时消耗更多内存.

      代价聚合公式:
      在这里插入图片描述
      图片的代价聚合是各个像素点代价聚合之和:
      在这里插入图片描述
      选取使代价聚合最小的视差值mindS[emb(q,d),d]min_dS[e_{mb}(q,d),d]即可.
      OpenCV-SGBM的参数如下:
    参数 参数含义
    1 minDisparity 最小视差值
    2 numDisparities 视差范围
    3 blockSize 匹配块大小,奇数,一般在[3,11]
    4 P1 P1是相邻像素视差变化±1的惩罚项
    5 P2 P2是视差变化超过1的惩罚项,P2要大于P1
    6 disp12MaxDiff 左右一致性检查的最大容差,小于等于0禁用检查
    7 preFilterCap 首先计算每个像素的x导数,并将其截断到 [-preFilterCap,preFilterCap],再将结果值传递给BT代价函数
    8 uniquenessRatio 最小代价为次小代价的百分比,常取值为[5,15]
    9 speckleWindowSize 散斑滤波窗口大小,0禁用,取值[50,200]
    10 speckleRange 连通区域的最大容差,内部乘16起效,整数,常[1,2]
    11 mode MODE_SGBM_3WAY:据我观察是微软贡献的一种HAL硬件抽象层加速的SGBM计算方法[29],其他的几种模式都调用默认的计算函数computeDisparitySGBM,包括MODE_SGBM 和MODE_HH
    1. 其他算法:OpenCV-GC,LibElas,ADCensus
      GC首发于2012年的ICPR[10].OpenCV-GC算法在3.0以后已经彻底移除,此外,在OpenCV2.4.13中曾经有过一个OpenCV-VAR算法,后来也移除了,在3.4.11中又出现了binnary_bm和binnary_sgbm算法,由于3.4.11是3的最后一个版本,在4之后,这两个算法也移除了,在此不过多介绍.

      然而不得不说的是,曾经出现的GC算法效果还是很不错的,能够生成非常接近Groundtruth的边缘锐利且平直的视差图,Var算法生成的的视差图很稀疏,不过其中的物体轮廓非常明显,效果类似于canny算子,这两个方法速度很慢且无法进行并行化.

      LibElas效果很不错,由于提供了编写良好的代码[12],应用极其广泛.首发于2011年ACCV[13].LibElas使用了自适应视差范围和基于稳健匹配点的先验知识,且易于并行化,作者给出的实现的确很快,然而效果却不尽如人意.我的观察是作者为了比赛效果使用了很多的后处理,似乎有一些填充过程造成了相反的效果.

      首发于2011年ICCV[11],第三方代码[21]实现效率奇差.创新点是使用了组合的代价函数和cross-based代价聚合,代价函数本身并没有什么亮点,但是一个简单的归一化和权重处理却能优势互补;cross-based代价聚合非常利于GPU实现,整个方法基本没有什么复杂的计算,就效果来看ADCensus的效果非常好,而且参数的适应性很强.

    2. 基于深度学习的视差算法
      一些代表性成果:[14], AdaptWeight[15],次年稍加改动[16],深度学习[17],光流法[18], DERS[19],DERS in pndn[22].此外,在KITTI等比赛上,视差算法的竞争依然在进行着,目前比较多的工作集中在单目视差或深度预测上.

    根据

    测量,三维重建

    此处我们仅在双目视觉的框架下讨论测量和三维重建的问题.

    测量指的是通过双目系统获知图像上一点到相机的距离或者图像上两点之间的距离,三维重建指的获取三维点云.实际上,当我们有了精确的视差值之后,这两个问题就变的容易了.

    在这里插入图片描述
    经过校正后的双目相机模型如上图,其中PleftP_{left}的坐标是(xl,yl)(x_l,y_l),PrightP_{right}的坐标是(xr,yr)(x_r,y_r),视差定义为d=xlxrd=x_l-x_r,ff是焦距,TT是基线长度,则有:
    TZ=TdZf\frac{T} {Z}=\frac{T-d}{Z-f},即Z=fTdZ=\frac{fT}{d}

    示例程序

    双目视觉程序: https://github.com/jiafeng5513/Evision
    新版程序演示视频:https://www.bilibili.com/video/av46024738
    旧版程序演示视频:https://www.bilibili.com/video/av8862669

    参考文献

    1. KANADE T, KANO H, KIMURA S, et al. Development of a video-rate stereo machine: Ieee/rsj International Conference on Intelligent Robots and Systems 95. ‘human Robot Interaction and Cooperative Robots’, Proceedings, 2002[C].
    2. KIM J, KOLMOGOROV V, ZABIH R. Visual Correspondence Using Energy Minimization and Mutual Information: IEEE International Conference on Computer Vision, 2003. Proceedings, 2008[C].
    3. EGNAL G. Mutual Information as a Stereo Correspondence Measure[J]. Technical Reports, 2000.
    4. MA L, LI J, MA J, et al. A Modified Census Transform Based on the Neighborhood Information for Stereo Matching Algorithm: Seventh International Conference on Image and Graphics, 2013[C].
    5. BAIK Y K, JO J H, LEE K M. Fast Census Transform-based Stereo Algorithm using SSE2: The 12th Korea-Japan Joint Workshop on Frontiers of Computer Vision, Tokushima, Japan, 2006[C].
    6. GU Z, SU X, LIU Y, et al. Local stereo matching with adaptive support-weight, rank transform and disparity calibration[J]. Pattern Recognition Letters, 2008,29(9):1230-1235.
    7. BANKS J, BENNAMOUN M, KUBIK K, et al. A constraint to improve the reliability of stereo matching using the rank transform: Acoustics, Speech, and Signal Processing, 1999. on 1999 IEEE International Conference, 1999[C].
    8. BIRCHFIELD S, TOMASI C. A Pixel Dissimilarity Measure That Is Insensitive to Image Sampling[J]. IEEE Transactions on Pattern Analysis & Machine Intelligence, 1998,20(4):401-406.
    9. Stereo Processing by Semiglobal Matching and Mutual Information
    10. Realistic CG Stereo Image Dataset with Ground Truth Disparity Maps
    11. On Building an Accurate Stereo Matching System on Graphics Hardware
    12. http://www.cvlibs.net/software/libelas/
    13. Efficient Large-Scale Stereo Matching
    14. Cross-Scale Cost Aggregation for Stereo Matching, https://github.com/rookiepig/CrossScaleStereo#GF
    15. Locally Adaptive Support-Weight Approach for Visual Correspondence Search
    16. Adaptive Support-Weight Approach for Correspondence Search
    17. Improved Stereo Matching With Constant Highway Networksand Reflective Confidence Learning
    18. Secrets of optical flow estimation and their principles
    19. http://www.fujii.nuee.nagoya-u.ac.jp/multiview-data/mpeg/DE.htm
    20. Stereo Matching by Training a Convolutional Neural Network to Compare Image Patches
    21. https://github.com/DLuensch/StereoVision-ADCensus
    22. http://www.pudn.com/downloads532/sourcecode/graph/texture_mapping/detail2201211.html
    23. http://vision.middlebury.edu/flow/data/
    24. http://www.cvlibs.net/datasets/kitti/eval_object.php
    25. https://yuiwong.org/2017/12/29/ocvbm/
    26. https://zhuanlan.zhihu.com/p/50801189
    27. Visual Correspondence Using Energy Minimization and Mutual Information
    28. Stan Birchfield and Carlo Tomasi. A pixel dissimilarity measure that is insensitive to image sampling. Pattern Analysis and Machine Intelligence, IEEE Transactions on, 20(4):401–406, 1998.
    29. https://github.com/Microsoft/opencv/commit/aea4157340b3136e67c02a5565b152dbe5c73f2e
    30. https://blog.csdn.net/vampireshj/article/details/53813773
    31. 邹宇华StereoVision:https://blog.csdn.net/chenyusiyuan/article/details/8131496
    展开全文
  • 很好的源码,思路清晰,利用双目实现的,大家有兴趣使用的可以直接下载,可以应用到比赛和论文内的,工具是matlab,代码也很好理解
  • 写在前面的话: 一个机器视觉的课程作业,是自行采集一组双目图像,完成立体视觉相关流程:包括相机标定(内参和外参)、畸变校正、基本矩阵估算、视差图计算(需要先进行图像矫正)、恢复并画出3D点坐标。...

    python双目标定及重建

    写在前面的话: 一个机器视觉的课程作业,是自行采集一组双目图像,完成立体视觉相关流程:包括相机标定(内参和外参)、畸变校正、基本矩阵估算、视差图计算(需要先进行图像矫正)、恢复并画出3D点坐标。网上的代码基本上都是基于棋盘格的,初始不懂,当你抄多了,自然就懂了。原理不多做详细介绍,简单介绍一下实现过程以及遇到的问题,还有最后一部分的重建不能确保准确,因为重建出来的图像单纯的不好看!文中的代码引用的都已在参考链接里标注,部分进行了更改和补充。

    仅供学习使用!

    1. 采集图像及预实验

    1.1 图像

    采集一组双目图像,该步骤由师兄本色出镜,本人用两个手机固定位置同时拍摄,以此得到10对图像。也就是左右手机对在一起,然后同时按快门,该操作定有不小误差。

    1.2 matlab 预实验

    初始使用 cv2.findChessboardCorners 并未检测到任何角点,而且检测时间漫长。因此先使用matlab的工具箱进行预实验(可参考链接3)。得到如下结果,说明数据还是可以被检测出来的。图中的数据已经经过筛选,去除了误差较大的几幅图像。观察其结果并查资料得知,该函数只能检测内角点,因此我的图像给的棋盘大小应为 (7,5) ,因此注意更改参数。

    1.3 python 预实验

    读取 left 和 right 两个文件夹内的所有图像,对其进行角点检测,记录下来能够检测到的图像,并获取能对应起来的图像对,最终可以得到4对图像(10对中才能获取4对,幸好师兄爱上镜)。下面代码中 truth 查看的是能够检测到角点的文件名。

    all_images = glob('./datas/*/*.jpg')
    
    truth = []
    def test1(image):
        img_l = cv2.imread(image)
        gray_l = cv2.cvtColor(img_l, cv2.COLOR_BGR2GRAY)
        ret_l, corners_l = cv2.findChessboardCorners(gray_l, (7, 5), None)
        if ret_l:
            truth.append(image)
    
    for i in all_images:
        print(f'-------------processing {i}--------------')
        test1(i)
    print(truth)
    

    2. 相机标定与参数求解

    2.1 确定坐标系以及坐标

    1. 定义世界坐标系,即三维坐标
    2. 对左右相机的两幅图像进行角点检测并记录其坐标
    3. 进行相机标定

    关于其原理简单来说就是世界坐标系、相机坐标系、图像坐标系的相互转化,转换之间就是用矩阵进行求解,如何求解就是通过一些计算方式,如通过对应点使用SVD等方式求解。

    标定的过程分为两部分:一是从世界坐标系转换为相机坐标系,二是从相机坐标系转为图像坐标系(可参考链接7)。如果知道相机坐标系中的一个点X(现实三维世界中的点),在像平面坐标系对应的点是X’,要求从相机坐标系转为像平面坐标系的转换,可以由一个变换矩阵得到,这个就是相机的内参和偏移量。从世界坐标系转换到相机坐标系是三维空间到三维空间的变换,一般来说需要一个平移操作和一个旋转操作就可以完成这个转换。

    畸变(distortion)是对直线投影(rectilinear projection)的一种偏移。简单来说直线投影是场景内的一条直线投影到图片上也保持为一条直线。畸变一般可以分为两大类,包括径向畸变和切向畸变。对于畸变校正,我们可以采用两个相机坐标系的而求得的变换矩阵来进行校正。

    通过该步骤我们可以得到二维相机坐标系下的坐标以及距离,包括两个相机内参,基本矩阵,旋转平移矩阵。

    求得参数结果如下:

    • M1: 第一个摄像机的相机矩阵
    • d1: 第一个摄像机的畸变向量
    • M2: 第二个摄像机的相机矩阵
    • d2: 第二个摄像机的畸变向量
    • R: 第一和第二个摄像机之间的旋转矩阵
    • T: 平移矩阵
    • E, F: 基本矩阵

    2.2 求解参数

    1. 通过步骤2.1中得到的参数,计算得到两个相机的内参、畸变向量和基本矩阵,两个相机间的旋转矩阵与平移矩阵
    2. 通过上一步求得的参数可以求得两个相机的校正变换矩阵,在新坐标系下的投影矩阵,以及深度差异映射矩阵

    在使用 remap 重构后,可能得到的图像是黑色的,这是因为双目校正函数 stereoRectify 中的参数 alpha 对矩阵的结果有影响,需要经过调试得到较高质量的重构图像。

    求得参数结果如下:

    • R1: 第一个摄像机的校正变换矩阵(旋转矩阵)
    • R2: 第二个摄像机的校正变换矩阵(旋转矩阵)
    • P1: 第一个摄像机在新坐标系下的投影矩阵
    • P2: 第二个摄像机在新坐标系下的投影矩阵
    • Q: 深度差异映射矩阵

    3. 视差图计算

    这部分直接借用了 链接3. 中的代码,讲解也比较清楚易懂。简单来说就是经过步骤2.2可以得到相机A与相机B之间相互转换的矩阵,以及投影到二维平面(新坐标系)的投影矩阵,通过矩阵乘法即可得到新的坐标,进行像素赋值即可。函数: cv2.remap

    得到的视差图可能是白色一片(投影特别浅),需要调节参数,我将其 sigma 调大了。

    # 视差图计算
    def depth_map(imgL, imgR, sigma=1.3):
        """ Depth map calculation. Works with SGBM and WLS. Need rectified images, returns depth map ( left to right disparity ) """
        # SGBM Parameters -----------------
        window_size = 3  # wsize default 3; 5; 7 for SGBM reduced size image; 15 for SGBM full size image (1300px and above); 5 Works nicely
    
        left_matcher = cv2.StereoSGBM_create(
            minDisparity=-1,
            numDisparities=5*16,  # max_disp has to be dividable by 16 f. E. HH 192, 256
            blockSize=window_size,
            P1=8 * 3 * window_size,
            # wsize default 3; 5; 7 for SGBM reduced size image; 15 for SGBM full size image (1300px and above); 5 Works nicely
            P2=32 * 3 * window_size,
            disp12MaxDiff=12,
            uniquenessRatio=10,
            speckleWindowSize=50,
            speckleRange=32,
            preFilterCap=63,
            mode=cv2.STEREO_SGBM_MODE_SGBM_3WAY
        )
        right_matcher = cv2.ximgproc.createRightMatcher(left_matcher)
        # FILTER Parameters
        lmbda = 80000
        visual_multiplier = 6
    
        wls_filter = cv2.ximgproc.createDisparityWLSFilter(matcher_left=left_matcher)
        wls_filter.setLambda(lmbda)
    
        wls_filter.setSigmaColor(sigma)
        displ = left_matcher.compute(imgL, imgR)  # .astype(np.float32)/16
        dispr = right_matcher.compute(imgR, imgL)  # .astype(np.float32)/16
        displ = np.int16(displ)
        dispr = np.int16(dispr)
        filteredImg = wls_filter.filter(displ, imgL, None, dispr)  # important to put "imgL" here!!!
    
        filteredImg = cv2.normalize(src=filteredImg, dst=filteredImg, beta=0, alpha=255, norm_type=cv2.NORM_MINMAX);
        filteredImg = np.uint8(filteredImg)
    
        return filteredImg
    

    将校对后每对图像合并,并画线查看其结果,与视差图一并显示。其中右图是视差图,仅截屏了一副图像,代码是将其都显示。


    4. 重构三维点

    该部分采用了 cv2.triangulatePoints 函数,triangulatePoints 函数得到3d点时,它们处于齐次坐标, points[x,y,z,w][x, y, z, w]ww 代表齐次的比值,其结果应为 [x/w,y/w,z/w][x/w, y/w, z/w]

    5. 参考链接

    1. matlab双目标定(详细过程)
    2. https://github.com/bvnayak/stereo_calibration (在该代码的基础上做了部分更改以及补充)
    3. 双目立体视觉 I:标定和校正 (代码借鉴较多的部分,值得学习)
    4. Python计算机视觉-张氏相机内参标定
    5. 【Python计算机视觉编程】第五章 多视几何
    6. 还有课程的PPT
    7. 相机标定(Camera calibration)(原理学习较多的部分)

    6. CODE

    直接放入计算过程的代码

    1. 基本函数
    import os
    import cv2
    import numpy as np
    from glob import glob
    import matplotlib.pyplot as plt
    
    class StereoCalibration():
        def __init__(self, filepath):
            
            self.criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 30, 0.001)
            
            # prepare object points, like (0,0,0), (1,0,0), (2,0,0) ....,(6,5,0)
            self.objp = np.zeros((7*5, 3), np.float32)
            self.objp[:, :2] = np.mgrid[0:7, 0:5].T.reshape(-1, 2)
    
            # Arrays to store object points and image points from all the images.
            self.objpoints = []  # 3d point in real world space
            self.imgpoints_l = []  # 2d points in image plane.
            self.imgpoints_r = []  # 2d points in image plane.
    
            self.filepath = filepath
            self.read_images()
    
        def read_images(self):
            images_right = glob(self.filepath + 'right/[4|6|8].jpg') + ['./datas/right/11.jpg']
            images_left = glob(self.filepath + 'left/[4|6|8].jpg') + ['./datas/left/11.jpg']
    
            for i, fname in enumerate(images_right):
                img_l = cv2.imread(images_left[i])
                img_r = cv2.imread(images_right[i])
    
                gray_l = cv2.cvtColor(img_l, cv2.COLOR_BGR2GRAY)
                gray_r = cv2.cvtColor(img_r, cv2.COLOR_BGR2GRAY)
    
                # Find the chess board corners
                ret_l, corners_l = cv2.findChessboardCorners(gray_l, (7, 5), None)
                ret_r, corners_r = cv2.findChessboardCorners(gray_r, (7, 5), None)
    
                # If found, add object points, image points (after refining them)
                
    
                if ret_l and ret_r:
                    
                    self.objpoints.append(self.objp)
                    
                    rt = cv2.cornerSubPix(gray_l, corners_l, (5, 5), (-1, -1), self.criteria)
                    self.imgpoints_l.append(corners_l)
    
                    rt = cv2.cornerSubPix(gray_r, corners_r, (5, 5), (-1, -1), self.criteria)
                    self.imgpoints_r.append(corners_r)
                else:
                    print('Couldn\'t be found')
    
                img_shape = gray_l.shape[::-1]
    
            rt, self.M1, self.d1, self.r1, self.t1 = cv2.calibrateCamera(self.objpoints, self.imgpoints_l, img_shape, None, None)
            rt, self.M2, self.d2, self.r2, self.t2 = cv2.calibrateCamera(self.objpoints, self.imgpoints_r, img_shape, None, None)
    
            self.camera_model = self.stereo_calibrate(img_shape)
        
        def stereo_calibrate(self, dims):
            
            h, w = dims
            
            flags = 0
            # 如果该标志被设置,那么就会固定输入的cameraMatrix和distCoeffs不变,只求解R,T,E,F.
            flags |= cv2.CALIB_FIX_INTRINSIC
            # 根据用户提供的cameraMatrix和distCoeffs为初始值开始迭代
            flags |= cv2.CALIB_USE_INTRINSIC_GUESS
            # 迭代过程中不会改变焦距
            flags |= cv2.CALIB_FIX_FOCAL_LENGTH
            # 切向畸变保持为零
            flags |= cv2.CALIB_ZERO_TANGENT_DIST
    
    
            stereocalib_criteria = (cv2.TERM_CRITERIA_MAX_ITER + cv2.TERM_CRITERIA_EPS, 100, 1e-5)
            ret, M1, d1, M2, d2, R, T, E, F = cv2.stereoCalibrate(
                self.objpoints, 
                self.imgpoints_l,
                self.imgpoints_r, 
                self.M1, 
                self.d1, 
                self.M2,
                self.d2, 
                dims,
                criteria=stereocalib_criteria,
                flags=flags
            )
    
            camera_model = dict([
                ('M1', M1), 
                ('M2', M2), 
                ('dist1', d1),
                ('dist2', d2), 
                ('rvecs1', self.r1),
                ('rvecs2', self.r2), 
                ('R', R), 
                ('T', T),
                ('E', E), 
                ('F', F), 
                ('dim', (h, w))
            ])
    
            return camera_model
        
        def rectify(self, camera_model):
            
            M1, d1, M2, d2, R, T =  camera_model.get('M1'), camera_model.get('d1'), camera_model.get('M2'), camera_model.get('d2'), camera_model.get('R'), camera_model.get('T')
            
            # 双目矫正 alpha=-1, 0, 0.9
            R1, R2, P1, P2, Q, validPixROI1, validPixROI2 = cv2.stereoRectify(M1, d1, M2, d2, camera_model.get('dim'), R, T, alpha=-1)
            print('stereo rectify done...')
    
            # 得到映射变换
            stereo_left_mapx, stereo_left_mapy = cv2.initUndistortRectifyMap(M1, d1, R1, P1, camera_model.get('dim'), cv2.CV_32FC1)
            stereo_right_mapx, stereo_right_mapy = cv2.initUndistortRectifyMap(M2, d2, R2, P2, camera_model.get('dim'), cv2.CV_32FC1)
            print('initUndistortRectifyMap done...')
            
            rectify_model = dict([
                ('R1', R1),
                ('R2', R2),
                ('P1', P1),
                ('P2', P2),
                ('Q', Q),
                ('stereo_left_mapx', stereo_left_mapx),
                ('stereo_left_mapy', stereo_left_mapy),
                ('stereo_right_mapx', stereo_right_mapx),
                ('stereo_right_mapy', stereo_right_mapy)
            ])
            
            return rectify_model
    
    # 视差图计算
    # 略,复制第3部分
    
    # 画线函数
    def drawLine(img, num=16):
        h, w, _ = img.shape
        for i in range(0, h, h // num):
            cv2.line(img, (0, i), (w, i), (0, 255, 0), 1, 8)
        return img
    
    1. 创建对象,得到基本参数
    filepath = 'X:/Notebook/HW/CV/datas/'
    cal_data = StereoCalibration(filepath)
    # 得到相机参数
    camera_model = cal_data.camera_model
    # 得到修正参数
    rectify_model = cal_data.rectify(camera_model)
    
    1. 可视化结果
    stereo_left_mapx, stereo_left_mapy = rectify_model.get('stereo_left_mapx'), rectify_model.get('stereo_left_mapy')
    stereo_right_mapx, stereo_right_mapy = rectify_model.get('stereo_right_mapx'), rectify_model.get('stereo_right_mapy')
    
    # 看网格是否对齐, 并且显示视差图
    plt.figure(figsize=(20,40))
    num = 1
    for i in range(len(left_images)):
        left_img = cv2.imread(left_images[i])
        right_img = cv2.imread(right_images[i])
        
        frame0 = cv2.remap(right_img, stereo_right_mapx, stereo_right_mapy, cv2.INTER_LINEAR, cv2.BORDER_CONSTANT)
        frame1 = cv2.remap(left_img, stereo_left_mapx, stereo_left_mapy, cv2.INTER_LINEAR, cv2.BORDER_CONSTANT)
        gray_left = cv2.cvtColor(frame0, cv2.COLOR_BGR2GRAY)
        gray_right = cv2.cvtColor(frame1, cv2.COLOR_BGR2GRAY)
        img = np.concatenate((frame0, frame1), axis=1).copy()
        img = drawLine(img, 32)
        plt.subplot(len(left_images), 2, num)
        plt.axis('off')
        plt.imshow(img[...,::-1])
        num += 1
        
        disparity_image = depth_map(gray_left, gray_right, sigma=2.0)
        plt.subplot(len(left_images), 2, num)
        plt.axis('off')
        plt.imshow(disparity_image, cmap='gray')
        num += 1
    #     break
    
    1. 重构
    pls = [np.array(i) for i in cal_data.imgpoints_l]
    prs = [np.array(i) for i in cal_data.imgpoints_r]
    
    # 取一对图像的二维坐标点
    ptl, ptr = pls[0], prs[0]
    P1 = rectify_model.get('P1')
    P2 = rectify_model.get('P2')
    
    points = cv2.triangulatePoints(P1,P2,ptl,ptr)
    coodinates_x = points[0] / points[-1]
    coodinates_y = points[1] / points[-1]
    coodinates_z = points[2] / points[-1]
    
    # 转换成三维坐标的形式,直接reshape结果不对
    coodinates = np.array([[i,j,k] for i,j,k in zip(coodinates_x, coodinates_y, coodinates_z)])
    
    fig = plt.figure(figsize=(12,12))
    ax = fig.add_subplot(111, projection='3d')
    
    plt.cla()
    coodinates_x, coodinates_y, coodinates_z
    ax.scatter(coodinates[:, 0], coodinates[:, 1], coodinates[:, 2], color='red')
    plt.draw()
    plt.show()
    
    展开全文
  • 基于双目视觉的深度计算和三维重建代码我自己用过,绝对没问题
  • 测量,三维重建 示例程序 参考文献 关于双目视觉的一些总结 笔者2013年进入吉林大学软件学院,2014年开始写自己的第一个完整的程序,期间受到过无数前辈的帮助,正是这个程序的完成给了我极大的信心,也让我喜欢上...

    关于双目视觉的一些总结

    笔者2013年进入吉林大学软件学院,2014年开始写自己的第一个完整的程序,期间受到过无数前辈的帮助,正是这个程序的完成给了我极大的信心,也让我喜欢上编程.这个程序是"基于OpenCV的双目测距",他的主要代码来自于邹宇华老师的OpenCV例程,我只不过进行了一些小小的修改,然后做了一个界面,在收获了许多经验的同时,也发现了双目视觉,乃至机器视觉的中文社区环境的一个问题,那就是有效资源非常少,很多人的博客只不过是互相转载,很多人提供的代码要么语焉不详,要么就是其他示例程序甚至是官方实例的修修补补,这对于初学者非常的不友好,一知半解的教程会让新手走很多弯路.在这期间,我的双目视觉程序仍然维持着更新,我决定把相关的资料全部公开,如果您发现错误,请及时留言评论,同时欢迎转载.文末附有程序代码链接,如果觉得有帮助,请加star,谢谢.

    双目视觉程序: https://github.com/jiafeng5513/Evision
    新版程序演示视频:https://www.bilibili.com/video/av46024738
    旧版程序演示视频:https://www.bilibili.com/video/av8862669

    相机模型

    通常关于相机模型的文章只会提到小孔成像模型,这里要提醒的是,小孔成像模型并不能代表所有种类的相机,为了彻底搞清楚计算机视觉的基础,我们必须把相机模型弄明白.
    1. 相机的投影模型
    图1.相机投影

    • 球面投影模型可以概括所有的相机模型,或者说所有的相机模型都是球面投影模型的一种特殊情况.
    • 我们可以在光心和成像平面之间假想一个球面,物体先在球面上成像,再从球面上投影到成像平面上.
    • 透镜的光学性质都是中心对称的,所以球面假设很合理.
    • 穿过光心的光线是直线传播的,不同模型的角度差异的来源在于球面向平面投影的不同方法.
    • 从球面向平面的投影,最常见的用处是绘制世界地图(可以百度:墨卡托投影,圆锥投影,等积投影等),
      从球面向平面的投影有很多种方法,相机在制造时,镜头由很多片透镜组成,这些透镜的左右就是逼近某一种投影方法,使得透过镜头的光线直接透射在传感器平面上的情况等效于先投影在成像球面,再向成像平面投影,物体与光轴的夹角θ变化时,像点与图像中心的距离R也随之变化,按照变化关系的不同,经典相机模型被分为5类,其中透视投影模型又称为小孔模型.
    计算公式 英文名 中文名
    R=f*tan(θ) Perspective 透射投影
    R=f*θ Equidistant 等距投影
    R=2fsin(θ/2) Equisolid angle 等立体角投影(等积投影)
    R=f*sin(θ) Orthographic 正交投影
    R=2ftan(θ/2) Stereographic 球极投影(体视投影,保角投影)

    图2 五种投影模型
    关于鱼眼镜头的投影特点,推荐大家看这篇知乎专栏,里面有一个视频直观的描述了这些投影的成像特点.

    2. 不同类型的镜头

    • 由于镜头制造上的限制,无法制造出满足理想情况的镜头,所以设计镜头时采用一个多项式逼近理想情况,其中K称为径向畸变系数.
    • 由于实际镜头在拟合理想投影模型略有不同,这种性质是相机镜头的一种属性,用K描述.
    • 径向畸变是最主要的畸变来源.此外,相机还有其他畸变,例如切向畸变,薄棱镜畸变,传感器倾斜畸变.
      在这里插入图片描述
      通过观察R和θ的曲线,我们能获知这些投影方式的一些性质:
      图4 R和θ的关系
      根据图像我们能得到一些结论:
    • 标准镜头(使用透射投影)的视角无法达到180°
    • 鱼眼镜头能够实现180°的视角.
    • 虽然理论上四种投影类型的鱼眼镜头都能实现180°视角,但是其最大视角对应的r是不同的,在传感器尺寸固定的情况下,不同的鱼眼相机实际能实现的最大视角不相同.
    • 标准镜头在被摄物体到光心与光轴之间的夹角接近90°的过程中,单位角度变化所对应的成像点位移迅速增大,可以预见,图片边缘的情况一定和图片中心有所区别.
      关于标准模型,这里提醒大家:
      在一些涉及到相机模型的场景,要关注其中使用的相机模型是否符合预期,小孔模型使用的比较多,但是在使用广角相机和鱼眼相机时,小孔模型未必是有效的,此外在有些场景下,小孔模型过多的畸变系数可能会引起麻烦(例如训练深度神经网络时,过多的参数可能会引入复杂的可训练性的证明).
      3. 小孔成像模型
      在这里插入图片描述
      在这里插入图片描述
    • (X,Y,Z) 是被摄物体的世界坐标
    • (u,v) 是投影点在像素坐标系下的坐标
    • A是相机矩阵,也叫内参矩阵,描述从相机坐标系到图像坐标系的变换
    • (cx,cy)是主点坐标,通常位于图像中心
    • fx,fy是以像素单位表示的焦距
    • S是缩放系数
    • [R|t] 也叫外参矩阵,描述由世界坐标系变换到相机坐标系

    外参数把世界坐标变换到相机坐标,内参数把相机坐标变换到图像坐标
    在这里插入图片描述
    考虑一些畸变:
    在这里插入图片描述
    其中k是径向畸变系数,P是切向畸变系数,S是薄棱镜畸变系数,在大多数情况下,只考虑径向和切向畸变的低阶系数即可,在实验精度有限的前提下,使用高阶系数可能会引入很大的误差.此外,薄棱镜畸变和倾斜传感器畸变在大多数相机上都非常小
    在这里插入图片描述
    径向畸变K1的正负会影响畸变的视觉效果,负数时一般称为桶形畸变,正数时称为枕形畸变

    4. 极线约束,本质矩阵和基本矩阵
    在这里插入图片描述

    • 相机从两个位姿O1和O2拍摄同一个点P,分别成像p1和p2,这里p1和p2是相机坐标系下的坐标.PO1O2共面,两个位姿的成像面为
    • I1I1I1I1I1 I_1Pbc(T×(RPac))=0

    5. 单应矩阵(Homography matrix)

    • Homography matrix表达两个平面之间的投影关系,例如从三维空间中的实际平面到像平面,或者一个像平面到另一个像平面.
    • 两张图片上的对应点的齐次坐标分别是(x’,y’,1)和(x,y,1),则有以下等式,其中H为单应矩阵:
    标定方法 优点 缺点 常用方法
    传统相机标定法 可使用于任意的相机模型、 精度高 需要标定物 Tsai两步法、张氏标定法
    主动视觉相机标定法 不需要标定物、算法简单、鲁棒性高 成本高、设备昂贵 主动系统控制相机做特定运动
    相机自标定法 灵活性强、可在线标定 精度低、鲁棒性差 分层逐步标定、基于Kruppa方程

    最常用的标定方法是张氏标定,下面简要介绍张氏标定的原理:
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

    视差算法:立体匹配

    立体匹配算法也称为视差算法,目的是求解两幅图片之间的视差值.首先我们来定义视差值这个概念:

    定义1:三维世界中有某一定点P,相机在两个位置对P成像,分别获取像点p1p1p1p1p1 p_1p2同名点
    定义2:如果两张图片经过了畸变矫正,并且同名点的Y坐标在各自的像素坐标系下相等,我们称这两张图片为已校正图片对
    定义3:在双目立体视觉中,两个相机的视角方向朝大致相同的方向,此时观察者站在相机的后方,视线方向与两台相机一致,他左手边的相机我们称之为 左相机,拍摄的图片称为 左视图
    定义4:已校正图片对上,同名点的是X坐标之差,为视差值,逐像素计算式差值后组成的图片,称为 视差图 .特别地,在双目立体视觉中,以左视图为基准,通过同名点在右图上的坐标减去其在左图上的坐标得到的,称为左视差(图),相对应的称为右视差图

    通过定义我们发现,立体匹配问题实际上是求同名点问题,因为只要找到准确的同名点,视差值不过是一个坐标差.可以说绝大多数视差算法的核心都是寻找同名点(一些基于深度学习的视差方法并不显式的求同名点).更具体点来说,当要求左视差图的时候,就是从左视图中依次取出每个像素,然后到右视图的某个范围内找 这个像素的同名点,这样当我们找完左视图上所有像素的同名点时,就得到了一张和左视图一样大的视差图了.

    上面是通俗化的解释,接下来我们进行形式化的说明.视差算法是一个机器视觉的传统领域,经过多年的发展已经相当成熟,根据Schrstein和Szeliski的总结,双目立体匹配可划分为四个步骤:匹配代价计算、代价聚合、视差计算和视差优化。

    1. 匹配代价计算
      匹配代价计算的目的是衡量待匹配像素与候选像素之间的相关性。两个像素无论是否为同名点,都可以通过匹配代价函数计算匹配代价,代价越小则说明相关性越大,是同名点的概率也越大。

      每个像素在搜索同名点之前,往往会指定一个视差搜索范围D(Dmin ~ Dmax),视差搜索时将范围限定在D内,用一个大小为W×H×D(W为影像宽度,H为影像高度)的三维矩阵C来存储每个像素在视差范围内每个视差下的匹配代价值。矩阵C通常称为DSI(Disparity Space Image)。

      匹配代价计算的方法有很多,传统的摄影测量中,使用灰度绝对值差(AD,Absolute Differences)[1]、灰度绝对值差之和(SAD,Sum of Absolute Differences)、归一化相关系数(NCC,Normalized Cross-correlation)等方法来计算两个像素的匹配代价;计算机视觉中,多使用互信息(MI,Mutual Information)法[2,3]、Census变换(CT,Census Transform)法[4,5]、Rank变换(RT, Rank Transform)法[6,7]、BT(Birchfield and Tomasi)法[8]等作为匹配代价的计算方法。不同的代价计算算法都有各自的特点,对各类数据的表现也不尽相同,选择合适的匹配代价计算函数是立体匹配中不可忽视的关键步骤。

      无论是什么样的代价计算方法,都可以近似的看作这样的一个过程:首先在左图中确定一个像素或者一块像素区域,然后在右图中用某种形状的模板或者策略选定一块匹配区域,这个区域的大小和滑动步长都是超参数控制的,然后计算代价,移动模板窗口,再计算代价值,直到耗尽搜索范围,然后左图中选定第二块或者第二个像素,再重复这个过程。

    2. 代价聚合
      在这里插入图片描述
      代价聚合的根本目的是让代价值能够准确的反映像素之间的相关性。上一步匹配代价的计算往往只会考虑局部信息,通过两个像素邻域内一定大小的窗口内的像素信息来计算代价值,这很容易受到影像噪声的影响,而且当影像处于弱纹理或重复纹理区域,这个代价值极有可能无法准确的反映像素之间的相关性,直接表现就是真实同名点的代价值非最小。

      而代价聚合则是建立邻接像素之间的联系,以一定的准则,如相邻像素应该具有连续的视差值,来对代价矩阵进行优化,这种优化往往是全局的,每个像素在某个视差下的新代价值都会根据其相邻像素在同一视差值或者附近视差值下的代价值来重新计算,得到新的DSI,用矩阵S来表示。

      实际上代价聚合类似于一种视差传播步骤,信噪比高的区域匹配效果好,初始代价能够很好的反映相关性,可以更准确的得到最优视差值,通过代价聚合传播至信噪比低、匹配效果不好的区域,最终使所有影像的代价值都能够准确反映真实相关性。常用的代价聚合方法有扫描线法、动态规划法、SGM算法中的路径聚合法等。

    3. 视差计算
      视差计算即通过代价聚合之后的代价矩阵S来确定每个像素的最优视差值,例如使用赢家通吃算法(WTA,Winner-Takes-All)来计算,如图2所示,即某个像素的所有视差下的代价值中,选择最小代价值所对应的视差作为最优视差。这一步非常简单,这意味着聚合代价矩阵S的值必须能够准确的反映像素之间的相关性,也表明上一步的代价聚合步骤是立体匹配中极为关键的步骤,直接决定了算法的准确性。
      在这里插入图片描述

    4. 视差优化
      视差优化的目的是对上一步得到的视差图进行进一步优化,是一些后处理方法的组合。改善视差图的质量,包括剔除错误视差、适当平滑以及子像素精度优化等步骤,一般采用左右一致性检查(Left-Right Check)算法剔除因为遮挡和噪声而导致的错误视差;采用剔除小连通区域算法来剔除孤立异常点;采用中值滤波(Median Filter)、双边滤波(Bilateral Filter)等平滑算法对视差图进行平滑;另外还有一些有效提高视差图质量的方法如鲁棒平面拟合(Robust Plane Fitting)、亮度一致性约束(Intensity Consistent)、局部一致性约束(Locally Consistent)等也常被使用。

      由于WTA算法所得到的视差值是整像素精度,为了获得更高的子像素精度,需要对视差值进行进一步的子像素细化,常用的子像素细化方法是一元二次曲线拟合法,通过最优视差下的代价值以及左右两个视差下的代价值拟合一条一元二次曲线,取二次曲线的极小值点所代表的视差值为子像素视差值。如图3所示。

      在编程实践中这里实际上最容易出现问题。OpenCV实现的BM和SGBM并没有亚像素增强的步骤,所以它输出的视差数据是CV8U也就是unsingned char类型的,他的视差等级非常有限,比如说设置的最大视差是50,那他的视差就是0-49这50个整数,精度就很差了;如果进行亚像素增强,输出的数据格式就是float或者double的,这是不能直接在屏幕上显示的,必须要进行归一化,所以如果视差数据只是整个程序的中间结果,那么就不可以保存成图片再读取,这样就会损失精度。

    5. 代价函数
      代价函数的想法是,确定了同名点实际上视差就确定了,那么我们输入左点和视差,在给定的视差范围内计算代价值,代价最低的那个情况就最可能是同名点,此时的视差也最可信.下面列举几种常见的代价函数:
      在这里插入图片描述
      在这里插入图片描述
      在这里插入图片描述
      在这里插入图片描述
      在这里插入图片描述

    6. OpenCV-BM
      OpenCV-BM算法使用SAD作为代价函数,在源码中我们看到该方法被称为“基于SAD的快速立体匹配算法”.根据对源码的解析,我把OpenCV-BM的流程分为四个部分:校验、预处理、视差计算和后处理.

      校验:校验算法的输入和超参数的选取是否合理.左右视图大小相等,格式均为CV_8UC1,接收视差图的Mat为CV_16SC1或者CV_32FC1, preFilterSize是[5,255]之间的奇数, preFilterCap在[1,63], SADWindowSize是[5,255]之间的奇数且不能大于图片的宽或高, numDisparities必须能被16整除, textureThreshold必须是非负数, uniquenessRatio必须是非负数.

    超参数 说明
    1 preFilterType 前处理可以是NORMALIZED或XSOBEL
    2 preFilterSize 归一化窗口大小,对于XSobel没有意义
    3 preFilterCap 截断值.对于XSobel(X方向上的Sobel),如果滤波后结果在[-prefilterCap,preFilterCap]之间,对应取[0, 2*preFilterCap]
    4 SADWindowSize SAD窗口大小
    5 minDisparity 视差搜索起点
    6 numDisparities 视差窗口,一个推荐值是((width / 8) + 15) & (~0xfl)
    7 textureThreshold 低纹理区域的判断阈值. 如果当前 SAD 窗口内所有邻居像素点的 x-导数绝对值之和小于指定阈值,则意味着当前窗口内的图像变化率太低也就是纹理很差,则该窗口对应的像素点的视差值设为0.
    8 uniquenessRatio 视差唯一性百分比. 视差窗口范围内最低代价是次低代价的 (1 + uniquenessRatio / 100) 倍时,最低代价对应的视差值才是该像素点的视差, 否则该像素点的视差为 0
    9 speckleRange 散斑窗口内允许的最大波动值
    10 speckleWindowSize 散斑滤波窗口大小,<=0时不进行散斑滤波
    11 disp12MaxDiff 左右视差图的最大容许差异,超标清零,左右一致性检查的目的是找到遮挡点,OpenCV默认的方法直接去除掉了遮挡点的视差,小于0时跳过
    12 dispType 视差图数据类型

    上表是OpenCV-BM的超参数,可以看到OpenCV实现的BM算法有很多地方还是很暴力的,它的前处理和后处理都很简单,再后来的很多论文中着墨甚多的后处理步骤例如遮挡点和miss点的处理在OpenCV中都极其简陋.然而OpenCV的源码水平非常高,对于内存的使用和优化,变量类型的使用等问题绝不含糊,非常值得学习.有人对源码进行了注释[26].
    预处理:X方向的sobel边缘检测和归一化
    视差计算:SAD带价函数
    后处理:左右视差一致性检查,散斑滤波

    1. OpenCV-SGBM
      也称为SGM、半全局匹配算法等[9]。
      观察OpenCV-BM的实现可知,在计算某一点的视差(代价)时,只考虑一个窗口内的像素,窗口外面的图像还有很多像素,但是并不会参与到本点的视差计算当中.事实上有一些视差算法会考虑所有像素的影响,这种算法就叫全局算法,BM是局部算法,而SGM是半全局的.需要注意的是,OpenCV-SGBM(以下简称SGBM)对于原有的SGM(以下简称SGM)做出了一些调整.

      首先看SGM算法.SGM使用分层互信息代价函数,基于动态规划的代价聚合,与OpenCV-BM类似的前处理和后处理.互信息的计算如下:
      在这里插入图片描述
      其中, MII1,I2MII1,I2MII1,I2MII1,I2MII1,I2 MI_{I_1,I_2}
    展开全文
  • 这里特别感谢博主shiter的原创文章:OpenCV+OpenGL 双目立体视觉三维重建 本博文参考了该博主的的核心代码,并针对该博主博文中声明的一些BUG进行了修正: 本文代码下载地址(已修正相关问题问题):...
  • 机器视觉、ZNCC、视差图计算,含图片。可以直接使用的代码。matlab代码。希望能帮到大家。
  • opencv1写的双目视觉摄像机标定和三维重建代码——opencv1写的双目视觉摄像机标定和三维重建代码
  • 文件名称: correlCorresp下载 收藏√ [5 4 3 2 1]开发工具: matlab文件大小: 759 KB上传时间: 2015-05-26下载次数: 67提 供 者: 何苗苗详细说明:matlab 实现双目视觉三维重建 利用两张图片重建三维信息-matlab ...
  • 双目立体视觉代码

    热门讨论 2012-02-11 14:08:48
    双目立体视觉代码,包括标定,匹配,三维重建
  • 这个是计算机视觉三维重建的领域。该代码实现特征点提取和立体匹配的功能。三维重建模拟人眼,进行双目拍摄。这里在经典特征点特区SURF算法基础上,加入极线约束的思想,去除噪声和匹配错误的杂点。有很好的匹配效果...
  •  目前获取三维数据的方法很多,比如:双目视觉技术,深度相机等等。  但是存在一个问题:数据量巨大!如果要做在线检测,那速度很可能达不到要求。 总思路:  用二维的数据来表示三维数据(数据量从n3降到n2)。
  • 8.结构光:双目视觉(基于视差)

    千次阅读 2017-04-13 21:33:34
    Matlab基于视差进行三维重建代码如下: %% % 清理空间 clc; clear; close all; %% % 导入立体标定参数 load stereoParams.mat % 立体参数的可视化 figure; showExtrinsics(stereoParams); %% % 导入数据 ...
  • 成立近两年来,工坊深挖3D视觉的各个领域,主要涉及计算机视觉与深度学习、点云处理、SLAM、三维重建、结构光、双目视觉、深度估计、3D检测、自动驾驶、多传感器融合等,在校的童鞋和已经工作...
  • 大家好,计算机视觉life公众号推出的视频课程《全网最详细的ORB-SLAM2精讲:原理推导+逐行代码分析》正式上线啦!...适用于自动驾驶、增强现实、机器人、三维重建等领域。 我们的系列课程以详细注释的
  • 最近想要实现基于双目视觉三维重建,其中一些重要的函数,在opencv3.0及其以后的版本中才有,而且还不在Main modules中,是在Extra modules中。即opencv3的版本,分为两部分,比如opencv3.1.0标准版和opencv_...
  • 理论知识参考于基于图像的三维重建——对极几何 在这里用的opencv是440版本的,接下来就是相应的代码 /* 测试8点法求取基础矩阵F * * [直接线性变换法] * 双目视觉中相机之间存在对极约束 * * p2'Fp1=0, * * ...

空空如也

空空如也

1 2
收藏数 33
精华内容 13
关键字:

双目视觉三维重建代码