简介基于C与OpenCV实现的增量式三维重建算法工程资料覆盖计算机视觉与数字摄影测量课程设计内容适合作为毕业设计、课程设计或初期项目参考面向具有一定编程基础、希望深入理解SFM流程的学习者。资料共100个文件压缩包约68.9MB主体为2个cpp源码与1个h头文件另含1个md说明文档、66个txt配置或数据文件以及28张测试图像目录结构与代码组织清晰便于根据README快速上手。算法实现重点解决航迹恢复问题通过暴力匹配与最小生成树建立图像关联同时提供经过检校的相机内参和背景处理后的模型数据可基于OpenCV、Ceres与PCL环境复现完整的增量式重建流程。资源目前已有73人学习主要价值在于可运行的代码骨架、环境依赖说明和真实数据样例不过资源声明仅作为参考资料需要读者自行调试排错并扩展功能不适合直接照搬提交。1. 增量式三维重建的工程切面从航迹恢复说起拿一套带检校内参的相机拍摄的模型照片直接丢给colmap一键出点云并不难难的是把重建过程拆成一块块能讲清楚原理、能自己写的代码。这套基于C/OpenCV的增量式三维重建项目正好覆盖了这条完整链路特征匹配、航迹恢复、初始像对选择、增量式位姿估计、三角化、BA平差最后用PCL做点云后处理。项目文件里main.cpp负责重建主流程bundle_adjustment.h封装了ceres优化cloudpointprocess.cpp单独处理点云一上来就是工程结构而不是算法玩具。真正作业时最花时间的是航迹恢复——判断哪些图相邻、以什么顺序加入重建。这里有个反直觉的结论照片只有十来张时手工排序比最小生成树更可靠照片上百张时MST才是被验证过的方案。本文用这个项目的数据和代码把从匹配到BA的参数逐条过一遍cams_1里的检校内参怎么用、哪些阈值值得调都会说清楚。作为计算机视觉大作业或数字摄影测量课程设计这套代码的模块拆分可以直接复用。2. 特征匹配与初始外极几何OpenCV 4.5.5下的选型和参数特征匹配是整个重建质量的上限。后面的PNP和BA都是在对应关系上做数值优化匹配错了BA再强也拉不回来。项目README里标明的OpenCV版本是4.5.5这个版本的特征模块接口和3.x相比变化不大但细分模块之间有明显的性能差异值得先讲清楚。2.1 ORB还是SIFT匹配距离、专利与工程取舍项目没有把SIFT作为主特征原因很实际OpenCV 4.5.5主库里SIFT位于opencv_contrib的xfeatures2d模块需要自己编译完整版OpenCV课程设计和毕设环境下不必要的依赖能少则少。ORB是二进制描述子用Hamming距离做暴力匹配速度比SIFT浮点描述子的L2距离快一个数量级对旋转和尺度缩放也有基本不变性。对本项目使用的模型测试数据来说ORB提取2000个点的匹配数量足够支撑后续的PnP和BA。特征描述子类型匹配距离典型提取速度项目适用性ORB二进制Hamming快主库自带课程设计首选SIFT浮点128维L2慢需contrib精度高但编译链复杂SURF浮点64维L2中需contrib不推荐新项目基本弃用如果换用SIFT描述子维度变高BFMatcher的耗时和内存都会显著上升匹配比例测试的0.75阈值依然有效但计算距离的度量要从NORM_HAMMING改成NORM_L2。对增量式重建来说特征点数量比特征质量更影响稳定性ORB的2000个点正好是速度与覆盖率之间的平衡点。2.2 BFMatcher、比例测试与RANSAC剔除的代码落地读取两张图进行特征提取和匹配的代码可以直接复制到main.cpp的初始化部分我一般封装成一个matchTwoImages函数返回内外点索引和基础矩阵F。#include opencv2/features2d.hpp #include opencv2/calib3d.hpp cv::Ptrcv::Feature2D detector cv::ORB::create( 2000, // nfeatures单张图最多保留的特征点数 1.2f, // scaleFactor金字塔每层尺度缩小比例 8, // nlevels金字塔层数 31, // edgeThreshold靠近图像边界31像素内的点不提取 0 // firstLevel从第0层开始 ); std::vectorcv::KeyPoint kp1, kp2; cv::Mat desc1, desc2; detector-detectAndCompute(img1, cv::noArray(), kp1, desc1); detector-detectAndCompute(img2, cv::noArray(), kp2, desc2); cv::BFMatcher matcher(cv::NORM_HAMMING, false); std::vectorstd::vectorcv::DMatch knnMatches; matcher.knnMatch(desc1, desc2, knnMatches, 2); std::vectorcv::DMatch goodMatches; for (const auto m : knnMatches) { if (m.size() 2 m[0].distance 0.75f * m[1].distance) { goodMatches.push_back(m[0]); } }BFMatcher构造第一个参数NORM_HAMMING表示用汉明距离比较ORB的二进制描述子第二个参数crossCheck设为false因为后面通过knnMatch取每个特征点的最近邻和次近邻。比例测试的0.75来自Lowe的SIFT论文最近邻距离与次近邻距离之比小于0.75说明特征足够独特。这个阈值在ORB上可以放宽到0.8但放到0.9时误匹配会明显增多。项目数据做过背景处理物体边缘和纹理清晰0.75能留下约一半匹配后续RANSAC收敛很快。匹配完必须过外极几何约束否则误匹配会直接进入PnPstd::vectorcv::Point2f pts1, pts2; for (const auto m : goodMatches) { pts1.push_back(kp1[m.queryIdx].pt); pts2.push_back(kp2[m.trainIdx].pt); } cv::Mat mask; cv::Mat F cv::findFundamentalMat(pts1, pts2, cv::FM_RANSAC, 1.0, 0.99, mask); int inlierCount cv::countNonZero(mask); double inlierRatio static_castdouble(inlierCount) / goodMatches.size();findFundamentalMat的第三个参数FM_RANSAC表示用RANSAC估计基础矩阵1.0是Sampson距离阈值单位像素0.99是期望置信度。mask中非零元素对应的就是内点。如果inlierRatio小于40%这对图的匹配质量就很差直接跳过不进初始对候选池。注意pts1、pts2拷贝成Point2f数组后原始索引关系丢失我在实际代码中额外维护一个索引对容器因为后面增量式加帧时还需要queryIdx和trainIdx来回溯三维点与图像点的对应关系。3. 航迹恢复与初始像对最小生成树和本质矩阵的四种解航迹恢复是SFM里最容易被忽略但直接影响全局质量的一环。正文里提到“SFM最重要的问题就是航迹恢复”这句话很真实航迹排错了后面的重建就是在一个错误的帧序列上做局部优化误差不可逆。3.1 把航迹恢复建模成图优化问题把所有图片两两匹配匹配点数量作为边的权重目标是找到一个让所有图连通且权重之和最大的边集合这就是最大生成树问题。把权重取倒数转为代价后等价于最小生成树问题。代码里排序时我不用原始匹配数而用RANSAC后的内点数作为权重因为原始匹配数里可能混入纹理重复区域的误匹配。struct Edge { int i, j; int matchesInliers; int rawMatches; }; // 按内点数从大到小排序Kruskal最大生成树 std::sort(edges.begin(), edges.end(), [](const Edge a, const Edge b) { return a.matchesInliers b.matchesInliers; }); std::vectorint uf(n); std::iota(uf.begin(), uf.end(), 0); // 滑动窗口维护并查集连通性判断省略 for (const auto e : edges) { if (findRoot(e.i) ! findRoot(e.j)) { unionRoot(e.i, e.j); mst.push_back(e); } }实际工程中我用“内点占比×内点数”作为排序权重在质量和数量之间取平衡。照片只有十几张时肉眼能直接看出拍摄顺序代码里可以先跑两两匹配打印匹配矩阵人工确认后再指定初始链这比盲目相信MST更稳。真实航测任务几百张图时暴力匹配O(n^2)对图的开销太大常见做法是先做词袋检索候选对再对候选对做特征匹配最后用MST恢复航迹。项目自带的测试数据经过背景处理四组图片005、010、029、079之间的匹配数差异明显MST选出的初始对通常就是纹理最丰富的那一帧。3.2 本质矩阵分解与四种解的消歧义确定初始像对后用cams_1文件夹中检校过的相机内参K1、K2通过E K2.t() * F * K1算出本质矩阵再做SVD分解得到候选旋转和平移。这是外极几何的标准步骤数字摄影测量课程里对应的是相对定向元素的计算。cv::Mat E K2.t() * F * K1; cv::SVD svd(E, cv::SVD::FULL_UV); // 本质矩阵的奇异值理论上为(1,1,0)实际有噪声需要强制归一化 cv::Mat W (cv::Mat_double(3,3) 0,-1,0, 1,0,0, 0,0,1); cv::Mat R1 svd.u * W * svd.vt; cv::Mat R2 svd.u * W.t() * svd.vt; cv::Mat t1 svd.u.col(2); cv::Mat t2 -svd.u.col(2);需要注意噪声环境下svd.w并不严格等于(1,1,0)比较严谨的做法是用E U * diag(1,1,0) * V.t()对E做投影后再分解。R1、R2与t1、t2组合成四种相机姿态只有一种是真实解。判定标准是三角化后的三维点在两个相机坐标系下的深度都为正。实际代码中我对四种组合分别做一次三角化统计正深度点数。组合旋转平移判定原则AR1t1正深度点数最多通常为正确解BR1t2点在两相机后方或左右矛盾CR2t1与A共享平移旋转错配DR2t2镜像翻转直接排除注意正深度点数的统计必须左右相机同时满足z0不能只看左相机。只统计单侧深度的做法在相机朝向相近时会把错误解也选进去。如果最优解的左右正深度占比都小于80%说明这个初始像对本身不可靠我一般直接放弃换下一对匹配内点数最多的图重试。这一步选错后续所有增量式位姿估计都在错误的坐标框架下累积漂移而且不可逆。cams_1的内参在求F之前先用undistortPoints把像素坐标转成去畸变后的归一化坐标比先求F再单独处理畸变更稳定。4. 增量式位姿估计与三角化一张一张把地图长出来初始对重建出第一批三维点后系统进入增量式循环。每加入一张新图分成两个步骤先用solvePnPRansac求新帧位姿再对新增可见的匹配点做三角化。这个循环是bundle_adjustment.h前面主循环的核心代码结构不复杂但参数和顺序很讲究。4.1 solvePnPRansac利用已重建三维点约束新帧位姿当已有的三维点云中有一部分能在新帧中找到对应2D投影时新帧位姿就是一个标准的2D-3D配准问题。std::vectorcv::Point3f objPts; std::vectorcv::Point2f imgPts; cv::Mat rvec, tvec; bool ok cv::solvePnPRansac( objPts, imgPts, K, distCoeffs, rvec, tvec, true, // useExtrinsicGuess用上一帧的R/t作为初值 1000, // 迭代次数 8.0, // 重投影误差阈值单位像素 0.99, // 置信度 cv::SOLVEPNP_ITERATIVE);SOLVEPNP_ITERATIVE内部用Levenberg-Marquardt做非线性优化useExtrinsicGuess为true时把上一帧的姿态作为初值连续帧之间运动小收敛快。8.0像素的重投影误差阈值在这个测试数据上偏宽松原因是模型表面纹理重复度高少量误匹配进入RANSAC后足够高的迭代次数能保证内点率稳定。如果实际图集中特征更干净阈值收到3像素效果更好。这里有一个经常踩的坑objPts和imgPts的顺序必须一一对应生成二者时如果处理不当比如一个对应第一张图的匹配索引而另一个对应第二张图的索引整个RANSAC结果都是错的。我维护一个带三个字段的结构体matchIndex、point3DIndex、point2DIndex从匹配链建立对应关系不靠两个vector的隐式顺序。求解出rvec和tvec之后常见做法是把新帧的可见三维点都做一次重投影删掉残差过大的观测再用ceres对位姿做一次单独的精修最后才进入三角化阶段。不做这步精修而直接三角化新插入的三维点会带着位姿误差进入BA增加后续优化的负担。4.2 triangulatePoints的输入约束与三重过滤新帧与已有关键帧之间的新增匹配对通过两个已知位姿恢复三维坐标。cv::Mat proj1 K * cv::Mat(rt1); // rt1为3x4的[R|t] cv::Mat proj2 K * cv::Mat(rt2); cv::Mat points4D; cv::triangulatePoints(proj1, proj2, pts1, pts2, points4D); // 齐次坐标转非齐次 std::vectorcv::Point3f points3D; for (int i 0; i points4D.cols; i) { double w points4D.atdouble(3, i); points3D.push_back(cv::Point3f( points4D.atdouble(0, i) / w, points4D.atdouble(1, i) / w, points4D.atdouble(2, i) / w)); }triangulatePoints的输入是3x4投影矩阵输出是4xN的齐次坐标每个点的w分量做归一化后才能得到非齐次XYZ。最容易出错的是投影矩阵的构造顺序cv::Mat拼接R和t时要确认是[R|t]而不是[t|R]这个错误不会报错但会把所有三角化点投到错误位置。distCoeffs参与PnP但不参与triangulatePoints正确做法是先undistortPoints把像素坐标转成去畸变坐标系再传给三角化接口。三角化之后必须过滤我按三个条件依次执行重投影误差把新三维点投到两帧上与原始2D点距离超过2像素的删除。三角化角度三维点与两相机光心的夹角在2度到60度之间角度太小深度不确定角度太大匹配误差放大。深度符号三维点在两帧相机坐标系下的z值必须为正。三个条件分开过滤并且分别打印删除了多少点。如果某一轮角度过滤删除比例突然超过一半说明新帧与参考帧基线过短视差不足。这种情况下PnP能算出姿态但三角化出来的点质量很差。我一般用平均视差作为插入关键帧的条件新帧与最近关键帧的平均像素视差大于30像素才允许三角化否则只跟踪位姿不新增点。这样可以显著减少无效计算避免点云在同一区域堆叠过于密集。5. 用ceres做Bundle Adjustment重投影误差最小化的边界增量式位姿估计是局部最优化误差随帧数累积。Bundle Adjustment把所有相机位姿和三维点放在一起重新优化把累计误差重新分布到所有变量上。项目里bundle_adjustment.h封装的是ceres 2.0.0与OpenCV的solvePnP不同ceres没有内置PnP但它给BA提供了更自由的代价函数定义。5.1 9维相机参数残差块的写法BA的相机参数一般有两种表达方式一种是旋转矩阵加平移向量另一种是旋转向量加平移向量。ceres里用旋转向量配合AngleAxisRotatePoint更稳妥因为它避免了欧拉角的万向锁同时让雅可比计算保持连续。struct ReprojectionError { ReprojectionError(double fx_, double fy_, double cx_, double cy_, const Eigen::Vector2d obs) : fx_(fx_), fy_(fy_), cx_(cx_), cy_(cy_), obs_(obs) {} template typename T bool operator()(const T* const camera, const T* const point, T* residuals) const { // camera: [rx, ry, rz, tx, ty, tz, fx, fy, cx, cy] T p[3]; ceres::AngleAxisRotatePoint(camera, point, p); p[0] camera[3]; p[1] camera[4]; p[2] camera[5]; T xp p[0] / p[2]; T yp p[1] / p[2]; residuals[0] camera[6] * xp camera[8] - obs_(0); residuals[1] camera[7] * yp camera[9] - obs_(1); return true; } private: double fx_, fy_, cx_, cy_; Eigen::Vector2d obs_; };这个残差块里camera数组的前三维是旋转向量中间三维是平移最后四维是fx、fy、cx、cy。ceres的AutoDiffCostFunction通过operator()模板自动计算雅可比残差值是归一化坐标投影后与观测值的差值。ceres 2.0.0对自动求导做了不少优化9维相机参数的雅可比矩阵计算开销比手写解析雅可比低不少写起来也省事。这里要强调一个边界问题如果fx、fy、cx、cy全部放开优化它们会和相机平移产生尺度耦合。单目序列本身没有绝对尺度内参如果参与优化优化器可能把fx调到极小来吸收深度误差导致重建的点云整体变形。常见做法是固定cx、cy只优化fx、fy或者干脆全部固定只优化位姿和三维点。5.2 损失函数与内参是否放开优化cams_1文件夹中提供了经过检校的内参这说明相机内参是可信的。实际代码里我倾向于把fx、fy也设成常量只保留位姿和三维点作为优化变量。这样BA的变量个数减少迭代收敛更快而且避免了内参与平移的耦合。如果你的图集没有标定文件内参未知才需要把fx、fy加入优化但要把cx、cy固定为图像中心。损失函数推荐HuberLoss而不是TrivialLossTrivialLoss等于最小二乘没有鲁棒性几个误匹配就能把优化结果拉偏。HuberLoss(1.0)在残差小于1像素时按二次惩罚大于1像素时按线性惩罚既能平滑中等噪声又不会过度压制真实的大残差。CauchyLoss(0.5)在极端离群值上更平滑但它会把本属于正常范围的重投影误差也压平导致优化提前收敛最终的RMSE反而偏高。BA的触发时机比损失函数的选择更重要。每加入一到两帧新关键帧做一次local BA只优化最近N个关键帧以及它们共同观测到的三维点。全局BA放到所有帧加完后连续跑三轮每轮迭代不少于50次。ceres的options里num_threads4linear_solver_type选SPARSE_SCHUR因为BA问题的信息矩阵是稀疏块结构SCHUR消元利用了这一特性。6. PCL点云后处理与重建质量的三个验证指标cloudpointprocess.cpp是独立于重建主流程的处理模块输入是BA后的稀疏三维点输出是经过滤波和降采样的PLY点云。这个模块的存在说明项目作者把重建和点云后处理解耦了替换算法时不需要动主流程。6.1 统计滤波与体素滤波的调参BA输出的稀疏点云带有离群点和重复观测点直接可视化会很乱。统计滤波按每个点与近邻的平均距离剔除离群点体素滤波再做均匀降采样。pcl::StatisticalOutlierRemovalpcl::PointXYZ sor; sor.setInputCloud(cloud); sor.setMeanK(30); sor.setStddevMulThresh(1.0); pcl::PointCloudpcl::PointXYZ::Ptr cloudFiltered(new pcl::PointCloudpcl::PointXYZ); sor.filter(*cloudFiltered); pcl::VoxelGridpcl::PointXYZ voxel; voxel.setInputCloud(cloudFiltered); voxel.setLeafSize(0.005f, 0.005f, 0.005f); pcl::PointCloudpcl::PointXYZ::Ptr cloudDownsampled(new pcl::PointCloudpcl::PointXYZ); voxel.filter(*cloudDownsampled);setMeanK30表示统计每个点最近30个邻居的平均距离与全局平均距离比较超过1.0倍标准差就删除。稀疏点云中每个点的邻居本来就少MeanK可以降到10到20否则邻近统计会把合法的边界点也当成离群点。体素滤波的leaf size取决于被重建物体的物理尺度。如果重建的是建筑或场景0.005m会让点云被抽稀得面目全非至少要放大到0.05m以上。6.2 验证重建质量的三个维度重建质量不能只看可视化效果要有量化指标。第一个是重投影误差RMSE。BA完成后对所有观测计算残差均方根正常应该在0.5到1.5像素之间。超过2像素说明有大量误匹配进入了优化或相机内参与实际不符。第二个是三维点的track长度分布。每个三维点被多少个相机观测到被4帧以上观测的点占比高说明三角化稳定如果大量点只有单帧观测说明增量式循环里有断链通常是新帧与参考帧基线太短或匹配被过滤过严。第三个是绝对尺度验证。单目SFM无法恢复真实尺度本质矩阵分解出的t是归一化后的用模型上已知长度的线段比如标定板格子边长或物体上的刻度对比重建结果计算一个全局scale系数乘到点云上点云才有测量的意义。最后一个实用技巧从cams_1直接拷入内参后注意cv::Mat的类型必须是CV_64FsolvePnPRansac和triangulatePoints对矩阵类型敏感类型不对会报type error。保存点云时用PLY二进制格式检查重建结果用CloudCompare打开比pcl_viewer的交互方便很多。本文还有配套的精品资源点击获取