Python实战：单目三维重建从原理到实现

1. 单目三维重建用一张照片“猜”出三维世界你有没有想过为什么我们看一张普通的照片就能大概知道里面物体的远近和形状比如一张桌子的照片你一眼就能看出哪个杯子在前面哪个花瓶在后面。我们的大脑天生就擅长从二维图像中“脑补”出三维信息。单目三维重建就是让计算机学会这个“脑补”技能。简单来说单目三维重建就是只用一个普通摄像头比如你的手机摄像头通过拍摄物体在不同位置的多张照片来恢复出这个物体在真实空间中的三维形状和位置。这听起来有点像魔法但背后的原理其实很直观想象一下你闭上一只眼睛只用一只眼睛看世界然后稍微移动一下头部你看到的物体位置会发生变化。你的大脑就是通过分析这些微小的位置变化来判断物体的深度和距离的。计算机做单目三维重建本质上就是在模仿这个过程。这个过程非常适合有一定Python基础并且对计算机视觉感兴趣的开发者来学习。你不需要昂贵的深度摄像头或激光雷达只需要一个普通的摄像头和你的代码就能开启三维世界的大门。学完这个你不仅能理解手机“人像模式”虚化背景的原理还能为自己的机器人项目添加“视觉感知”甚至尝试做一些简单的三维扫描。我自己在刚开始接触这个领域时就被这种“无中生有”的能力深深吸引了虽然踩了不少坑但一步步实现出来的成就感是无与伦比的。2. 环境搭建与相机标定给你的摄像头做“体检”万事开头难但第一步往往是最关键的。在开始写代码“重建”世界之前我们得先把工具准备好并且最重要的是了解我们手中“武器”——摄像头的特性。2.1 搭建你的Python视觉开发环境我强烈建议你使用Anaconda来管理你的Python环境它能很好地解决不同项目之间库版本冲突的噩梦。创建一个新的环境比如叫3d_reconstruction然后安装我们需要的核心库。conda create -n 3d_reconstruction python3.8 conda activate 3d_reconstruction pip install opencv-python3.4.2.16 opencv-contrib-python3.4.2.16 pip install numpy matplotlib pillow这里有个小坑我踩过OpenCV的版本要特别注意。我们用的3.4.2.16是一个比较经典的版本它包含了我们后面会用到的SURF和SIFT这些专利算法在新版本中这些算法被移到了opencv-contrib中且可能需要单独编译。直接用这个版本可以省去很多麻烦。如果你安装的是更新的OpenCV 4.x代码中的cv2.xfeatures2d_SURF等模块可能需要不同的调用方式对新手不友好。2.2 理解并完成相机标定相机标定听起来很高大上其实就像是给摄像头做一次全面的“体检”。我们得知道它的“视力”怎么样——有没有近视焦距是多少、有没有散光畸变参数。任何摄像头拍出的照片都不是完美的透视投影镜头本身会引入畸变比如直线拍出来变弯了桶形畸变或枕形畸变。标定的核心思想是用一个我们已知精确三维尺寸的物体通常是棋盘格作为“标尺”让摄像头从不同角度拍摄它。通过检测图片中棋盘格角点的二维像素位置和我们已知的角点三维世界坐标建立对应关系从而反推出摄像头的内部参数焦距、光心和畸变系数。我们来详细拆解一下代码里的camera_calibration函数。首先你需要打印一张棋盘格图比如11x8的内角点贴在一个平整的板子上然后用你的摄像头从不同角度、不同距离拍摄至少15-20张照片保存在Checkerboard_Image文件夹里。照片要拍得清晰棋盘格要完整出现在画面中并且角度要多样。def camera_calibration(ImagePath): # 设置角点精细化迭代的终止条件最多迭代30次或精度达到0.001 criteria (cv2.TERM_CRITERIA_EPS cv2.TERM_CRITERIA_MAX_ITER, 30, 0.001) # 棋盘格尺寸注意这里指的是内部角点的数量比如11行8列 row 11 column 8 # 准备世界坐标系中的角点坐标假设棋盘格在Z0的平面上每个格子边长假设为1个单位 objpoint np.zeros((row * column, 3), np.float32) objpoint[:, :2] np.mgrid[0:row, 0:column].T.reshape(-1, 2) objpoints [] # 存储所有图片的世界坐标 imgpoints [] # 存储所有图片的像素坐标 batch_images glob.glob(ImagePath /*.jpg) for i, fname in enumerate(batch_images): img cv2.imread(batch_images[i]) imgGray cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) # 在灰度图中寻找棋盘格角点 ret, corners cv2.findChessboardCorners(imgGray, (row, column), None) if ret: objpoints.append(objpoint) # 将角点位置精确到亚像素级别提高标定精度 corners2 cv2.cornerSubPix(imgGray, corners, (11, 11), (-1, -1), criteria) imgpoints.append(corners2) # 可视化角点并保存方便检查 img cv2.drawChessboardCorners(img, (row, column), corners2, ret) cv2.imwrite(Checkerboard_Image/Temp_JPG/Temp_ str(i) .jpg, img) print(f成功提取了 {len(batch_images)} 张图片的角点) # 核心标定函数 ret, mtx, dist, rvecs, tvecs cv2.calibrateCamera(objpoints, imgpoints, imgGray.shape[::-1], None, None) return mtx, dist运行完这段代码你会得到两个最重要的输出相机内参矩阵mtx和畸变系数dist。mtx矩阵通常长这样[[fx, 0, cx], [ 0, fy, cy], [ 0, 0, 1]]其中fx,fy是焦距像素单位(cx, cy)是主点通常接近图像中心。这个mtx就是我们后续所有三维计算的基础它描述了摄像头本身的成像几何特性。记得把这个矩阵保存下来比如用numpy.save下次就不用重新标定了。3. 特征提取与匹配让计算机找到“相同”的点拿到了相机的“体检报告”我们就可以开始处理我们想要重建的物体照片了。假设我们对着一个玩具小车水平移动摄像头拍了两张照片。现在计算机的挑战是它怎么知道第一张照片里的车灯和第二张照片里的车灯是同一个点3.1 特征点图像的“关键地标”这就是特征提取与匹配要解决的问题。我们不可能比较每一个像素那样计算量太大而且像素本身很容易受光照、角度影响。我们需要找到一些稳定、独特的“关键地标”也就是特征点。好的特征点就像一座城市里的标志性建筑比如埃菲尔铁塔无论从哪个角度看你都能认出来。原始文章提到了SIFT、SURF和ORB三种算法。我这里多说几句我的使用感受SIFT精度最高最稳定但速度最慢SURF可以看作是SIFT的加速版在保持较好稳定性的同时速度快了很多ORB最快而且是免费的SIFT和SURF过去有专利问题但有时在纹理不明显或视角变化大时匹配效果会差一些。对于入门和大多数情况我推荐使用SURF它在速度和精度上取得了很好的平衡。不过由于专利原因在OpenCV的默认安装中SURF可能无法直接使用这也是为什么我们选择了特定版本。3.2 动手实现特征匹配让我们看看代码具体是怎么做的。我们读取两张待重建的图片然后让SURF探测器去“寻找地标”。def epipolar_geometric(Images_Path, K): IMG glob.glob(Images_Path) img1, img2 cv2.imread(IMG[0]), cv2.imread(IMG[1]) img1_gray cv2.cvtColor(img1, cv2.COLOR_BGR2GRAY) img2_gray cv2.cvtColor(img2, cv2.COLOR_BGR2GRAY) # 创建SURF检测器 # 注意在新版OpenCV中可能需要使用 cv2.SURF_create()并设置Hessian阈值 SURF cv2.xfeatures2d_SURF.create() # 同时检测关键点并计算描述子。描述子是一个向量可以理解为该特征点的“指纹” keypoint1, descriptor1 SURF.detectAndCompute(img1_gray, None) keypoint2, descriptor2 SURF.detectAndCompute(img2_gray, None) print(f第一张图找到 {len(keypoint1)} 个特征点第二张图找到 {len(keypoint2)} 个) # 使用暴力匹配器Brute-Force Matcher进行匹配 # cv2.NORM_L2 是计算描述子之间距离的度量方式适用于SURF/SIFT bf cv2.BFMatcher(cv2.NORM_L2, crossCheckTrue) matches bf.match(descriptor1, descriptor2) # 按距离排序距离越小匹配越可靠 matches sorted(matches, keylambda x: x.distance) print(f初步匹配到 {len(matches)} 对点)运行到这里计算机会给你找出成百上千对“可能”的匹配点。但这里面有很多是错的就像把东方明珠误认成了广州塔。所以我们需要一个“去伪存真”的步骤。3.3 用RANSAC剔除误匹配误匹配是三维重建失败的主要原因之一。我们使用RANSAC随机抽样一致算法和单应性矩阵来过滤掉那些“不合群”的坏点。RANSAC的思想很聪明它随机选几对匹配点计算一个变换模型比如单应性矩阵然后看有多少其他点符合这个模型。符合的点叫“内点”不符合的叫“外点”。重复这个过程很多次最后采用内点最多的那个模型。# 将匹配点对转换为NumPy数组格式 src_pts np.asarray([keypoint1[m.queryIdx].pt for m in matches]) dst_pts np.asarray([keypoint2[m.trainIdx].pt for m in matches]) # 使用RANSAC算法和单应性矩阵来寻找内点掩码 # 参数100.0是判断点是否为内点的阈值像素距离 H, mask cv2.findHomography(src_pts, dst_pts, cv2.RANSAC, 100.0) # mask是一个数组值为1代表是内点好匹配0代表外点误匹配 print(f经过RANSAC过滤后保留 {mask.sum()} 对可靠匹配点) # 只保留内点 points1 src_pts[mask.ravel() 1] points2 dst_pts[mask.ravel() 1]经过这一步我们得到了一组相对干净、可靠的匹配点对。这些点对是连接两张二维图像的桥梁也是我们推算三维空间的唯一依据。你可以把匹配结果画出来看看好的匹配应该看起来整齐、有规律大部分匹配点都落在物体的同一特征部位上。4. 核心三角测量从二维到三维的“一跃”这是最激动人心也最需要数学理解的一步。我们现在有了1. 两张图片上对应点的像素坐标points1,points22. 相机的内参矩阵K。我们的目标是求出每一个匹配点在世界坐标系中的三维坐标。4.1 对极几何与本质矩阵两张图片之间的几何关系由对极几何描述。想象一下空间中的一个点X它在两个相机位置成像在两个像平面上的投影点分别是x1和x2。这三个点X,相机1中心,相机2中心构成一个平面称为极平面。这个约束关系可以用一个神奇的矩阵——本质矩阵 E来表示x2^T * E * x1 0。本质矩阵E包含了两个相机之间的旋转和平移信息E [t]_x * R其中[t]_x是平移向量的叉乘矩阵。我们的代码里compute_essential_normalized函数就是用归一化的8点法来求解这个E矩阵。求解之前需要先用相机内参K的逆矩阵将像素坐标points1和points2转换到归一化相机坐标系下的坐标points1n和points2n消除内参的影响。# 将像素坐标转换到归一化相机坐标系相当于使用了一个理想相机 points1n np.dot(np.linalg.inv(K), points1.T) # 注意转置使形状为 3 x n points2n np.dot(np.linalg.inv(K), points2.T) # 计算本质矩阵E E compute_essential_normalized(points1n, points2n) print(计算得到的本质矩阵 E:) print(E)4.2 从本质矩阵恢复相机姿态得到E矩阵后我们可以通过奇异值分解SVD从中恢复出两个相机之间的相对旋转R和平移t。但是这里有一个尺度不确定性和四重歧义的问题。简单说我们只知道平移的方向不知道具体移动了多远尺度而且R和t有四种可能的组合。compute_P_from_essential(E)这个函数就列出了这四种可能的第二个相机的投影矩阵P2假设第一个相机矩阵P1 [I | 0]即位于世界坐标系原点。我们需要从这四种可能中选出正确的那一个。判断标准很直观正确的相机姿态应该使得三维点位于两个相机的前方即Z坐标为正。代码里通过reconstruct_one_point函数试探性地三角化一个点然后检查它在两个相机坐标系下的深度从而选出正确的P2。P1 np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0]]) # 第一个相机矩阵 [I | 0] P2s compute_P_from_essential(E) # 四个可能的第二个相机矩阵 correct_index -1 for i, P2 in enumerate(P2s): # 三角化一个点来测试 d1 reconstruct_one_point(points1n[:, 0], points2n[:, 0], P1, P2) # 将点变换到第二个相机坐标系下 P2_homogenous np.linalg.inv(np.vstack([P2, [0, 0, 0, 1]])) d2 np.dot(P2_homogenous[:3, :4], d1) # 如果该点在两个相机坐标系下的深度Z值都为正则是正确的姿态 if d1[2] 0 and d2[2] 0: correct_index i break # 得到正确的第二个相机矩阵并调整格式 P2 np.linalg.inv(np.vstack([P2s[correct_index], [0, 0, 0, 1]]))[:3, :4]4.3 线性三角测量得到三维点云一旦确定了两个相机的投影矩阵P1和P2对于每一对匹配点x1和x2我们就可以通过线性三角测量来求解它们对应的三维空间点X。原理是x1 P1 * X,x2 P2 * X。这构成了一个超定线性方程组AX 0可以通过SVD求解最小二乘解。linear_triangulation函数就是为所有匹配点对完成这个计算。Points3D linear_triangulation(points1n, points2n, P1, P2) # Points3D 的形状是 (4, n)是齐次坐标需要转换为非齐次 (3, n) Points3D Points3D / Points3D[3, :] # 归一化 Points3D Points3D[:3, :] # 取前三维最后我们用 Matplotlib 的 3D 绘图功能将Points3D画出来你就能看到从二维图片中“飞”出来的三维点云了虽然这个点云可能稀疏、有噪声但它确确实实是三维的你可以旋转视角从各个方向观察它。5. 效果优化与实战经验分享第一次跑通代码看到三维点云时你可能会有点失望点云稀疏、形状扭曲、甚至是一团乱麻。别灰心单目重建的效果非常依赖于输入图片的质量和拍摄方式。根据我多次实验的经验以下几点能极大提升成功率第一拍摄素材是关键。这是最重要的一步。务必保持相机纯平移运动最好是水平或垂直移动尽量避免旋转和倾斜。物体和背景的纹理要丰富一个纯色的光滑物体比如白色石膏像对单目重建来说是“隐形”的。光照要均匀稳定避免闪光灯和移动的阴影。我建议先从有明显纹理的物体开始比如一本封面丰富的书、一个毛绒玩具。第二参数调优有技巧。在特征匹配阶段cv2.findHomography中的ransacReprojThreshold参数代码中的100.0很关键。它决定了判断一个匹配点是否为内点的距离阈值。如果重建结果点很少可以适当调大这个值比如到200-300如果点很多但很杂乱可以调小这个值比如到50。SURF检测器的Hessian阈值也可以调整阈值越高检测到的特征点越少但越稳定。第三尝试不同的特征探测器。如果SURF效果不好可以试试ORB。虽然ORB是二值描述子用cv2.NORM_HAMMING距离进行匹配但在某些场景下速度极快且效果不错。也可以尝试AKAZE。把代码里的SURF部分换成下面这样# 使用ORB orb cv2.ORB_create(nfeatures5000) # 增加特征点数量 keypoint1, descriptor1 orb.detectAndCompute(img1_gray, None) keypoint2, descriptor2 orb.detectAndCompute(img2_gray, None) # 匹配器也要改用汉明距离 bf cv2.BFMatcher(cv2.NORM_HAMMING, crossCheckTrue)第四后续处理提升观感。得到原始点云后你可以使用像Open3D或PCL这样的库进行点云后处理。比如用统计滤波或半径滤波去除离群噪点用泊松重建等表面重建算法将点云转换成网格模型。这能让你的重建结果从“一堆散点”升级成“一个模型”。单目三维重建是一个入门计算机视觉三维感知的绝佳项目。它让你亲手触摸到从图像到几何的核心管道。虽然它的精度和鲁棒性比不上双目或多目视觉更无法与激光雷达相比但它的简单性和低成本赋予了它巨大的魅力。当你用自己的代码把手机拍的照片变成可以旋转的三维物体时那种感觉就像掌握了某种创造世界的魔法。我建议你把整个流程跑通后试着重建你桌面的小物件不断调整拍摄方法和参数观察点云的变化。这个过程里积累的直觉对你理解更复杂的视觉SLAM、NeRF等前沿技术都会有莫大的帮助。