目标

• 学习摄像机畸变以及摄像机的内部参数和外部参数
• 根据摄像机相关参数对畸变图像进行修复

基础说明 

        今天的低价单孔摄像机(照相机)会给图像带来很多畸变。畸变主要有两种:径向畸变和切向畸变。如下图所示用红色直线将棋盘的两个边标注出来，但是你会发现棋盘的边界并不和红线重合。所有我们认为应该是直线的也都凸出来了。

这种畸变可以通过下面的方程组进行纠正:

        于此相似，另外一个畸变是切向畸变，这是由于透镜与成像平面不可能绝对平行造成的。这种畸变会造成图像中的某些点看上去的位置会比我们认为的位置要近一些。它可以通过下列方程组进行校正：

        简单来说，如果我们想对畸变的图像进行校正就必须找到五个造成畸变的系数:

        除此之外，我们还需要再找到一些信息，比如摄像机的内部和外部参数。内部参数是摄像机特异的。它包括的信息有焦距(fx, fy)，光学中心（cx, cy）等。这也被称为摄像机矩阵。它完全取决于摄像机自本身，只需要计算算一次以后就可以已知使用了。可以用下面的3x3 的矩阵表示：

        外部参数与旋转和变换向量相对应，它可以将3D点的坐标转换到坐标系系统中。 

        在3D 相关应用中，必须先校正这些畸变。为了找到这些参数，我们必须先提供一些包含明显图案模式的样本图片(比如说棋盘)。我们可以在上面找到一些特殊点(如棋盘的四个角点）。我们起到这些特殊点在图片中的位置以及它们的真实位置。有了这些信息，我们就可以使用数学方法求解畸变系数。这就是整个故事的摘要了。为了得到更好的结果我们至少需要10 个这样的图案模式。

代码 

         如上所述，我们至少需要10 图案模式来进行摄像机标定。OpenCV 自带了一些棋盘图像(/sample/cpp/left001.jpg--left14.jpg), 所以我们可以使用它们。为了便于理解我们可以认为仅有一张棋盘图像。重要的是在进行摄像机标定时我们需要输入一组3D 真实世界中的点以及与它们对应2D 图像中的点。2D 图像的点可以在图像中很容易的找到。(这些点在图像中的位置是棋盘上两个黑色方块相互接触的地方）

        那么真实世界中的3D 的点呢？这些图像来源与静态摄像机和棋盘不同的摆放位置和朝向。所以我们需要知道(X,Y,Z)的值。但是为了简单,我们可以说棋盘在XY 平面是静止的(所以Z 总是等于0)摄像机在围着棋盘移动。这种假设我们只需要知到(X,Y) 的值就可以了。现在为了求X,Y 的值,我们只需要传入这些点（0,0）,(1,0),(2,0)...它们代表了点的位置。在这个例子中，我们的结果的单位就是棋盘(单个)方块的大小。但是如果我们知道单个方块的大小(假如􄖣30mm），我们输入的值就可以是(0,0),(30,0),（60,0）...结果的单位就是mm。在本例中我们不知道方块的大小(因为不是我们拍的)所以只能用前一种方法了。3D 点被称为对象点，2D 图像点被称为图像点。

 设置

        为了找到棋盘的图案，我们要使用函数cv2.findChessboardCorners()。我们还需要传入图案的类型，比如 8x8 的格子或5x5 的格子等。在本例中我们使用的7x8 的格子。(常情况下棋盘都是8x8 或者7x7)。它会返回角点，如果得到图像的返回值类型(Retval)就会是True。这些角点会按顺序排列(从左到右,从上到下)。

        在找到这些角点之后我们可以使用函数cv2.cornerSubPix() 增加准确度。我们使用函数cv2.drawChessboardCorners() 绘制图案。所有的这些步骤都被包含在下面的代码中了：

import numpy as np
import cv2
import glob

# termination criteria
criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 30, 0.001)

# prepare object points, like (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)

# Arrays to store object points and image points from all the images.
objpoints = [] # 3d point in real world space
imgpoints = [] # 2d points in image plane.

images = glob.glob('*.jpg')

for fname in images:
    img = cv2.imread(fname)
    gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
    # Find the chess board corners
    ret, corners = cv2.findChessboardCorners(gray, (7,6),None)
    # If found, add object points, image points (after refining them)
    if ret == True:
        objpoints.append(objp)
        corners2 = cv2.cornerSubPix(gray,corners,(11,11),(-1,-1),criteria)
        imgpoints.append(corners2)
        # Draw and display the corners
        img = cv2.drawChessboardCorners(img, (7,6), corners2,ret)
        cv2.imshow('img',img)
        cv2.waitKey(500)
cv2.destroyAllWindows()

一副图像和被绘制在上边的图案：

标定 

         在得到了这些对象点和图像点之后，我们已经准备好来做摄像机标定了。我们要使用的函数是cv2.calibrateCamera()。它会返回摄像机矩阵，畸变系数、旋转和变换向量等。

ret, mtx, dist, rvecs, tvecs = cv2.calibrateCamera(objpoints, imgpoints, gray.shape[::-1],None,None)

畸变校正 

        现在我们找到我们想要的东西了，我们可以找到一幅图像来对它进行校正。OpenCV 提供两种方法供我们都学习一下。不过在那之前我们可以使用从函数cv2.getOptimalNewCameraMatrix() 得到的自由缩放系数对摄像机矩阵进行优化。如果缩放系数alpha = 0，返回的非畸变图像会带有最少量的不想要的像素。它甚至有可能在图像角点去除一些像素。如果alpha = 1，所有的像素都会被返回，还有一些黑图像。它还会返回一个ROI 图像，我们可以用来对结果进行裁剪。我们读取一个新的图像(left2.ipg):

img = cv2.imread('left12.jpg')
h, w = img.shape[:2]
newcameramtx, roi=cv2.getOptimalNewCameraMatrix(mtx,dist,(w,h),1,(w,h))

        使用cv2.undistort() 这是最简单的方法。只需使用这个函数和上面得到的ROI 对结果进行裁剪。 

# undistort
dst = cv2.undistort(img, mtx, dist, None, newcameramtx)
# crop the image
x,y,w,h = roi
dst = dst[y:y+h, x:x+w]
cv2.imwrite('calibresult.png',dst)

        使用remapping 这应该属于“曲线救国”了。首先我们要找到从畸变图像到非畸变图像的映射方程。再使用映射方程。 

# undistort
mapx,mapy = cv2.initUndistortRectifyMap(mtx,dist,None,newcameramtx,(w,h),5)
dst = cv2.remap(img,mapx,mapy,cv2.INTER_LINEAR)
# crop the image
x,y,w,h = roi
dst = dst[y:y+h, x:x+w]
cv2.imwrite('calibresult.png',dst)

这两中方法给出的结果是相同的。结果如下所示：

        你会发现结果图像中所有的边界都变直了。
        现在我们可以使用Numpy 提供写函数(np.savez,np.savetxt 等)将摄像机矩阵和畸变系数保存以便以后使用。 

反向投影误差

        我们可以利用反向投影误差对我们找到的参数的准确性进行估计。得到的结果越接近0越好。有了内部参数、畸变参数和旋转变换矩，我们就可以使用cv2.projectPoints() 将对象点转换到图像点。然后就可以计算变换得到图像与角点检测算法的绝对差了。然后我们计算所有标定图像的误差平均值。

mean_error = 0
for i in xrange(len(objpoints)):
    imgpoints2, _ = cv2.projectPoints(objpoints[i], rvecs[i], tvecs[i], mtx, dist)
    error = cv2.norm(imgpoints[i],imgpoints2, cv2.NORM_L2)/len(imgpoints2)
    tot_error += error
print ("total error: ", mean_error/len(objpoints))