【立体匹配】双目相机外参自标定方法介绍

双目相机外参自标定方法是一种无需固定标定板，在拍摄实际场景的两张图像时，通过计算两幅图像之间的匹配特征点对，结合相机的内参矩阵，来实时求解两个相机之间相对位置（即外参）的方法。

原理

双目相机外参自标定基于双目立体视觉原理，通过匹配两幅图像中的对应点，利用这些匹配点对以及相机的内参矩阵，计算出两个相机之间的旋转矩阵和平移向量，即外参矩阵。这种方法避免了使用固定标定板的繁琐过程，提高了标定的灵活性和实时性。

具体步骤：

1. 畸变校正：
   首先，利用已知的相机内参矩阵对两幅原始图像进行畸变校正。这是因为在相机成像过程中，由于镜头制造和安装误差等原因，图像会产生畸变，影响后续的特征点匹配和参数计算。

2. 特征点提取与匹配：
   使用特征点检测算法（如SIFT、SURF、ORB等）在两幅校正后的图像中分别提取特征点。
   然后，利用特征点匹配算法（如FLANN、BFMatcher等）找出两幅图像中的匹配点对。这些匹配点对是两幅图像中同一物点在两个相机成像平面上的投影。

3. 本质矩阵求解：
   利用匹配点对和相机内参矩阵，通过数学方法求解本质矩阵（Essential Matrix）。本质矩阵描述了两个相机坐标系之间的旋转和平移关系，但不包含相机的内参信息。
   在OpenCV中，可以使用cv::findEssentialMat函数来求解本质矩阵。该函数通常结合随机抽样一致性算法（RANSAC）来提高求解的鲁棒性，减少误匹配点对的影响。

4. 外参矩阵分解：
   通过对本质矩阵进行奇异值分解（SVD）等操作，可以分解出两个相机之间的旋转矩阵和平移向量，即外参矩阵。
   在OpenCV中，可以使用cv::recoverPose函数从本质矩阵中恢复出相机的旋转矩阵和平移向量。该函数同样可以结合RANSAC算法来提高结果的准确性。

实践

外参自标定步骤：

第一步：使用左右相机的对应匹配点。

第二步：依据双目极线约束条件构造本征矩阵关于关键点位置的代价函数，将匹配关键点对带入代价函数，所求的代价函数之和称之为能量函数。

第三步：使用非线性优化方法优化本征矩阵的值，使得能量函数最小，此时估计的本征矩阵为最终结果，将得到的本征矩阵分解，即可得到双目相机的相对位姿。

左右相机的匹配点：

下图依次为原标定参数极线校正标定板视差图、恢复参数的标定板视差图：

下图依次为原标定参数极线校正人脸深度图、使用一段时间后校正深度图、恢复参数后校正深度图：

以下是著名的八点法求解本质矩阵：

# 读取左右图像  
imgL = cv2.imread('left.jpg',0) # 左图像  
imgR = cv2.imread('right.jpg',0) # 右图像  
  
# 初始化SIFT检测器  
sift = cv2.xfeatures2d.SIFT_create()  
  
# 检测并计算关键点和描述符  
kpL, desL = sift.detectAndCompute(imgL,None)  
kpR, desR = sift.detectAndCompute(imgR,None)  
  
# 匹配描述符  
bf = cv2.BFMatcher()  
matches = bf.knnMatch(desL,desR,k=2)  
  
# 选择好的匹配点  
good = []  
for m,n in matches:  
    if m.distance < 0.75*n.distance:  
        good.append(m)  
  
# 如果匹配点数量不够，则退出程序  
if len(good)<4:  
    print("Not enough matches are found - {}/{}".format(len(good),4))  
    exit()  
  
# 提取匹配点的坐标  
src_pts = np.float32([ kpL[m.queryIdx].pt for m in good ]).reshape(-1,1,2)  
dst_pts = np.float32([ kpR[m.trainIdx].pt for m in good ]).reshape(-1,1,2)  
  
# 利用非线性优化求解本质矩阵和旋转矩阵  
E, mask = cv2.findEssentialMat(src_pts, dst_pts, None, method=cv2.RANSAC, prob=0.999, maxIters=1000)  
R, t, mask = cv2.recoverPose(E, src_pts, dst_pts)  
print("Rotation matrix: \n" + str(R))  
print("Translation vector: \n" + str(t))

但求解结果随缘，待优化。

以下是利用Scipy库进行非线性优化：

def skew_matrix