相机校准和3D重建1

相机校准

基础

一些针孔相机会给图像带来明显的失真，两种主要的变形是径向变形和切向变形。
径向变形会导致直线出现弯曲。距图像中心越远，径向畸变越大。例如，下面显示一个图像，其中棋盘的两个边缘用红线标记；但是，可以看到棋盘的边框不是直线，并且与红线不匹配，所有预期的直线都凸出。

径向变形可以表示成如下：
$$\begin{gathered} x_{distorted}=x(1+k_1{r}^2+k_2{r}^4+k_3{r}^6) \\ {y}_{distorted}=y(1+k_1{r}^2+k_2{r}^4+k_3{r}^6) \end{gathered}$$
同样，由于摄像镜头未完全平行于成像平面对齐，因此会发生切向畸变。因此，图像中的某些区域看起来可能比预期的要近。切向畸变的量可以表示为：
$$\begin{gathered} x_{distorted}=x+[2p_{1}xy+p_2(r^2+2x^2)] \\ {y}_{distorted}=y+[p_1(r^2+2y^2)+2p_{2}xy] \end{gathered}$$
简而言之，需要找到五个参数，称为失真系数，公式如下：
$$Distortioncoefficients=(k_1{k}_2{p}_1{p}_2{k}_3)$$
除此之外，还需要其他一些信息，例如相机的内在和外在参数。内部参数特定于摄像机，它们包括诸如焦距$(f_x,f_y)$和光学中心$(c_x,c_y)$之类的信息。焦距和光学中心可用于创建相机矩阵，该相机矩阵可用于消除由于特定相机镜头而引起的畸变。相机矩阵对于特定相机而言是唯一的，因此一旦计算出，就可以在同一相机拍摄的其他图像上重复使用。它表示为3x3矩阵：
$$cameramatrix=\left[\begin{array}{ccc}{f}_x& 0& c_x\\ 0& f_y& c_y\\ 0& 0& 1\end{array}\right]$$
外在参数对应于旋转和平移矢量，其将3D点的坐标平移为坐标系。
对于立体声应用，首先需要纠正这些失真。要找到这些参数，必须提供一些定义良好的图案的示例图像（例如国际象棋棋盘）。找到一些已经知道其相对位置的特定点（例如棋盘上的四角）。知道现实世界空间中这些点的坐标，也知道图像中的坐标，因此可以求解失真系数。为了获得更好的结果，至少需要10个测试模式。

处理流程及代码

如上所述，相机校准至少需要10个测试图案，OpenCV附带了一些国际象棋棋盘的图像（请参见samples/data/left01.jpg–left14.jpg），将利用这些图像。考虑棋盘的图像，相机校准所需的重要输入数据是3D现实世界点集以及图像中这些点的相应2D坐标，可以从图像中轻松找到2D图像点（这些图像点是国际象棋棋盘中两个黑色正方形相互接触的位置）。
真实世界中的3D点如何处理？这些图像是从静态相机拍摄的，而国际象棋棋盘放置在不同的位置和方向。因此，需要知道$(X,Y,Z)$值，但是为简单起见，可以说棋盘在XY平面上保持静止（因此Z始终为0），并且照相机也相应地移动了。这种考虑有助于仅找到X，Y值。现在对于X，Y值，可以简单地将点传递为（0,0），（1,0），（2,0），…，这表示点的位置。在这种情况下，得到的结果将是棋盘正方形的大小比例。但是，如果知道正方形大小（例如30毫米），则可以将值传递为（0,0），（30,0），（60,0），…。因此，得到的结果以毫米为单位（在这种情况下，不知道正方形的大小，因为没有拍摄那些图像，因此以正方形的大小进行传递）。
3D点称为对象点，而2D图像点称为图像点。

开始

要在国际象棋棋盘中查找图案，可以使用函数cv.findChessboardCorners()。还需要传递所需的图案，例如8x8网格，5x5网格等。在下面示例中，使用7x6网格（通常，棋盘有8x8的正方形和7x7的内部角）。它返回角点和retval，如果获得图案，则为True，这些角将按顺序放置（从左到右，从上到下）。
此功能可能无法在所有图像中找到所需的图案。因此，一个不错的选择是编写代码，使它启动相机并检查每帧所需的图案，获得图案后，找到角并将其存储在列表中。另外，在阅读下一帧之前请提供一些时间间隔，以便我们可以在不同方向上调整棋盘，继续此过程，直到获得所需数量的良好图案为止。即使在此处提供的示例中，也不确定给出的14张图像中有多少张是好的。
因此，必须阅读所有图像并仅拍摄好图像。 除了棋盘，还可以使用圆形网格，在这种情况下，必须使用函数cv.findCirclesGrid()来找到模式。 较少的图像足以使用圆形网格执行相机校准。
一旦找到拐角，就可以使用cv.cornerSubPix()来提高其精度，还可以使用cv.drawChessboardCorners()绘制图案。

校准

有了目标点和图像点，现在可以进行校准了。可以使用函数cv.calibrateCamera()返回相机矩阵，失真系数，旋转和平移矢量等。

不失真

可以拍摄图像并对其进行扭曲，OpenCV提供了两种方法来执行此操作。但是首先，可以使用cv.getOptimalNewCameraMatrix()基于自由缩放参数来优化相机矩阵。如果缩放参数alpha = 0，则返回具有最少不需要像素的未失真图像；因此，它甚至可能会删除图像角落的一些像素。如果alpha = 1，则所有像素都保留有一些额外的黑色图像。此函数还返回可用于裁剪结果的图像ROI。

1.使用cv.undistort()

这是最简单的方法，只需调用该函数并使用上面获得的ROI裁剪结果即可。

2.使用remapping

该方式有点困难。首先，找到从扭曲图像到未扭曲图像的映射函数，然后使用重映射功能。

两种方法都给出相同的结果，使所有边缘都是笔直的。可以使用NumPy中的写入功能(np.savez，np.savetxt等)存储相机矩阵和失真系数，以备将来使用。

重投影误差

重投影误差可以很好地估计找到的参数的精确程度，重投影误差越接近零，发现的参数越准确。给定固有，失真，旋转和平移矩阵，必须首先使用cv.projectPoints()将对象点转换为图像点，然后，可以计算出通过变换得到的绝对值和拐角发现算法之间的绝对值范数。为了找到平均误差，计算为所有校准图像计算的误差的算术平均值。
参见下面完整示例代码cam_calibration.py

import cv2 as cv
import numpy as np
import glob

# 终止条件
criteria = (cv.TERM_CRITERIA_EPS + cv.TermCriteria_MAX_ITER, 30, 0.001)
# 准备对象点,如(0, 0, 0), (1, 0, 0), (2, 0, 0) ..., (6, 5, 0)
objp = np.zeros((6 * 7, 3), np.float32)
objp[:,:2] = np.mgrid[0:7, 0:6].T.reshape(-1, 2)
# 用于存储所有图像的对象点和图像点的数组
objpoints = []  #真实世界中的3d点
imgpoints = []  #图像中的2d点
images = glob.glob('./OpenCV/data/left*.jpg')
for fname in images:
    img = cv.imread(fname)
    gray = cv.cvtColor(img, cv.COLOR_BGR2GRAY)
    # 找到棋盘角落
    ret, corners = cv.findChessboardCorners(gray, (7, 6), None)
    # 如果找到,添加对象点,图像点(细化之后)
    if ret == True:
        objpoints.append(objp)
        corners2 = cv.cornerSubPix(gray, corners, (11, 11), (-1, -1), criteria)
        imgpoints.append(corners)
        # 绘制并显示拐角
        cv.drawChessboardCorners(img, (7, 6), corners2, ret)
        cv.imshow('img', img)
        cv.waitKey(500)
cv.destroyAllWindows()
ret, mtx, dist, rvecs, tvecs = cv.calibrateCamera(objpoints, imgpoints,
                                                    gray.shape[::-1], None, None)
np.savez('./OpenCV/data/B.npz', mtx = mtx, dist = dist, rvecs = rvecs, tvecs = tvecs)                                                    
img = cv.imread('./OpenCV/data/left12.jpg')
h, w = img.shape[:2]
newcamerax, roi = cv.getOptimalNewCameraMatrix(mtx, dist, (w, h), 1, (w, h))
# undistort
dst1 = cv.undistort(img, mtx, dist, None, newcamerax)
#剪裁图像
x, y, w, h = roi
dst1 = dst1[y:y + h, x:x + w]
cv.imwrite('./OpenCV/data/calibresult1.png', dst1)
# remapping
mapx, mapy = cv.initUndistortRectifyMap(mtx, dist, None, newcamerax, (w, h), 5)
dst2 = cv.remap(img, mapx, mapy, cv.INTER_LINEAR)
#剪裁图像
x, y, w, h = roi
dst2 = dst2[y:y + h, x:x + w]
cv.imwrite('./OpenCV/data/calibresult2.png', dst2)
mean_error = 0
for i in range(len(objpoints)):
    imgpoints2, _ = cv.projectPoints(objpoints[i], rvecs[i], tvecs[i], mtx, dist)
    error = cv.norm(imgpoints[i], imgpoints2, cv.NORM_L2)/len(imgpoints2)
    mean_error += error
print("total error: {}".format(mean_error/len(objpoints)))

total error: 0.023686000375385676

一张上面画有图案的图像如下所示：

calibresult1.png

calibresult2.png

姿态估计

基础及流程

上面相机校准中，发现了相机矩阵，失真系数等。给定图案图像，可以利用以上信息来计算其姿势或物体在空间中的位置，例如其旋转方式， 对于平面物体，可以假设Z = 0，这样，问题就变成了如何将相机放置在空间中以查看图案图像。 因此，如果知道对象在空间中的位置，则可以在其中绘制一些2D图以模拟3D效果。
看看如何做，问题是，想在棋盘的第一个角上绘制3D坐标轴(X,Y,Z)。 X轴为蓝色，Y轴为绿色，Z轴为红色。 因此，实际上Z轴应该感觉像它垂直于棋盘平面。
首先，从先前的校准结果中加载相机矩阵和失真系数。创建一个函数draw，该函数将棋盘上的角（使用cv.findChessboardCorners()获得）和轴点绘制为3D轴。
然后，与前面的情况一样，创建终止条件，对象点（棋盘上角的3D点）和轴点。 轴点是3D空间中用于绘制轴的点。 绘制长度为3的轴（由于根据该棋盘尺寸进行了校准，因此单位将以国际象棋正方形的尺寸为单位）。因此我们的X轴从(0,0,0)绘制为(3,0,0)，对于Y轴、Z轴，从(0,0,0)绘制为(0,0,-3)，负号表示它被拉向相机。
之后加载每个图像。搜索7x6网格，如果找到，将使用子角像素对其进行优化。然后使用函数cv.solvePnPRansac()计算旋转和平移，一旦有了这些变换矩阵，就可以使用它们将轴点投影到图像平面上。简而言之，在图像平面上找到与3D空间中(3,0,0),(0,3,0),(0,0,3)中的每一个相对应的点，一旦获得它们，就可以使用draw()函数从第一个角到这些点中的每个点绘制线条。

绘制立方体

如果要绘制立方体，可参照draw()函数创建drawcube()函数和修改轴点为如下3D空间中多维数据集的8个角。

axis = np.float32([[0,0,0], [0,3,0], [3,3,0], [3,0,0], [0,0,-3], [0,3,-3],[3,3,-3],[3,0,-3] ])

完整示例代码如下：
pose_estimate.py

import cv2 as cv
import numpy as np
import glob

# 加载先前保存的数据
with np.load('./OpenCV/data/B.npz') as X:
    mtx, dist, _, _ = [X[i] for i in ('mtx', 'dist', 'rvecs', 'tvecs')]
def draw(img, corners, imgpts):
    corner = tuple(corners[0].ravel())
    corner = tuple(map(int, corner))
    img = cv.line(img, corner, tuple(map(int, tuple(imgpts[0].ravel()))), (255, 0, 0), 5)
    img = cv.line(img, corner, tuple(map(int, tuple(imgpts[1].ravel()))), (0, 255, 0), 5)
    img = cv.line(img, corner, tuple(map(int, tuple(imgpts[2].ravel()))), (0, 0, 255), 5)
    return img
def drawcube(img, corners, imgpts):
    imgpts = np.int32(imgpts).reshape(-1, 2)
    # 用绿色绘制底层
    img = cv.drawContours(img, [imgpts[:4]], -1, (0, 255, 0), -3)
    # 用蓝色绘制高
    for i, j in zip(range(4), range(4, 8)):
        img = cv.line(img, tuple(imgpts[i]), tuple(imgpts[j]), (255), 3)
    # 用红色绘制顶层
    img = cv.drawContours(img, [imgpts[4:]], -1, (0, 0, 255), 3)
    return img
criteria = (cv.TermCriteria_EPS + cv.TermCriteria_MAX_ITER, 30, 0.001)
objp = np.zeros((6 * 7, 3), np.float32)
objp[:, :2] = np.mgrid[0:7, 0:6].T.reshape(-1, 2)
axis = np.float32([[3,0,0], [0,3,0], [0,0,-3]]).reshape(-1,3)
axis2 = np.float32([[0,0,0], [0,3,0], [3,3,0], [3,0,0], [0,0,-3],
                        [0,3,-3],[3,3,-3],[3,0,-3]])
for fname in glob.glob('./Opencv/data/left*.jpg'):
    img = cv.imread(fname)
    gray = cv.cvtColor(img,cv.COLOR_BGR2GRAY)
    ret, corners = cv.findChessboardCorners(gray, (7,6),None)
    if ret == True:
        corners2 = cv.cornerSubPix(gray,corners,(11,11),(-1,-1),criteria)
        # 找到旋转和平移矢量。
        ret,rvecs, tvecs = cv.solvePnP(objp, corners2, mtx, dist)
        # 将3D点投影到图像平面
        imgpts, jac = cv.projectPoints(axis, rvecs, tvecs, mtx, dist)
        img = draw(img, corners2, imgpts)
        cv.imshow('img', img)
        imgpts2, jac2 = cv.projectPoints(axis2, rvecs, tvecs, mtx, dist)
        img2 = drawcube(img, corners2, imgpts2)
        cv.imshow('img2', img2)
        k = cv.waitKey(0) & 0xFF
        if k == ord('s'):
            cv.imwrite('./OpenCV/' + fname[:6] + '.png', img)
            cv.imwrite('./OpenCV/c' + fname[:6] + '.png', img2)
cv.destroyAllWindows()

当使用corner = tuple(corners[0].ravel())和img = cv.line(img, corner, tuple(imgpts[0].ravel()), (255,0,0), 5)时出现

cv2.error: OpenCV(4.8.1) :-1: error: (-5:Bad argument) in function 'line'
> Overload resolution failed:
>  - Can't parse 'pt1'. Sequence item with index 0 has a wrong type
>  - Can't parse 'pt1'. Sequence item with index 0 has a wrong type

错误，网上参考方法是将cv.line()中使用的两个端点转化为int，打印corner和tuple(imgpts[1].ravel())的值发现确实为浮点型小数。使用int(corner)、int(tuple(imgpts[1].ravel()))提示

TypeError: int() argument must be a string, a bytes-like object or a number, not 'tuple'

可使用示例中代码corner = tuple(map(int, corner))，将tuple中的值转换为int型。


学习来源：OpenCV-Python中文文档