仅供参考,还未运行程序,理解部分有误,请参考英文原版。html
绿色部分非文章内容,是我的理解。算法
转载请注明:http://blog.csdn.net/raby_gyl/article/details/17471617编程
在这一章,咱们将讨论来至运动结构(Structure from Motion,SfM)的概念,或者从一个运动的相机拍摄到的图像中更好的推测提取出来的几何结构,使用OpenCV的API函数能够帮助咱们完成这个任务。首先,让咱们将咱们使用的冗长的方法约束为使用单个摄像机,一般称为单目方法,而且是一组分离的和稀疏的视频帧而不是连续的视频流。这两个约束在很大程度上简化了系统,这个系统咱们将在接下来的页码中进行描述,而且帮助咱们理解任何SfM方法的原理。为了实现咱们的方法,咱们将跟随Hartley和Zisserman的脚步(后面称做H和Z),伴随着他们有创意的书——《计算机视觉中的多视觉几何》的第9章到第12章的记录。数组
在这一章,咱们将涉及到如下内容:数据结构
一、来至运动结构的概念app
二、从一对图像中估计摄像机的运动框架
三、重构场景iview
四、从多个视图中重构机器学习
五、重构提纯(Refinement of the reconstruction)函数
六、可视化3D点云
本章自始自终假定使用一个标定过的相机——预先标定的相机。在计算机视觉中标定是一个广泛存在的操做,在OpenCV中获得了很好的支持,咱们可使用在前一章讨论的命令行工具来完成标定。所以,咱们假定摄像机内参数存在,而且具体化到K矩阵中,K矩阵为摄像机标定过程的一个结果输出。
为了在语言方面表达的清晰,从如今开始咱们将单个摄像机称为场景的单个视图而不是光学和获取图像的硬件。一个摄像机在空间中有一个位置,一个观察的方向。在两个摄像机中,有一个平移成分(空间运动)和一个观察方向的旋转。咱们一样对场景中的点统一一下术语,世界(world),现实(real),或者3D,意思是同样的,即咱们真实世界存在的一个点。这一样适用于图像中的点或者2D,这些点在图像坐标系中,是在这个位置和时间上一些现实的3D点投影到摄像机传感器上造成的。
在这一章的代码部分,你将注意参考《计算视觉中的多视觉几何》(Multiple View Geometry in Computer Vision),例如,//HZ9.12,这指的是这本书的第9章第12个等式。一样的,文本仅包括代码的摘录,而完整的代码包含在伴随着这本书的材料中。
第一个咱们应当区别的是立体(Stereo or indeed any multiview),使用标准平台的三维重建和SfM之间的差别。在两个或多个摄像机的平台中,假定咱们已经知道了两个摄像机之间的运动,而在SfM中,咱们实际上不知道这个运动而且咱们但愿找到它。标准平台,来至观察的一个过度简单的点,能够获得一个更加精确的3D几何的重构,由于在估计多个摄像机间距离和旋转时没有偏差(距离和旋转已知)。实现一个SfM系统的第一步是找到相机之间的运动。OpenCV能够帮助咱们在许多方式上得到这个运动,特别地,使用findFundamentalMat函数。
让咱们想一下选择一个SfM算法背后的目的。在不少状况下,咱们但愿得到场景的几何。例如,目标相对于相机的位置和他们的形状是什么。假定咱们知道了捕获同一场景的摄像机之间的运动,从观察的一个合理的类似点,如今咱们想重构这个几何。在计算视觉术语中称为三角测量(triangulation),而且有不少方法能够作这件事。它能够经过 射线相交的方法完成,这里咱们构造两个射线:来至于每一个摄像机投影中心的点和每个图像平面上的点。理想上,这两个射线在空间中将相交于真实世界的一个3D点(这个3D点在每一个摄像机中成像),以下面图表展现:
实际上,射线相交很是不可靠。H和Z推荐不使用它。这是由于两个射线一般不能相交,让咱们回到使用链接两个射线的最短线段上的中间点。相反,H和Z建议一些方法来三角化3D点(trianglulate 3D points,三角化3D点就是计算3D点的坐标,能够经过后面的内容加以理解),这些方法中咱们将在重构场景部分讨论两个。OpenCV如今的版本没有包含三角测量(triangulation)的API,所以,这部分咱们将本身编写代码。
学习完如何从两个视图恢复3D几何以后,咱们将看到咱们是怎么样加入更多的同一个场景的视图来得到更丰富的重构。在那时,大部分的SfM方法试图依靠束调整(Bundle Adjustment)来优化咱们的摄像机和3D点一束估计位置(the bundle of estimated positon),这部份内容在重构提纯部分(Refinement of the reconstruction section)。OpenCV在新的图像拼接的工具箱内包含了用来束调整的方法。然而,使OpenCV和C++完美结合的工做是丰富的外部工具,这些工具能够很容易地整合到一块儿。所以,咱们将看到如何如何整合一个外部的束调节器——灵巧的SSBA库。
既然咱们已经描述了使用OpenCV实现咱们的SfM方法的一个归纳,咱们将看到每一个分红是如何实现的。
事实上,在咱们开始着手两个摄像机之间的运动以前,让咱们检查一下咱们手边用来执行这个操做的输入和工具。首先,咱们有来至(但愿并非很是)不一样空间位置的同一场景的两个图像。这是一个强大的资产,而且咱们将确保使用它。如今工具都有了,咱们应当看一下数学对象,来对咱们的图像,相机和场景施加约束。
两个很是有用的数学对象是基础矩阵(用F表示)和本征矩阵(用E表示)。除了本征矩阵是假设使用的标定的相机,它们很是类似。这种状况针对咱们,因此咱们将选着它。OpenCV函数仅容许咱们经过findFundamentalMat函数找到基础矩阵。然而,咱们很是简单地使用标定矩阵(calibration matrix)K从本征矩阵中得到基础矩阵,以下:
Mat_<double> E = K.t() * F * K; //according to HZ (9.12)
本征矩阵,是一个3×3大小的矩阵,使用x’Ex=0在图像中的一点和另一个图像中的一点之间施加了一个约束,这里x是图像一中的的一点,x’是图像二中与之相对应的一点。这很是有用,由于咱们将要看到。咱们使用的另外一个重要的事实是本征矩阵是咱们用来为咱们的图像恢复两个相机的全部须要,尽管只有尺度,可是咱们之后会获得。所以,若是咱们得到了本征矩阵,咱们知道每个相机在空间中的位置,而且知道它们的观察方向,若是咱们有足够这样的约束等式,那么咱们能够简单地计算出这个矩阵。简单的由于每个等式能够用来解决矩阵的一小部分。事实上,OpenCV容许咱们仅使用7个点对来计算它,可是咱们但愿得到更多的点对来获得一个鲁棒性的解。
如今咱们将使用咱们的约束等式来计算本征矩阵。为了得到咱们的约束,记住对于图像A中的每个点咱们必须在图像B中找到一个与之相对应的点。咱们怎样完成这个匹配呢?简单地经过使用OpenCV的普遍的特征匹配框架,这个框架在过去的几年了已经很是成熟。
在计算机视觉中,特征提取和描述子匹配是一个基础的过程,而且用在许多方法中来执行各类各样的操做。例如,检测图像中一个目标的位置和方向,或者经过给出一个查询图像在大数据图像中找到类似的图像。从本质上讲,提取意味着在图像中选择点,使得得到好的特征,而且为它们计算一个描述子。一个描述子是含有多个数据的向量,用来描述在一个图像中围绕着特征点的周围环境。不一样的方法有不一样的长度和数据类型来表示描述子矢量。匹配是使用它的描述子从另一个图像中找到一组与之对应的特征。OpenCV提供了很是简单和有效的方法支持特征提取和匹配。关于特征匹配的更多信息能够在Chapter 3少许(无)标记加强现实中找到。
让咱们检查一个很是简单特征提取和匹配方案:
// detectingkeypoints SurfFeatureDetectordetector(); vector<KeyPoint> keypoints1, keypoints2; detector.detect(img1, keypoints1); detector.detect(img2, keypoints2); // computing descriptors SurfDescriptorExtractor extractor; Mat descriptors1, descriptors2; extractor.compute(img1, keypoints1, descriptors1); extractor.compute(img2, keypoints2, descriptors2); // matching descriptors BruteForceMatcher<L2<float>> matcher; vector<DMatch> matches; matcher.match(descriptors1, descriptors2, matches);
你可能已经看到相似的OpenCV代码,可是让咱们快速的回顾一下它。咱们的目的是得到三个元素:两个图像的特征点,特征点的描述子,和两个特征集的匹配。OpenCV提供了一组特征检测器、描述子提取器和匹配器。在这个简单例子中,咱们使用SurfFeatureDetector函数来得到SURF(Speeded-Up-Robust Features)特征的2D位置,而且使用SurfDescriptorExtractor函数来得到SURF描述子。咱们使用一个brute-force匹配器来得到匹配,这是一个最直接方式来匹配两个特征集,该方法是经过比较第一个集中的每个特征和第一个集的每个特征而且得到最好的匹配来实现的。
在下一个图像中,咱们将看到两个图像上的特征点的匹配:(这两个图像来至于the Fountain-P11 sequence,能够在网址中找到:http://cvlab.epfl.ch/~strecha/multiview/denseMVS.html.)
实际上,像咱们执行的原始匹配(raw matching),只有到达某一程度时执行效果才比较好而且许多匹配多是错误的。所以,大多数SfM方法对原始匹配进行一些滤波方式来确保正确和减小错误。一种滤波方式叫作交叉检验滤波,它内置于OpenCV的brute-force匹配器中。也就是说,若是第一个图像的一个特征匹配第二个图像的一个特征,而且经过反向检查,第二个图像的特征也和第一个图像的特征匹配,那么这个匹配被认为是正确的。另外一个常见的滤波机制(该机制用在提供的代码中)是基于同一场景的两个图像而且在他们之间存在某一种立体视觉关系,这样的一个事实基础之上来滤波的。在实践中,这个滤波器尝试采用鲁棒性的算法来计算这个基础矩阵,咱们将会在寻找相机矩阵(Finding camera matrices)部分学习这种计算方法,而且保留由该计算获得的带有小偏差的特征对。
一个替代使用丰富特征(如SURF)匹配的方法,是使用optical flow(OF光流)进行匹配。下面的信息框提供了光流的一个简短的概述。最近的OpenCV为从两个图像得到流场扩展了API,而且如今更快,更强大。咱们将尝试使用它做为匹配特征的一个替代品。
【注释:
光流是匹配来至一幅图像选择的点到另一幅图像选着点的过程,假定这两个图像是一个视频序列的一部分而且它们彼此很是相近。大多数的光流方法比较一个小的区域,称为搜索窗口或者块,这些块围绕着图像A中的每一点和一样区域的图像B中的每一点。遵循计算机视觉中一个很是普通的规则,称为亮度恒定约束(brightness constancy constraint)(和其余名字),图像中的这些小块从一个图像到另一个图像不会有太大的变化,所以,他们的幅值差接近于0。除了匹配块,更新的光流方法使用一些额外的方法来得到更好的结果。其中一个方法就是使用图像金字塔,它是图像愈来愈小的尺寸(大小)版本,这考虑到了工做的从粗糙到精致——计算机视觉中一个很是有用的技巧。另一个方法是定义一个流场上的全局约束,假定这些点相互靠近,向同一方向一块儿运动。在OpenCV中,一个更加深刻的光流方法能够在Chapter Developing Fluid Wall Using the Microsoft Kinect中找到,这个书能够在出版社网站上访问到。】
在OpenCV中,使用光流至关的简单,能够经过调用calcOpticalFlowPyrLK函数来实现。咱们想要保存光流的结果匹配,像使用丰富的特征那样,在将来,咱们但愿两种方法能够互换。为了这个目标,咱们必须安装一个特殊的匹配方法——能够与基于特征的方法互换,这将在下面的代码中看到:
Vector<KeyPoint>left_keypoints,right_keypoints; // Detect keypoints in the left and right images FastFeatureDetectorffd; ffd.detect(img1, left_keypoints); ffd.detect(img2, right_keypoints); vector<Point2f>left_points; KeyPointsToPoints(left_keypoints,left_points); vector<Point2f>right_points(left_points.size()); // making sure images are grayscale Mat prevgray,gray; if (img1.channels() == 3) { cvtColor(img1,prevgray,CV_RGB2GRAY); cvtColor(img2,gray,CV_RGB2GRAY); } else { prevgray = img1; gray = img2; } // Calculate the optical flow field: // how each left_point moved across the 2 images vector<uchar>vstatus; vector<float>verror; calcOpticalFlowPyrLK(prevgray, gray, left_points, right_points, vstatus, verror); // First, filter out the points with high error vector<Point2f>right_points_to_find; vector<int>right_points_to_find_back_index; for (unsigned inti=0; i<vstatus.size(); i++) { if (vstatus[i] &&verror[i] < 12.0) { // Keep the original index of the point in the // optical flow array, for future use right_points_to_find_back_index.push_back(i); // Keep the feature point itself right_points_to_find.push_back(j_pts[i]); } else { vstatus[i] = 0; // a bad flow } } // for each right_point see which detected feature it belongs to Mat right_points_to_find_flat = Mat(right_points_to_find). reshape(1,to_find.size()); //flatten array vector<Point2f>right_features; // detected features KeyPointsToPoints(right_keypoints,right_features); Mat right_features_flat = Mat(right_features).reshape(1,right_ features.size()); // Look around each OF point in the right image // for any features that were detected in its area // and make a match. BFMatchermatcher(CV_L2); vector<vector<DMatch>>nearest_neighbors; matcher.radiusMatch( right_points_to_find_flat, right_features_flat, nearest_neighbors, 2.0f); // Check that the found neighbors are unique (throw away neighbors // that are too close together, as they may be confusing) std::set<int>found_in_right_points; // for duplicate prevention for(inti=0;i<nearest_neighbors.size();i++) { DMatch _m; if(nearest_neighbors[i].size()==1) { _m = nearest_neighbors[i][0]; // only one neighbor } else if(nearest_neighbors[i].size()>1) { // 2 neighbors – check how close they are double ratio = nearest_neighbors[i][0].distance / nearest_neighbors[i][1].distance; if(ratio < 0.7) { // not too close // take the closest (first) one _m = nearest_neighbors[i][0]; } else { // too close – we cannot tell which is better continue; // did not pass ratio test – throw away } } else { continue; // no neighbors... :( } // prevent duplicates if (found_in_right_points.find(_m.trainIdx) == found_in_right_points. end()) { // The found neighbor was not yet used: // We should match it with the original indexing // ofthe left point _m.queryIdx = right_points_to_find_back_index[_m.queryIdx]; matches->push_back(_m); // add this match found_in_right_points.insert(_m.trainIdx); } } cout<<"pruned "<< matches->size() <<" / "<<nearest_neighbors.size() <<" matches"<<endl;
函数KeyPointsToPoints和PointsToKeyPoints是用来进行cv::Point2f和cv::KeyPoint结构体之间相互转换的简单方便的函数。从先前的代码片断咱们能够看到一些有趣的事情。第一个注意的事情是,当咱们使用光流时,咱们的结果代表,一个特征从左手边图像的一个位置移动到右手边图像的另一个位置。可是咱们有一组在右手边图像中检测到的新的特征,在光流中从这个图像到左手边图像的特征不必定是对齐的。咱们必须使它们对齐。为了找到这些丢失的特征,咱们使用一个k邻近(KNN)半径搜索,这给出了咱们两个特征,即感兴趣的点落入了2个像素半径范围内。
咱们能够看获得另一个事情是用来测试KNN的比例的实现,在SfM中这是一种常见的减小错误的方法。实质上,当咱们对左手边图像上的一个特征和右手边图像上的一个特征进行匹配时,它做为一个滤波器,用来移除混淆的匹配。若是右手边图像中两个特征太靠近,或者它们之间这个比例(the rate)太大(接近于1.0),咱们认为它们混淆了而且不使用它们。咱们也安装一个双重防护滤波器来进一步修剪匹配。
下面的图像显示了从一幅图像到另一幅图像的流场。左手边图像中的分成色箭头表示了块从左手边图像到右手边图像的运动。在左边的第二个图像中,咱们看到一个放大了的小的流场区域。粉红色箭头再一次代表了块的运动,而且咱们能够经过右手边的两个原图像的片断看到它的意义。在左手边图像上可视的特征正在向左移动穿过图像,粉红色箭头的方向在下面的图像中展现:
使用光流法代替丰富特征的优势是这个过程一般更快而且能够适应更多的点,使重构更加稠密。在许多光流方法中也有一个块总体运动的统一模型,在这个模型中,丰富的特征匹配一般不考虑。使用光流要注意的是对于从同一个硬件获取的连续图像,它处理的很快,然而丰富的特征一般不可知。它们之间的差别源于这样的一个事实:光流法一般使用很是基础的特征,像围绕着一个关键点的图像块,然而,高阶丰富特征(例如,SURF)考虑每个特征点的较高层次的信息。使用光流或者丰富的特征是设计师根据应用程序的输入所作的决定。
既然咱们得到了两个关键点之间的匹配,咱们能够计算基础矩阵而且从中得到本征矩阵。然而,咱们必须首先调整咱们的匹配点到两个数组,其中数组中的索引对应于另一个数组中一样的索引。这须要经过findFundametalMat函数得到。咱们可能也须要转换KeyPoint 结构到Point2f结构。咱们必须特别注意DMatch的queryIdx和trainIdx成员变量,在OpenCV中DMatch是容纳两个关键点匹配的结构,由于它们必须匹配咱们使用matche.match()函数的方式。下面的代码部分展现了如何调整一个匹配到两个相应的二维点集,以及如何使用它们来找到基础矩阵:
vector<Point2f>imgpts1,imgpts2; for( unsigned inti = 0; i<matches.size(); i++ ) { // queryIdx is the "left" image imgpts1.push_back(keypoints1[matches[i].queryIdx].pt); // trainIdx is the "right" image imgpts2.push_back(keypoints2[matches[i].trainIdx].pt); } Mat F = findFundamentalMat(imgpts1, imgpts2, FM_RANSAC, 0.1, 0.99, status); Mat_<double> E = K.t() * F * K; //according to HZ (9.12)
稍后咱们可能使用status二值向量来修剪这些点,使这些点和恢复的基础矩阵匹配。看下面的图像:用来阐述使用基础矩阵修剪后的点匹配。红色箭头表示特征匹配在寻找基础矩阵的过程当中被移除了,绿色箭头表示被保留的特征匹配。
如今咱们已经准备好寻找相机矩阵。这个过程在H和Z书的第九章进行了详细的描述。然而,咱们将使用一个直接和简单的方法来实现它,而且OpenCV很容易的为咱们作这件事。可是首先,咱们将简短地检查咱们将要使用的相机矩阵的结构:
咱们的相机使用该模型,它由两个成分组成,旋转(表示为R)和平移(表示为t)。关于它的一个有趣的事情是它容纳一个很是基本的等式:x=PX,这里x是图像上的二维点,X是空间上的三维点。还有更多,可是这个矩阵给咱们一个很是重要的关系,即图像中点和场景中点之间的关系。所以,既然咱们有了寻找摄像机矩阵的动机,那么咱们将看到这个过程是怎么完成的。下面的的代码部分展现了如何将本征矩阵分解为旋转和平移成分。
SVD svd(E); Matx33d W(0,-1,0,//HZ 9.13 1,0,0, 0,0,1); Mat_<double> R = svd.u * Mat(W) * svd.vt; //HZ 9.19 Mat_<double> t = svd.u.col(2); //u3 Matx34d P1( R(0,0),R(0,1), R(0,2), t(0), R(1,0),R(1,1), R(1,2), t(1), R(2,0),R(2,1), R(2,2), t(2));
很是简单。咱们须要作的工做仅是对咱们先前得到的本征矩阵进行奇异值分解(SVD),而且乘以一个特殊的矩阵W。对于咱们所作的数学上的操做不进行过深地阐述,咱们能够说SVD操做将咱们的矩阵E分解成为两部分,一个旋转成分和一个平移成分。实时上,本征矩阵起初是经过这两个成分相乘构成的。彻底地是知足咱们的好奇心,咱们能够看下面的本征矩阵的等式,它在字面意义上表现为:E=[t]xR(x是t的下标)。我能够看到它是由一个平移成分和一个旋转成分组成。
咱们注意到,咱们所作的仅是获得了一个相机矩阵,所以另一个相机矩阵呢?好的,咱们在假定一个相机矩阵是固定的而且是标准的(没有旋转和平移)状况下进行这个操做。下一个矩阵(这里的下一个矩阵表示相对于上面的P)也是标准的:
咱们从本征矩阵恢复的另一个相机相对于固定的相机进行了移动和旋转。这一样意味着咱们从这两个相机矩阵中恢复的任何三维点都是拥有第一个相机在世界坐标系中的原点(0,0,0)。
然而,这不是一个彻底解。H和Z在他们的书中展现了如何和为何这样的分解有四个可能的相机矩阵,可是仅有一个是正确的。正确的相机矩阵将产生带有一个正Z值(点在摄像机的前面)的重构点。可是咱们仅有当学了下一步将要讨论的三角测量和3维重建以后才能理解。
还有一个咱们能够考虑添加到咱们的方法中的事情是错误检测。屡次从点的匹配中计算基础矩阵是错误的,并影响相机矩阵。带有错误相机矩阵的进行三角测量是毫无心义的。咱们能够安装一个检查来检测是否旋转成分是一个有效的旋转矩阵。牢记旋转矩阵必须是一个行列式值为1(或者-1),咱们能够简单地进行以下作法:
bool CheckCoherentRotation(cv::Mat_<double>& R) { if(fabsf(determinant(R))-1.0 > 1e-07) { cerr<<"det(R) != +-1.0, this is not a rotation matrix"<<endl; return false; } return true; } We can now see how all these elements combine into a function that recovers the P matrices, as follows: void FindCameraMatrices(const Mat& K, const Mat& Kinv, const vector<KeyPoint>& imgpts1, const vector<KeyPoint>& imgpts2, Matx34d& P, Matx34d& P1, vector<DMatch>& matches, vector<CloudPoint>& outCloud ) { //Find camera matrices //Get Fundamental Matrix Mat F = GetFundamentalMat(imgpts1,imgpts2,matches); //Essential matrix: compute then extract cameras [R|t] Mat_<double> E = K.t() * F * K; //according to HZ (9.12) //decompose E to P' , HZ (9.19) SVD svd(E,SVD::MODIFY_A); Mat svd_u = svd.u; Mat svd_vt = svd.vt; Mat svd_w = svd.w; Matx33d W(0,-1,0,//HZ 9.13 1,0,0, 0,0,1); Mat_<double> R = svd_u * Mat(W) * svd_vt; //HZ 9.19 Mat_<double> t = svd_u.col(2); //u3 if (!CheckCoherentRotation(R)) { cout<<"resulting rotation is not coherent\n"; P1 = 0; return; } P1 = Matx34d(R(0,0),R(0,1),R(0,2),t(0), R(1,0),R(1,1),R(1,2),t(1), R(2,0),R(2,1),R(2,2),t(2)); }
此时,咱们拥有两个咱们须要用来重建场景的相机。第一个相机是标准的,存储在P变量中,第二个相机是咱们计算获得的,构成基础矩阵,存储在P1变量中。下一部分咱们将揭示如何使用这些相机来得到场景的三维重建。
接下来咱们看一下从咱们已经得到的信息中恢复场景的3D结构的事情。像先前所作的,咱们应当看一下,用来完成这个事情咱们手边所拥有的工具和信息。在前一部分,咱们从本征矩阵和矩阵矩阵中得到了两个相机矩阵。咱们已经讨论了这些工具用来得到空间中一点的3D位置是如何的有用。那么,我能够返回咱们的匹配点来用数据填充等式。点对一样对于计算咱们得到的近似计算的偏差有用。
如今咱们看一些如何使用OpenCV执行三角测量。此次咱们将会按照Tartley和Strum的三角测量的文章的步骤,文章中他们实现和比较了一些三角剖分的方法。咱们将实现他们的线性方法的一种,由于使用OpenCv很是容易编程。
回忆一下,咱们有两个由2D点匹配和P矩阵产生的关键等式:x=PX和x’=P’X,x和x’是匹配的二维点,X是两个相机进行拍照的真实世界三维点。若是咱们重写这个等式,咱们能够公式化为一个线性方程系统,该系统能够解出X的值,X正是咱们所指望寻找的值。假定X=(x,y,z,1)t(一个合理的点的假设,这些点离相机的中心不太近或者不太远)产生一个形式为AX=B的非齐次线性方程系统。咱们能够编码和解决这个问题,以下:
Mat_<double> LinearLSTriangulation( Point3d u,//homogenous image point (u,v,1) Matx34d P,//camera 1 matrix Point3d u1,//homogenous image point in 2nd camera Matx34d P1//camera 2 matrix ) { //build A matrix Matx43d A(u.x*P(2,0)-P(0,0),u.x*P(2,1)-P(0,1),u.x*P(2,2)-P(0,2), u.y*P(2,0)-P(1,0),u.y*P(2,1)-P(1,1),u.y*P(2,2)-P(1,2), u1.x*P1(2,0)-P1(0,0), u1.x*P1(2,1)-P1(0,1),u1.x*P1(2,2)-P1(0,2), u1.y*P1(2,0)-P1(1,0), u1.y*P1(2,1)-P1(1,1),u1.y*P1(2,2)-P1(1,2) ); //build B vector Matx41d B(-(u.x*P(2,3)-P(0,3)), -(u.y*P(2,3)-P(1,3)), -(u1.x*P1(2,3)-P1(0,3)), -(u1.y*P1(2,3)-P1(1,3))); //solve for X Mat_<double> X; solve(A,B,X,DECOMP_SVD); return X; }
咱们将得到由两个2维点产生的3D点的一个近似。还有一个要注意的事情是,2D点用齐次坐标表示,意味着x和y的值后面追加一个1。咱们必须确保这些点在标准化的坐标系中,这意味着它们必须乘以先前的标定矩阵K.咱们可能注意到咱们简单地利用KP矩阵来替代k矩阵和每个点相乘(每一点乘以KP),就是H和Z遍布第9章的作法那样。咱们如今能够写一个关于点匹配的循环语句来得到一个完整的三角测量,以下:
double TriangulatePoints( const vector<KeyPoint>& pt_set1, const vector<KeyPoint>& pt_set2, const Mat&Kinv, const Matx34d& P, const Matx34d& P1, vector<Point3d>& pointcloud) { vector<double> reproj_error; for (unsigned int i=0; i<pts_size; i++) { //convert to normalized homogeneous coordinates Point2f kp = pt_set1[i].pt; Point3d u(kp.x,kp.y,1.0); Mat_<double> um = Kinv * Mat_<double>(u); u = um.at<Point3d>(0); Point2f kp1 = pt_set2[i].pt; Point3d u1(kp1.x,kp1.y,1.0); Mat_<double> um1 = Kinv * Mat_<double>(u1); u1 = um1.at<Point3d>(0); //triangulate Mat_<double> X = LinearLSTriangulation(u,P,u1,P1); //calculate reprojection error Mat_<double> xPt_img = K * Mat(P1) * X; Point2f xPt_img_(xPt_img(0)/xPt_img(2),xPt_img(1)/xPt_img(2)); reproj_error.push_back(norm(xPt_img_-kp1)); //store 3D point pointcloud.push_back(Point3d(X(0),X(1),X(2))); } //return mean reprojection error Scalar me = mean(reproj_error); return me[0]; }
在下面的图像中咱们将看到一个两个图像三角测量的结果。这两个图像来至于P-11序列:
http://cvlab.epfl.ch/~strecha/multiview/denseMVS.html.
上面的两个图像是原始场景的两个视图,下面的一对图像是从两个视图重构获得的点云视图,包含着估计相机朝向喷泉。咱们能够看到右手边红色砖块墙部分是如何重构的,而且也能够看到突出于墙的喷泉。
然而,像咱们前面提到的那样,关于重构存在着这样的一个问题:重构仅能达到尺度上的。咱们应当花一些时间来理解达到尺度(up-to-scale)的意思。咱们得到的两个摄像机之间的运动存在一个随意测量的单元,也就是说,它不是用厘米或者英寸,而是简单地给出尺度单位。咱们重构的相机将是尺度单元距离上的分开。这在很大程度上暗示咱们应当决定过会从新得到更多的相机,由于每一对相机都拥有各自的尺度单元,而不是一个通常的尺度。
如今咱们将讨论咱们创建的偏差测量是如何可能的帮助咱们来找到一个更加鲁棒性的重构。首先,咱们须要注意重投影意味着咱们简单地利用三角化的3D点而且将这些点重塑到相机上以得到一个重投影的2D点。若是这个距离很大,这意味着咱们在三角测量中存在一个偏差,所以,在最后的结果中,咱们可能不想包含这个点。咱们的全局测量是平均重投影距离而且可能提供一个暗示——咱们的三角剖分整体执行的怎么样。高的重投影率可能代表P矩阵存在问题,所以该问题多是存在于本征矩阵的计算或者匹配特征点中。
咱们应当简短地回顾一下在前一部分咱们讨论的相机矩阵。咱们提到相机矩阵P1的分解能够经过四个不一样的方式进行分解,可是仅有一个分解是正确的。既然,咱们知道如何三角化一个点,咱们能够增长一个检测来看四个相机矩阵中哪个是有效地的。咱们应当跳过在一点实现上的细节,由于它们是本书随书实例代码中的精选(专题)。
下一步,咱们将要看一下从新得到直视场景的更多的相机,而且组合这些3D重建的结果。
既然咱们知道了如何从两个相机中恢复运动和场景几何,简单地应用相同的过程得到额外相机的参数和更多场景点彷佛不重要。事实上,这件事并不简单由于咱们仅能够得到一个达到尺度的重构,而且每对图像给咱们一个不一样的尺度。
有一些方法能够正确的从多个视场中重构3D场景数据。一个方法是后方交会(resection)或者相机姿态估计(camera pose estimation),也被称为PNP(Perspective-N-Point),这里咱们将使用咱们已经找到的场景点来解决一个新相机的位置。另一个方法是三角化更多的点而且看它们是如何适应于咱们存在的场景几何的。凭借ICP(Iterative Closest Point )咱们能够得到新相机的位置。在这一章,咱们将讨论使用OpenCV的sovlePnP函数来完成第一个方法。
第一步咱们选择这样的重构类型即便用相机后方交会的增长的3D重构,来得到一个基线场景结构。由于咱们将基于一个已知的场景结构来寻找任何相机的位置,咱们须要找到要处理的一个初始化的结构和一条基线。咱们可使用先前讨论的方法——例如,在第一个视频帧和第二个视频帧之间,经过寻找相机矩阵(使用FindCameraMatrices函数)来得到一条基线而且三角化几何(使用TriangulatePoints函数)。
发现一个基础结构后,咱们能够继续。然而,咱们的方法须要至关多的数据记录。首先,咱们须要注意solvePnP函数须要两个对齐的3D和2D点的矢量。对齐的矢量意味着一个矢量的第i个位置与另一个矢量的第i位置对齐。为了得到这些矢量,咱们须要在咱们早前恢复的3D点中找到这些点,这些点与在咱们新视频帧下的2D点是对齐的。完成这个的一个简单的方式是,对于云中的每个3D点,附加一个来至2D点的矢量。而后咱们可使用特征匹配来得到一个匹配对。
让咱们为一个3D点引入一个新的结构,以下:
struct CloudPoint { cv::Point3d pt; std::vector<int>index_of_2d_origin; };
它容纳,3D点和一个容器,容器内的元素为每帧图像上2D点的索引值,这些2D点用来计算3D点。当三角化一个新的3D点时,index_of_2d_origin的信息必须被初始化,来记录在三角化中哪些相机涉及到。然而,咱们可使用它来从咱们的3D点云追溯到每一帧上的2D点,以下:
std::vector<CloudPoint> pcloud; //our global 3D point cloud //check for matches between i'th frame and 0'th frame (and thus the current cloud) std::vector<cv::Point3f> ppcloud; std::vector<cv::Point2f> imgPoints; vector<int> pcloud_status(pcloud.size(),0); //scan the views we already used (good_views) for (set<int>::iterator done_view = good_views.begin(); done_view != good_views.end(); ++done_view) { int old_view = *done_view; //a view we already used for reconstrcution //check for matches_from_old_to_working between <working_view>'th frame and <old_view>'th frame (and thus the current cloud) std::vector<cv::DMatch> matches_from_old_to_working = matches_matrix[std::make_pair(old_view,working_view)]; //scan the 2D-2D matched-points for (unsigned int match_from_old_view=0; match_from_old_view<matches_from_old_to_working.size(); match_from_old_view++) { // the index of the matching 2D point in <old_view> int idx_in_old_view = matches_from_old_to_working[match_from_old_view].queryIdx; //scan the existing cloud to see if this point from <old_view> exists for (unsigned int pcldp=0; pcldp<pcloud.size(); pcldp++) { // see if this 2D point from <old_view> contributed to this 3D point in the cloud if (idx_in_old_view == pcloud[pcldp].index_of_2d_origin[old_view] && pcloud_status[pcldp] == 0) //prevent duplicates { //3d point in cloud ppcloud.push_back(pcloud[pcldp].pt); //2d point in image <working_view> Point2d pt_ = imgpts[working_view][matches_from_old_to_ working[match_from_old_view].trainIdx].pt; imgPoints.push_back(pt_); pcloud_status[pcldp] = 1; break; } } } } cout<<"found "<<ppcloud.size() <<" 3d-2d point correspondences"<<endl;
如今,咱们有一个场景中3D点到一个新视频帧中2D点水平对齐对,咱们可使用他们从新获得相机的位置,以下:
cv::Mat_<double> t,rvec,R; cv::solvePnPRansac(ppcloud, imgPoints, K, distcoeff, rvec, t, false); //get rotation in 3x3 matrix form Rodrigues(rvec, R); P1 = cv::Matx34d(R(0,0),R(0,1),R(0,2),t(0), R(1,0),R(1,1),R(1,2),t(1), R(2,0),R(2,1),R(2,2),t(2));
既然咱们正在使用sovlePnPRansac函数而不是sovlePnP函数,由于它对于异常值有更好的鲁棒性。既然咱们得到了一个新的P1矩阵,咱们能够简单的再次使用咱们早先定义的TriangualtePoints函数而且用更多的3D点来填充咱们的3D点云。
在下面的图像中,咱们看到一个增长的喷泉场景的重构(访问:http://cvlab.epfl.ch/~strecha/multiview/denseMVS.html,),从第四个图像。左上角的图像是使用了4个图像的重构;参加拍摄的相机用带有白线的红色简单表示,箭头表示了方向。其余的图像展现了更多 的相机来添加更多的点到点云中。
SfM方法中最重要的一个部分是提纯和最优化重构场景,一般被称做BA过程(Bundle Adjustment)。这是一个优化步骤,这里咱们得到的全部数据适应于一个统一的模型。3D点的位置和相机的位置都获得了最优化,所以重投影偏差最小(也就是说估计的3D点重投影到图像上接近于起源的2D点的位置)。这个过程一般须要解决带有几十千个参数指令的巨大的线性方程。这个过程可能会有些费力,可是咱们前面采起的步骤容许带有一个束调节(bundle adjuster)的简单的整合。前面看起来奇怪的事情变的清晰了。例如,咱们保留为点云中的每个3D点保存原始的2D点的理由。
束调节的一个实现算法是简单稀疏束条调节SSBA(Simple Sparse Bundle Adjustment)库。咱们将选择它做为咱们的BA优化器, 由于它拥有简单的API。它仅须要少许的输入参数,这些输入参数咱们能够至关简单的从咱们的数据结构中建立。SSBA中咱们使用的关键对象是CommonInternasMetricBundleOptimizer函数,这个函数执行最优化。它须要相机参数,3D点云,点云中每个点相对应的2D图像点,以及直视场景的相机。到如今为止,利用这些参数应该很直接。咱们应当注意这个BA方法假定全部图像经过一样的硬件获取,所以共同的内部,其余操做模式可能不须要假定这样的状况。咱们能够执行束调节(Bundle Adjustment),以下:
voidBundleAdjuster::adjustBundle( vector<CloudPoint>&pointcloud, const Mat&cam_intrinsics, conststd::vector<std::vector<cv::KeyPoint>>&imgpts, std::map<int ,cv::Matx34d>&Pmats ) { int N = Pmats.size(), M = pointcloud.size(), K = -1; cout<<"N (cams) = "<< N <<" M (points) = "<< M <<" K (measurements) = "<< K <<endl; StdDistortionFunction distortion; // intrinsic parameters matrix Matrix3x3d KMat; makeIdentityMatrix(KMat); KMat[0][0] = cam_intrinsics.at<double>(0,0); KMat[0][1] = cam_intrinsics.at<double>(0,1); KMat[0][2] = cam_intrinsics.at<double>(0,2); KMat[1][1] = cam_intrinsics.at<double>(1,1); KMat[1][2] = cam_intrinsics.at<double>(1,2); ... // 3D point cloud vector<Vector3d >Xs(M); for (int j = 0; j < M; ++j) { Xs[j][0] = pointcloud[j].pt.x; Xs[j][1] = pointcloud[j].pt.y; Xs[j][2] = pointcloud[j].pt.z; } cout<<"Read the 3D points."<<endl; // convert cameras to BA datastructs vector<CameraMatrix> cams(N); for (inti = 0; i< N; ++i) { intcamId = i; Matrix3x3d R; Vector3d T; Matx34d& P = Pmats[i]; R[0][0] = P(0,0); R[0][1] = P(0,1); R[0][2] = P(0,2); T[0] = P(0,3); R[1][0] = P(1,0); R[1][1] = P(1,1); R[1][2] = P(1,2); T[1] = P(1,3); R[2][0] = P(2,0); R[2][1] = P(2,1); R[2][2] = P(2,2); T[2] = P(2,3); cams[i].setIntrinsic(Knorm); cams[i].setRotation(R); cams[i].setTranslation(T); } cout<<"Read the cameras."<<endl; vector<Vector2d > measurements; vector<int> correspondingView; vector<int> correspondingPoint; // 2D corresponding points for (unsigned int k = 0; k <pointcloud.size(); ++k) { for (unsigned int i=0; i<pointcloud[k].imgpt_for_img.size(); i++) { if (pointcloud[k].imgpt_for_img[i] >= 0) { int view = i, point = k; Vector3d p, np; Point cvp = imgpts[i][pointcloud[k].imgpt_for_img[i]].pt; p[0] = cvp.x; p[1] = cvp.y; p[2] = 1.0; // Normalize the measurements to match the unit focal length. scaleVectorIP(1.0/f0, p); measurements.push_back(Vector2d(p[0], p[1])); correspondingView.push_back(view); correspondingPoint.push_back(point); } } } // end for (k) K = measurements.size(); cout<<"Read "<< K <<" valid 2D measurements."<<endl; ... // perform the bundle adjustment { CommonInternalsMetricBundleOptimizeropt(V3D::FULL_BUNDLE_FOCAL_ LENGTH_PP, inlierThreshold, K0, distortion, cams, Xs, measurements, correspondingView, correspondingPoint); opt.tau = 1e-3; opt.maxIterations = 50; opt.minimize(); cout<<"optimizer status = "<<opt.status<<endl; } ... //extract 3D points for (unsigned int j = 0; j <Xs.size(); ++j) { pointcloud[j].pt.x = Xs[j][0]; pointcloud[j].pt.y = Xs[j][1]; pointcloud[j].pt.z = Xs[j][2]; } //extract adjusted cameras for (int i = 0; i< N; ++i) { Matrix3x3d R = cams[i].getRotation(); Vector3d T = cams[i].getTranslation(); Matx34d P; P(0,0) = R[0][0]; P(0,1) = R[0][1]; P(0,2) = R[0][2]; P(0,3) = T[0]; P(1,0) = R[1][0]; P(1,1) = R[1][1]; P(1,2) = R[1][2]; P(1,3) = T[1]; P(2,0) = R[2][0]; P(2,1) = R[2][1]; P(2,2) = R[2][2]; P(2,3) = T[2]; Pmats[i] = P; } }
这个代码,虽然很长,是主要的关于转换咱们的内部数据结构到和来至SSBA的数据结构,而且调用最优化的过程。
下面的图像展现了BA的效果。左边的两个图像调整前的点云中的点,来至两个视角的观察,而且右边的图像展现了优化后的点云。变化至关明显,而且从不一样视图获得三角化点之间的不重合如今大部分统一了。咱们一样能够注意到调整建立了一个更好的平面重建。
当操做3D数据时,经过简单地观察重投影偏差测量或原始点信息很难快速的理解结果是否正确。另外一方面,若是咱们观察点云(itself),咱们能够当即的检查这个点是否有意义或者存在偏差。为了可视化,咱们将使用一个颇有前途的OpenCV的姊妹工程,称为点云库(Point Cloud Librar)(PCL)。它带有许多可视化和分析点云的工具,例如找到一个平面,匹配点云,分割目标以及排除异常值。若是咱们的目标不是一个点云,而是一些高级信息例如3D模型,这些工具将很是有用。
首先,咱们须要在PCL的数据结构中表示咱们的点云(本质上是3D点的列表)。能够经过以下的作法实现:
pcl::PointCloud<pcl::PointXYZRGB>::Ptr cloud; void PopulatePCLPointCloud(const vector<Point3d>& pointcloud, const std::vector<cv::Vec3b>& pointcloud_RGB ) //Populate point cloud { cout<<"Creating point cloud..."; cloud.reset(new pcl::PointCloud<pcl::PointXYZRGB>); for (unsigned int i=0; i<pointcloud.size(); i++) { // get the RGB color value for the point Vec3b rgbv(255,255,255); if (pointcloud_RGB.size() >= i) { rgbv = pointcloud_RGB[i]; } // check for erroneous coordinates (NaN, Inf, etc.) if (pointcloud[i].x != pointcloud[i].x || isnan(pointcloud[i].x) || pointcloud[i].y != pointcloud[i].y || isnan(pointcloud[i].y) || pointcloud[i].z != pointcloud[i].z || isnan(pointcloud[i].z) || fabsf(pointcloud[i].x) > 10.0 || fabsf(pointcloud[i].y) > 10.0 || fabsf(pointcloud[i].z) > 10.0) { continue; } pcl::PointXYZRGB pclp; // 3D coordinates pclp.x = pointcloud[i].x; pclp.y = pointcloud[i].y; pclp.z = pointcloud[i].z; // RGB color, needs to be represented as an integer uint32_t rgb = ((uint32_t)rgbv[2] << 16 | (uint32_t)rgbv[1] << 8 | (uint32_t)rgbv[0]); pclp.rgb = *reinterpret_cast<float*>(&rgb); cloud->push_back(pclp); } cloud->width = (uint32_t) cloud->points.size(); // number of points cloud->height = 1; // a list of points, one row of data }
为了可视化有一个好的效果,咱们也能够提供彩色数据如同图像中的RGB值。咱们一样也能够对原始点云应用一个滤波器,这将消除多是异常值的点,使用统计移除移除(statistical outlier removal)(SOR)工具以下:
Void SORFilter() { pcl::PointCloud<pcl::PointXYZRGB>::Ptr cloud_filtered (new pcl::PointC loud<pcl::PointXYZRGB>); std::cerr<<"Cloud before SOR filtering: "<< cloud->width * cloud->height <<" data points"<<std::endl; // Create the filtering object pcl::StatisticalOutlierRemoval<pcl::PointXYZRGB>sor; sor.setInputCloud (cloud); sor.setMeanK (50); sor.setStddevMulThresh (1.0); sor.filter (*cloud_filtered); std::cerr<<"Cloud after SOR filtering: "<<cloud_filtered->width * cloud_filtered->height <<" data points "<<std::endl; copyPointCloud(*cloud_filtered,*cloud); }
而后,咱们可使用PCL的API来运行一个简单的点云的可视化器,以下:
Void RunVisualization(const vector<cv::Point3d>& pointcloud, const std::vector<cv::Vec3b>& pointcloud_RGB) { PopulatePCLPointCloud(pointcloud,pointcloud_RGB); SORFilter(); copyPointCloud(*cloud,*orig_cloud); pcl::visualization::CloudViewer viewer("Cloud Viewer"); // run the cloud viewer viewer.showCloud(orig_cloud,"orig"); while (!viewer.wasStopped ()) { // NOP } }
下面的图像展现了统计移除移除工具(statistical outlier removal tool)使用以后的输出结果。左手边的图像是SfM的原始结果点云。右手边的图像展现通过SOR操做滤波以后的点云。咱们可以注意到一些离群的点被移除了,剩下了一个更干净的点云。
咱们能够在这本书提供的材料中找到SfM的实例代码。咱们如今看一些怎么样编译,运行和利用它。代码使用CMake,一个交叉编译环境,相似于Maven或者SCons。咱们一样应当确保咱们有下面的全部前提条件来编译咱们的应用程序:
• OpenCV v2.3 or highe
• PCL v1.6 or higher
• SSBA v3.0 or higher
首先,咱们必须创建编译环境。为此,咱们可能建立一个文件夹,命名为build,咱们将全部编译相关的文件存储在这里。如今咱们将假定全部的命令行操做都是在build/文件夹内,虽然这个过程是相似的(取决于文件的位置),即便没有使用build文件夹。
咱们应当确保CMake能够找到SSBA和PCL,若是PCL正确的安装了,那没有问题。然而,咱们必须经过-DSSBA_LIBRARY_DIR=...编译参数设置正确的位置来找到SSBA的预编译库。若是咱们正在使用Windows操做系统,咱们可使用Microsoft Visual Studio来编译,所以,咱们应当运行下面的的命令:
cmake –G "Visual Studio 10" -DSSBA_LIBRARY_DIR=../3rdparty/SSBA-3.0/ build/ ..
若是咱们使用Linux,Mac Os,或者其余Unix-Like操做系统,咱们执行下面的命令:
cmake –G "Unix Makefiles" -DSSBA_LIBRARY_DIR=../3rdparty/SSBA-3.0/build/ ..
若是咱们喜欢使用MacOS上的XCode,执行下面的命令:
cmake –G Xcode -DSSBA_LIBRARY_DIR=../3rdparty/SSBA-3.0/build/ ..
CMake一样能够为Eclipse,Codeblocks,和更多的环境编译宏命令。CMake完成建立编译环境以后,咱们准备编译。若是咱们正在使用一个Unix-like系统,咱们能够简单地执行这个生成工具(the make utility),不然咱们应当使用咱们开发环境的编译过程。
编译完成以后,咱们应当得到了一个执行程序名为ExploringSfMExex,这用来运行SfM过程。不带参数运行这个程序会致使以下的显示: USAGE: ./ExploringSfMExec <path_to_images>
为了在图像集上执行这个过程,咱们应当提供位置做为驱动来找到图像文件。若是提供了有效的位置,过程开始,我应当看到这个进程和屏幕上的调试信息。程序的结束将显示源于图像的点云。按1和2键,能够切换到调整的(adjusted)点云和非调整的(non-adjusted)点云。
在这一章,咱们已经看到了OpenCV是怎样用一个既简单编码又好理解的方式来帮助咱们处理来至运动的结构。OpenCV的API包含了一些有用的函数和数据结构,这使得咱们生活的更加轻松,一样地协助咱们有一个更清洁的实现。
然而,最早进的SfM方法更复杂。那里存在不少问题,咱们选择忽略,喜欢简单化,以及在这些地方一般有更多的错误检查。对于不一样的SfM成分,咱们选择的方法一样能够再次访问。例如,H和Z提出了一个高精度的三角测量方法,甚至使用N-视图三角测量,曾经他们利用多个图像理解特征之间的关系。
若是咱们想延伸和深化熟悉SfM,固然,咱们将从观察其余开源SfM库中收益。一个特别感兴趣的工程是libMV,它实现了大量的SfM成分,经过互换,这可能得到最好的结果。华盛顿大学有一个伟大的做品,为不少类型SfM(Bundler and Visual SfM)提供了工具。这项做品的灵感来至于微软的在线产品,称做PhotoSynth。网上有不少容易访问到的SfM的实现,而且一我的仅需来找到更多的SfM的实现。
咱们没有深刻讨论的另一个重要的关系是SfM和可视化定位以及映射,更好的称为同步定位和映射(SLAM)方法。在这一章,咱们已经处理给出的图像集和一个视频序列,而且在这些状况下,使用SfM是可行的。然而,一些应用没有预记录的数据集而且必须动态时引导重建。这个过程被成为映射,而且当咱们正在使用特征匹配和2D跟踪以及三角测量后来建立世界的3D映射时,这个过程被完成。
在下一章,咱们将看到如何使用机器学习中的各类技术利用OpenCV从图像中提出车牌数字。
• Multiple View Geometry in Computer Vision,Richard Hartley and Andrew Zisserman,Cambridge University Press• Triangulation, Richard I. Hartley and Peter Sturm, Computer vision and image understanding, Vol. 68, pp. 146-157• http://cvlab.epfl.ch/~strecha/multiview/denseMVS.html• On Benchmarking Camera Calibration and Multi-View Stereo for High Resolution Imagery,C. Strecha, W. von Hansen, L. Van Gool, P. Fua,and U. Thoennessen, CVPR• http://www.inf.ethz.ch/personal/chzach/opensource.html• http://www.ics.forth.gr/~lourakis/sba/• http://code.google.com/p/libmv/• http://www.cs.washington.edu/homes/ccwu/vsfm/• http://phototour.cs.washington.edu/bundler/• http://photosynth.net/• http://en.wikipedia.org/wiki/Simultaneous_localization_and_mapping• http://pointclouds.org• http://www.cmake.org