刚开始学习关于双目立体视觉系统三维重建方面的知识，大都是采取了很多大佬的文章，所以写篇文章做个总结。目前只是找出一个可执行的方案，不是说去掌握其中的每一个细节，掌握完细节后，再去落地，关于精确掌握各个部分，是之后进一步的研究，而不是目前所做的事。在完成项目的过程中，若遇到问题，就带着问题去解决，这样效率高。
此篇文章会涉及到相关程序的编写，其中的相关细节，比如说变量的问题，在编写的过程中需要搞清楚，不然会出现问题。
首先进行双目标定之前，先看怎样能够在连接了双目相机的情况下驱动左右两个相机（通过usb），注意在下载opencv-python与opencv-contrib-python两个包时，一定下载与python的版本相对应的包版本，以下提供关于如何下载的文章以及两个包的下载地址：
(18条消息) Python3.6下安装opencv_考古学家lx(李玺)的博客-CSDN博客_python3.6对应的opencv版本
Links for opencv-python (tsinghua.edu.cn)
Links for opencv-contrib-python (tsinghua.edu.cn)
双目相机驱动程序如下：
# -*- coding: utf-8 -*-
import cv2 as cv
import time
AUTO = False # 自动拍照，或手动按s键拍照
INTERVAL = 2 # 自动拍照间隔
cv.namedWindow("left")
cv.namedWindow("right")
camera = cv.VideoCapture(1) #这里的0表示使用电脑内置摄像头，1表示使用双目相机
# 设置分辨率 左右摄像机同一频率，同一设备ID；左右摄像机总分辨率1280x480；分割为两个640x480、640x480
camera.set(cv.CAP_PROP_FRAME_WIDTH, 1280)
camera.set(cv.CAP_PROP_FRAME_HEIGHT, 480)
counter = 0
utc = time.time()
folder = "./Saveimage/" # 拍照文件目录
def shot(pos, frame):
global counter
path = folder + pos + "_" + str(counter) + ".jpg"
cv.imwrite(path, frame)
print("snapshot saved into: " + path)
while True:
ret, frame = camera.read()
# 裁剪坐标为[y0:y1, x0:x1] HEIGHT*WIDTH
left_frame = frame[0:480, 0:640]
right_frame = frame[0:480, 640:1280]
cv.imshow("left", left_frame)
cv.imshow("right", right_frame)
now = time.time()
if AUTO and now - utc >= INTERVAL:
shot("left", left_frame)
shot("right", right_frame)
counter += 1
utc = now
key = cv.waitKey(1)
if key == ord("q"):
break
elif key == ord("s"):
shot("left", left_frame)
shot("right", right_frame)
counter += 1
camera.release()
cv.destroyWindow("left")
cv.destroyWindow("right")拍照是输入法在英文状态下，按下s可以拍照，按下q退出，还有就是目前使用的双目相机不是鱼眼的，可能不需要进行畸变校正
对于双目立体视觉系统三维重建来说，分以下几个步骤：
一、双目标定
1.为什么要对双目相机进行标定？
参考链接：对双目相机进行标定的原因
2.双目标定有两种方式：
（1）opencv标定；
（2）matlab自带的标定工具箱（注意matlab里面有关于单目与双目的标定工具，别搞错了）。
在这里主要是研究关于matlab方面的标定，opencv标定没有研究，关于matlab标定的参考链接如下所示：
(18条消息) 双目摄像头Matlab参数定标_matlab双目相机标定_iNBC的博客-CSDN博客
注意：双目相机标定里，一定要删除一些误差比较大的图片，可以参考以下链接：
(18条消息) matlab双目标定（详细过程）_双目标定步骤_西海岸看日出的博客-CSDN博客
（误差很大：标定左下角有一个虚线，其中有些数据比虚线要高很多的要删掉，上述文章有说怎样删除误差大的操作）
二、图像校正
对双目相机进行标定之后，得到标定的相关数据，接下来就要进行图像的校正，图像校正方面主要是分为两个：
1.立体校正
在介绍立体校正的具体方法之前，让我们来看一下，为什么要进行立体校正？
答：双目相机系统主要的任务就是测距，而视差求距离公式是在双目系统处于理想情况下推导的，但是在现实的双目立体视觉系统中，是不存在完全的共面行对准的两个摄像机图像平面的。所以我们要进行立体校正。
立体校正的目的：把实际中非共面行对准的两幅图像，校正成共面行对准。
（共面行对准：两摄像机图像平面在同一平面上，且同一点投影到两个摄像机图像平面时，应该在两个像素坐标系的同一行）。
将实际的双目系统校正为理想的双目系统。
立体校正示例图如下所示：
2.畸变矫正
（1）为什么要进行畸变校正？（借用上面双目标定那个参考链接中一些话）
答：人最开始接触到的成像方面的知识应该是有关小孔成像的，但是由于这种成像方式只有小孔部分能透过光线就会导致物体的成像亮度很低，于是聪明的人类发明了透镜。虽然亮度问题解决了，但是新的问题又来了：由于透镜的制造工艺，会使成像产生多种形式的畸变，于是为了去除畸变（使成像后的图像与真实世界的景象保持一致），人们计算并利用畸变系数来矫正这种像差。虽然理论上可以设计出不产生畸变的透镜，但其制造工艺相对于球面透镜会复杂很多，所以相对于复杂且高成本的制造工艺，人们更喜欢用数学来解决问题。
（2）畸变类型与定义
参考链接：畸变类型与定义
3.双目相机的立体校正以及畸变矫正具体code的实现：
参考链接：Matlab标定_opencv立体校正 （这个文章的方法是基于C与C++写的）
Matlab双目标定与python-opencv配置标定参数
上述两个参考链接中，会涉及到三个函数，分别是立体校正code中的stereoRectify函数、畸变校正code中的initUndistortRectifyMap函数以及remap函数。这些函数都是针对于双目相机而言的，而对于单目相机来说，会有一些函数调用的不同。
上述三个函数中具体参数的含义的了解，提供如下三个参考链接：
（1）stereoRectify函数：https://baike.baidu.com/item/stereoRectify/1594417
（2）initUndistortRectifyMap函数：https://blog.csdn.net/u013341645/article/details/78710740
（3）remap函数：https://blog.csdn.net/yangfengman/article/details/52769716
三、视差图
一个完整的双目立体视觉系统三维重建过程，必须要做的一个部分就是画出视差图，所以需要学习双目立体匹配，这一块还在学习，有两篇文章可用于学习。
参考链接：双目立体匹配
其中第二篇文章，有一个关于整个双目立体视觉系统三维重建流程图，非常有用！
------------------------------2022年11月26日更新--------------------------------
关于立体匹配这点，在以下这篇文章中解释的非常的清楚，建议认真学习与理解：
参考链接： 立体匹配详解
-----------------------------2022年11月27日更新--------------------------------
PS：
1、建立视差图是方便点云的建立，而不是深度图的建立，要明白这一点，还有之前参考链接中的wls滤波可以看作是参考链接（立体匹配详解）中的视差后处理，其中视差后处理部分还包括二次插值的知识点。
2、视差图到深度图
参考链接：视差图到深度图的建立（不是基于python的，里面用到了C）
关于立体匹配算法有很多，考虑的方面一般有以下几点：
（1）计算量。计算量太大，程序处理就会比较长
（2）效果。视差图建立出的效果如何，比如说：所建立的视差图中是否噪声较多或者深度信息不明显（没有体现远暗近亮，物体与物体的相交处视差图没有连续光滑的变化等），一般速度快的算法，效果不太好，而速度较慢的，视差图的效果较好。
针对以上所考虑的方面，要选出适合的算法去进行立体匹配。一般的大多使用的是SGBM半全局立体匹配算法。它的特点：程序处理适中，效果也适中。 （此处存在算法创新点！）
对于SGBM算法来说，作者认为总体就包括两个步骤：一个是SGBM算法的参数设置，一个是视差计算（WTA）。
用一个例子来具体的解释立体匹配，比如说在立体匹配中使用SGBM算法，那在整个视差图的建立中，立体匹配=SGBM算法+视差后处理（各种滤波、插值等）。
四、点云重建
关于点云的三维重建，打算使用opencv自带的方法去重建，code如下：（摘取自上述视差图中参考链接中的内容）
# 将h×w×3数组转换为N×3的数组
def hw3ToN3(points):
height, width = points.shape[0:2]
points_1 = points[:, :, 0].reshape(height * width, 1)
points_2 = points[:, :, 1].reshape(height * width, 1)
points_3 = points[:, :, 2].reshape(height * width, 1)
points_ = np.hstack((points_1, points_2, points_3))
return points_
def DepthColor2Cloud(points_3d, colors):
rows, cols = points_3d.shape[0:2]
size = rows * cols
points_ = hw3ToN3(points_3d).astype(np.int16)
colors_ = hw3ToN3(colors).astype(np.int64)
# 颜色信息
blue = colors_[:, 0].reshape(size, 1)
green = colors_[:, 1].reshape(size, 1)
red = colors_[:, 2].reshape(size, 1)
# rgb = np.left_shift(blue, 0) + np.left_shift(green, 8) + np.left_shift(red, 16)
rgb = blue + green + red
# 将坐标+颜色叠加为点云数组
pointcloud = np.hstack((points_, red/255., green/255., blue/255.)).astype(np.float64)
# 删掉一些不合适的点
X = pointcloud[:, 0]
Y = pointcloud[:, 1]
Z = pointcloud[:, 2]
remove_idx1 = np.where(Z <= 0)
remove_idx2 = np.where(Z > 1000)
remove_idx3 = np.where(X > 1000)
remove_idx4 = np.where(X < -1000)
remove_idx5 = np.where(Y > 1000)
remove_idx6 = np.where(Y < -1000)
remove_idx = np.hstack((remove_idx1[0], remove_idx2[0], remove_idx3[0], remove_idx4[0], remove_idx5[0], remove_idx6[0]))
pointcloud_1 = np.delete(pointcloud, remove_idx, 0)
return pointcloud_1
以上这种方法比open3d重建速度快，处理中需要剪掉一些不合理的点云数据，其中remove就是剪掉操作。
threeD = cv2.reprojectImageTo3D(disp, camera_config.Q)
# 因为没有引入camera_config模块，所以根据自己的程序camera_config.Q改成Q
pointcloud = DepthColor2Cloud(threeD, left_remap)
# 转换为open3d的点云数据
pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(pointcloud[:,:3])
pcd.colors = o3d.utility.Vector3dVector(pointcloud[:,3:])
o3d.visualization.draw_geometries_with_editing([pcd], window_name="3D", width=1280, height=720)以上code中，camera_cofig.Q是reprojectImageTo3D参数设置中的透视变换矩阵Q，也就是深度视差映射
# 进行立体更正
R1, R2, P1, P2, Q, validPixROI1, validPixROI2 = cv2.stereoRectify(left_camera_matrix, left_distortion,
right_camera_matrix, right_distortion, size, R, T)
#上述Q-深度视差映射矩阵
矩阵，这个Q可以通过stereoRectify函数获得，也就是在立体校正中获得。
-----------2023.3.5-------------------------------
很久没有更新了，转了一圈发现这个点云重建还是得搞，所以这两天根据以上得内容以及网上其他文章，完成了双目相机关于点云重建的整个流程，基于python的程序如下：
import cv2
import numpy as np
import matplotlib.pyplot as plt
import open3d as o3d
print("---------------------------校正-----------------------------------------------------")
left_camera_matrix = np.array([[1216.225116298302, 3.659096104209817, 3.487070535684013e+02],
[0, 1.215952027618365e+03, 3.399275020161882e+02],
[0, 0, 1]])
# [1.216225116298302e+03,0,0;3.659096104209817,1.215952027618365e+03,0;3.487070535684013e+02,3.399275020161882e+02,1]
# 以上数据经过转置输入
left_distortion = np.array([[-0.912187031217956, 12.785022143809050, -0.014913596924274, 0.007976545147835, -69.494316124517500]])
# [-0.912187031217956,12.785022143809050,-69.494316124517500]: K1 K2 K3
# [-0.014913596924274,0.007976545147835]: P1 P2
# K1 K2 P1 P2 K3
right_camera_matrix = np.array([[1191.898960425957, 4.832857269551879, 3.322180942134622e+02],
[0, 1.190818757556084e+03, 3.102311949569156e+02],
[0, 0, 1]])
# [1.191898960425957e+03,0,0;4.832857269551879,1.190818757556084e+03,0;3.322180942134622e+02,3.102311949569156e+02,1]
# 以上数据经过转置输入
right_distortion = np.array([[0.207752483080191, -7.383107701152468, 0.004959717734527, 0.020783602492480, 48.655740553211850]])
# [0.207752483080191,-7.383107701152468,48.655740553211850]: K1 K2 K3
# [0.004959717734527,0.020783602492480]: P1 P2
# K1 K2 P1 P2 K3
R = np.array([
[0.999664922758816, 0.004908667067664, -0.025415491205739],
[-0.005460005450388, 0.999750284797761, -0.021669249808811],
[0.025302777438295, 0.021800757656874, 0.999442092579402],
])
# 旋转关系矩阵
# [0.999664922758816,0.004908667067664,-0.025415491205739;-0.005460005450388,0.999750284797761,-0.021669249808811;0.025302777438295,0.021800757656874,0.999442092579402]
T = np.array([-123.3666954274162, 0.846495617376818, -10.681118315039633]) # 平移关系向量
# [-1.233666954274162e+02,0.846495617376818,-10.681118315039633]
size = (640, 480) # 图像尺寸
# 进行立体更正
R1, R2, P1, P2, Q, validPixROI1, validPixROI2 = cv2.stereoRectify(left_camera_matrix, left_distortion,
right_camera_matrix, right_distortion, size, R,
T)
# 计算更正map
left_map1, left_map2 = cv2.initUndistortRectifyMap(left_camera_matrix, left_distortion, R1, P1, size, cv2.CV_16SC2)
right_map1, right_map2 = cv2.initUndistortRectifyMap(right_camera_matrix, right_distortion, R2, P2, size, cv2.CV_16SC2)
# 引入需要进行立体校正的图片
left_image_nre = cv2.imread("./Saveimage/left_0.jpg")
right_image_nre = cv2.imread("./Saveimage/right_0.jpg")
# 展示校正前的、需要进行立体校正的图片
cv2.imshow("Before Rectify1", left_image_nre)
cv2.imshow("Before Rectify2", right_image_nre)
cv2.waitKey(0)
# remap
left_image_re = cv2.remap(left_image_nre, left_map1, left_map2, cv2.INTER_LINEAR)
right_image_re = cv2.remap(right_image_nre, right_map1, right_map2, cv2.INTER_LINEAR)
cv2.imshow("After Rectify1", left_image_re)
cv2.imshow("After Rectify2", right_image_re)
cv2.imwrite("./Saveimage/left_image_re.jpg", left_image_re) # 存储图片
cv2.imwrite("./Saveimage/right_image_re.jpg", right_image_re) # 存储图片
cv2.waitKey(0) # 填0，填其他可能会卡住
"""
# 划线
for i in range(1,20):
len = 480/20
plt.axhline(y=i*len, color='r', linestyle='-')
plt.imshow(left_image_re)
plt.imshow(right_image_re)
plt.show()
"""
print("----------------------------立体匹配-----------------")
# 灰度化
imgL_gray = cv2.cvtColor(left_image_re, cv2.COLOR_BGR2GRAY)
imgR_gray = cv2.cvtColor(right_image_re, cv2.COLOR_BGR2GRAY)
# 设置参数，块大小必须为奇数(3-11)
blockSize = 5
img_channels = 2
num_disp = 16 * 8
param = {
'preFilterCap': 63, # 映射滤波器大小，默认15
"minDisparity" : 0, # 最小视差
"numDisparities" : num_disp, # 视差的搜索范围，16的整数倍
"blockSize" : blockSize,
"uniquenessRatio" : 10, # 唯一检测性参数，匹配区分度不够，则误匹配(5-15)
"speckleWindowSize" : 0, # 视差连通区域像素点个数的大小（噪声点）(50-200)或用0禁用斑点过滤
"speckleRange" : 1, # 认为不连通(1-2)
"disp12MaxDiff" : 2, # 左右一致性检测中最大容许误差值
"P1" : 8 * img_channels * blockSize** 2, # 值越大，视差越平滑，相邻像素视差+/-1的惩罚系数
"P2" : 32 * img_channels * blockSize** 2, # 同上，相邻像素视差变化值>1的惩罚系数
# 'mode': cv2.STEREO_SGBM_MODE_SGBM_3WAY
}
# # 不用wls，效果不太行
# left_matcher = cv2.StereoSGBM_create(**param)
# left_disp = left_matcher.compute(imgL_gray, imgR_gray)
# disp = cv2.normalize(left_disp, left_disp, alpha=0, beta=255, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_8U)
# cv2.imshow("depth", disp)
# cv2.waitKey(0)
# 立体匹配SGBM法
left_matcher = cv2.StereoSGBM_create(**param)
right_matcher = cv2.ximgproc.createRightMatcher(left_matcher) # 为了产生right_disp
# 视差图
left_disp = left_matcher.compute(imgL_gray, imgR_gray)
right_disp = right_matcher.compute(imgR_gray, imgL_gray) # 为了wls滤波
# wls
wls_filter = cv2.ximgproc.createDisparityWLSFilter(left_matcher)
# sigmaColor典型范围值为0.8-2.0
wls_filter.setLambda(8000.)
wls_filter.setSigmaColor(1.3)
wls_filter.setLRCthresh(24)
wls_filter.setDepthDiscontinuityRadius(3)
filtered_disp = wls_filter.filter(left_disp, imgL_gray, disparity_map_right=right_disp)
# wls深度图，下面用了公式，产生深度图
disp = cv2.normalize(filtered_disp, filtered_disp, alpha=0, beta=255, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_8U)
cv2.imshow("depth", disp)
cv2.waitKey(0)
print("-----------------------------点云重建（包含深度图到点云数据的转换）----------------------------------")
# 将h×w×3数组转换为N×3的数组
def hw3ToN3(points):
height, width = points.shape[0:2]
points_1 = points[:, :, 0].reshape(height * width, 1)
points_2 = points[:, :, 1].reshape(height * width, 1)
points_3 = points[:, :, 2].reshape(height * width, 1)
points_ = np.hstack((points_1, points_2, points_3))
return points_
def DepthColor2Cloud(points_3d, colors):
rows, cols = points_3d.shape[0:2]
size = rows * cols
points_ = hw3ToN3(points_3d).astype(np.int16)
colors_ = hw3ToN3(colors).astype(np.int64)
# 颜色信息
blue = colors_[:, 0].reshape(size, 1)
green = colors_[:, 1].reshape(size, 1)
red = colors_[:, 2].reshape(size, 1)
# rgb = np.left_shift(blue, 0) + np.left_shift(green, 8) + np.left_shift(red, 16)
rgb = blue + green + red
# 将坐标+颜色叠加为点云数组
pointcloud = np.hstack((points_, red/255., green/255., blue/255.)).astype(np.float64)
# 删掉一些不合适的点
X = pointcloud[:, 0]
Y = pointcloud[:, 1]
Z = pointcloud[:, 2]
remove_idx1 = np.where(Z <= 0)
remove_idx2 = np.where(Z > 1000)
remove_idx3 = np.where(X > 1000)
remove_idx4 = np.where(X < -1000)
remove_idx5 = np.where(Y > 1000)
remove_idx6 = np.where(Y < -1000)
remove_idx = np.hstack((remove_idx1[0], remove_idx2[0], remove_idx3[0], remove_idx4[0], remove_idx5[0], remove_idx6[0]))
pointcloud_1 = np.delete(pointcloud, remove_idx, 0)
return pointcloud_1
threeD = cv2.reprojectImageTo3D(disp, Q)
pointcloud = DepthColor2Cloud(threeD, left_image_re)
# 转换为open3d的点云数据，o3d需要加入模块open3d
# conda search xxx
# connected papers
# conda install -c open3d-admin open3d;pip install open3d-python
# 双目测距
# ！不更新open3d，报错，更新python版本的对应版本open3d
pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(pointcloud[:,:3])
pcd.colors = o3d.utility.Vector3dVector(pointcloud[:,3:])
o3d.visualization.draw_geometries_with_editing([pcd], window_name="3D", width=1280, height=720)
PS：上面程序使用的是open3d的方法进行点云重建，使用open3d的话，需要在虚拟环境里面下载open3d包，下载open3d的包的方法以及下载会遇到的问题，有一篇文章参考参考，文章如下：
open3d python版本安装及‘module‘ object has no attribute ‘read_point_cloud‘问题
上述文章还有一种下载open3d-python的方法：先打开清华大学的镜像源，在里面下载对应于虚拟环境中的python版本的oped3d-python版本的包，而清华大学关于oped3d-python的镜像源链接如下所示：
下载了包怎样在本地安装？方法：
Python3.6下安装opencv_python3.6对应的opencv版本_考古学家lx(李玺)的博客-CSDN博客
--------------------------------------------------------------------------------------------------------------------------------
2023.10.21
双目相机参数参考文章：
扫一扫
3263
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
