一、实验原理

1.1相机标定简介

在图像测量过程以及机器视觉应用中，为确定空间物体表面某点的三维几何位置与其在图像中对应点之间的相互关系，必须建立相机成像的几何模型，这些几何模型参数就是相机参数。在大多数条件下这些参数必须通过实验与计算才能得到，这个求解参数的过程就称之为相机标定（或摄像机标定）。无论是在图像测量或者机器视觉应用中，相机参数的标定都是非常关键的环节，其标定结果的精度及算法的稳定性直接影响相机工作产生结果的准确性。因此，做好相机标定是做好后续工作的前提，提高标定精度是科研工作的重点所在。在零失真相机标定法中，每一个像素矢量在空间都是独立标定的，无需知道相机内部的结构，无需建立几何模型。

1.2算法流程

（1）打印一张棋盘方格图并贴在一个平面上
（2）从不同角度拍摄若干张模板图像
（3）检测出图像中的特征点
（4）由检测到的特征点计算出每幅图像中的平面投影矩阵H
（5）确定出摄像机的参数

相机标定方法有：传统相机标定法、主动视觉相机标定方法、相机自标定法、零失真相机标定法。

传统相机标定法需要使用尺寸已知的标定物，通过建立标定物上坐标已知的点与其图像点之间的对应，利用一定的算法获得相机模型的内外参数。根据标定物的不同可分为三维标定物和平面型标定物。三维标定物可由单幅图像进行标定，标定精度较高，但高精密三维标定物的加工和维护较困难。平面型标定物比三维标定物制作简单，精度易保证，但标定时必须采用两幅或两幅以上的图像。传统相机标定法在标定过程中始终需要标定物，且标定物的制作精度会影响标定结果。同时有些场合不适合放置标定物也限制了传统相机标定法的应用。

目前出现的自标定算法中主要是利用相机运动的约束。相机的运动约束条件太强，因此使得其在实际中并不实用。利用场景约束主要是利用场景中的一些平行或者正交的信息。其中空间平行线在相机图像平面上的交点被称为消失点，它是射影几何中一个非常重要的特征，所以很多学者研究了基于消失点的相机自标定方法。自标定方法灵活性强，可对相机进行在线定标。但由于它是基于绝对二次曲线或曲面的方法，其算法鲁棒性差。

基于主动视觉的相机标定法是指已知相机的某些运动信息对相机进行标定。该方法不需要标定物，但需要控制相机做某些特殊运动，利用这种运动的特殊性可以计算出相机内部参数。基于主动视觉的相机标定法的优点是算法简单，往往能够获得线性解，故鲁棒性较高，缺点是系统的成本高、实验设备昂贵、实验条件要求高，而且不适合于运动参数未知或无法控制的场合。

零失真相机标定法是以LCD显示屏为参考基准，以相移光栅为媒介，建立LCD像素与相机传感器像素之间的映射关系，确定每个相机像素点在LCD上的视点位置。镜头使相机在LCD上的视场为非矩形。在这个有畸变的视场内，可以构造一个内接的虚拟传感器，并保持相同的像素数。这样每个虚拟像素点就一定落在某四个相邻视点构成的任意四边形之内。虚拟像素点的亮度将由这四点的亮度经加权插值确定，而与其它像素点无关。用这四个加权系数的集合对原始图像作重采样（四次乘法、四次加法），就可以得到零畸变的输出图。对彩色相机的RGB三通道分别处理，但选用一个公共的虚拟传感器，则合成的彩色图像将是零畸变、零色差的。每个像素点（物理视点和虚拟点）的位置误差都是零均值，均方差都可小于1/1,000像素点距。

1.3相机标定原理

世界坐标系(world coordinate system)
用户定义的三维世界的坐标系，为了描述目标物在真实世界里的位置而被引入，坐标 (xw, yw, zw)，单位为m。
相机坐标系(camera coordinate system)
在相机上建立的坐标系，为了从相机的角度描述物体位置而定义，作为沟通世界坐标系和图像/像素坐标系的中间一环，坐标 (xc, yc, zc)，单位为m。
图像坐标系(image coordinate system)
也有看到称作“像平面坐标系”的，为了描述成像过程中物体从相机坐标系到图像坐标系的投影透射关系而引入，方便进一步得到像素坐标系下的坐标，坐标 (x,y)，单位为m。
像素坐标系(pixel coordinate system)
为了描述物体成像后的像点在数字图像上的坐标而引入，是我们真正从相机内读取到的信息所在的坐标系，坐标 (u,v)，单位为pixels（像素数目）

针孔相机模型

图像径向畸变是图像像素点以畸变中心为中心点,沿着径向产生的位置偏差,从而导致图像中所成的像发生形变。图像径向畸变是成像过程中最主要的畸变,同时也是对成像效果影响最大的畸变。图像的径向畸变的存在不仅影响了图像视觉效果,而且增加了后续识别、跟踪等图像处理算法的度。

1.4矩阵计算

外参矩阵：

计算单应性矩阵H

设三维世界坐标的点为，二维相机平面像素坐标为，所以标定用的棋盘格平面到图像平面的单应性关系为： (其中，K为相机的内参矩阵，R为外部参数矩阵（旋转矩阵），T为平移向量。令：

设棋盘格位于Z=0的平面，定义旋转矩阵R的第i列为 ri, 则有：

于是空间到图像的映射可改为：H=λK[r1 r2 t]
其中H 是描述Homographic矩阵,可通过最小二乘，从角点世界坐标到图像坐标的关系求解。

计算内参数矩阵

根据步骤1中的式子，令 H 为 H = [h1 h2 h3]，则 [h1 h2 h3]=λK[r1 r2 t]，再根据正交和归一化的约束可以得到等式：

即每个单应性矩阵能提供两个方程，而内参数矩阵包含5个参数，要求解，至少需要3个单应性矩阵。为了得到三个不同的单应性矩阵，我们使用至少三幅棋盘格平面的图片进行标定。通过改变相机与标定板之间的相对位置来得到三个不同的图片。为了方便计算，我们定义：

B 中的未知量可表示为6D 向量 b，

设H中的第i列为 hi，，根据b的定义，可以推导出公式，

最后推导出：

通过上式，我们可知当观测平面 n ≥ 3 时，即至少3幅棋盘格图像，可以得到b的唯一解，求得相机内参数矩阵K。

1.5相机标定参数求解

二、代码及实现运行

2.1代码实现

生成棋盘:

import cv2
import sys
import numpy as np
image = np.ones([1080, 1920, 3], np.uint8) * 255
x_nums = 14                                             
y_nums = 7                                 
square_pixel = 120   # 1080/9 = 120 pixels
x0 = square_pixel
y0 = square_pixel
def DrawSquare():
    flag = -1                                        
    for i in range(y_nums):
        flag = 0 - flag
        for j in range(x_nums):
            if flag > 0:
                color = [0,0,0]                         
            else:
                color = [255,255,255]
            cv2.rectangle(image,(x0 + j*square_pixel,y0 + i*square_pixel),
                          (x0 + j*square_pixel+square_pixel,y0 + i*square_pixel+square_pixel),color,-1)
            flag = 0 - flag
    cv2.imwrite('D:/chess_map_14x7.bmp',image)
if __name__ == '__main__':
    DrawSquare()

相机标定：

import cv2
import numpy as np
import glob
#criteria:角点精准化迭代过程的终止条件
criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 30, 0.001)
#棋盘格模板规格
len = 18#黑白格长度
w = 6
h = 13
# 世界坐标系中的棋盘格点,例如(0,0,0), (1,0,0), (2,0,0) ....,(8,5,0)，去掉Z坐标，记为二维矩阵
world_point = np.zeros((w*h, 3), np.float32)#初始化一个6*13行3列的矩阵，类型为float
#把每一个世界坐标的xy赋值[:,:2]中，：表示列表的所有字列表全部，：2表示从1截止到第二个数字
world_point[:, :2] = np.mgrid[0:w*len:len, 0:h*len:len].T.reshape(-1, 2)
# 储存棋盘格角点的世界坐标和图像坐标对
world_points = []  # 在世界坐标系中的三维点
imgpoints = []  # 在图像平面的二维点
j=1
images = glob.glob('*.jpg')#读取所有jpg文件
for fname in images:
    img = cv2.imread(fname)
    #将图片缩小
    img = cv2.resize(img,None,fx=0.4, fy=0.4, interpolation = cv2.INTER_CUBIC)
    gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
    # 找到棋盘格角点
    #寻找角点，存入corners，ret是找到角点的flag（如果找到角点则为true）
    ret, corners = cv2.findChessboardCorners(gray, (w, h), None)
    # 如果找到足够点对，将其存储起来
    if ret is True:
        cv2.cornerSubPix(gray, corners, (11, 11), (-1, -1), criteria)
        #加入到世界坐标系
        world_points.append(world_point)
        #角点加入到图像坐标系
        imgpoints.append(corners)
        # 将角点在图像上显示
        cv2.drawChessboardCorners(img, (w, h), corners, ret)    
        cv2.imshow('findCorners', img)
        cv2.waitKey(1)
        print(str(j)+"完成\n")
        j+=1
    else:
        print('错误')
        j+=1
cv2.destroyAllWindows()
#求解摄像机的内在参数和外在参数。mtx 内参数矩阵，dist 畸变系数，rvecs 旋转向量，tvecs 平移向量 
ret, mtx, dist, rvecs, tvecs = cv2.calibrateCamera(world_points, imgpoints, gray.shape[::-1], None, None)
# 找棋盘格角点
# 去畸变
img2 = cv2.imread('1.jpg')
img2= cv2.resize(img2,None,fx=0.4, fy=0.4, interpolation = cv2.INTER_CUBIC)
cv2.imshow("畸变图像", img2)
cv2.waitKey(100)
cv2.destroyAllWindows()
h,  w = img2.shape[:2]
#内参数矩阵
newcameramtx, roi = cv2.getOptimalNewCameraMatrix(mtx, dist, (w, h), 0, (w, h))  # 自由比例参数
#校正后图像
dst = cv2.undistort(img2, mtx, dist, None, newcameramtx)
cv2.imwrite('calibresult.png', dst)
f = open('内参数矩阵.txt', 'w+')
f.write('内参数矩阵:\n'+str(newcameramtx)+'\n')
f.close()
print('内参数输出完毕')
f = open('外参数矩阵.txt', 'w+')
#得到旋转矩阵与平移
for t in range(0,j-1):#向量变为矩阵
    newrvecs=cv2.Rodrigues(rvecs[t].ravel())
    newtvecs=cv2.Rodrigues(tvecs[t].ravel())
    f.write('第'+str(t+1)+'张图像的外参数矩阵:\n'+'(1)旋转矩阵：\n'+str(newrvecs[0])+'\n'+'(2)平移矩阵：\n'+str(newtvecs[0])+'\n')
    t+=1
f.close()
print('外参数输出完毕')
f = open('畸变系数.txt', 'w+')
f.write('畸变系数:\n'+str(dist)+'\n')
f.close()
print('结果输出完毕')
np.savez('data.npz', mtx= mtx, dist= dist,rvecs=rvecs,tvecs=tvecs)
print('*.npz文件输出完毕')

三维立方体：

import cv2
import numpy as np
import glob
#加载数据
with np.load('data.npz') as X:
#加载上一部生成的参数
mtx, dist, _, _= [X[i] for i in ('mtx','dist','rvecs','tvecs')]
#criteria:角点精准化迭代过程的终止条件
criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 30, 0.001)
#棋盘格模板规格
len = 18#黑白格长度
w = 6
h = 13
#画坐标轴和立方体
def draw(img, corners, imgpts,imgpts2):
    corner=tuple(corners[0].ravel())#ravel()方法将数组维度拉成一维数组
   # img要画的图像，corner起点，tuple终点，颜色，粗细
    img= cv2.line(img, corner, tuple(imgpts2[0].ravel()), (255,0,0), 8)
    img= cv2.line(img, corner, tuple(imgpts2[1].ravel()), (0,255,0), 8)
    img= cv2.line(img, corner, tuple(imgpts2[2].ravel()), (0,0,255), 8)
    font = cv2.FONT_HERSHEY_SIMPLEX
    cv2.putText(img,'X',tuple(imgpts2[0].ravel()+2), font, 1,( 255,0,0),2,cv2.LINE_AA)
    cv2.putText(img,'Y',tuple(imgpts2[1].ravel()+2), font, 1,( 0,255,0),2,cv2.LINE_AA)
    cv2.putText(img,'Z',tuple(imgpts2[2].ravel()+2), font, 1,( 0,0,255),2,cv2.LINE_AA)
    imgpts= np.int32(imgpts).reshape(-1,2)#draw ground floor in green
    for i,j in zip(range(4),range(4,8)):#正方体顶点逐个连接
        img= cv2.line(img, tuple(imgpts[i]), tuple(imgpts[j]),(255,215,0),3)#draw top layer in red color  
    #imgpts[4:]是八个顶点中上面四个顶点
    #imgpts[:4]是八个顶点中下面四个顶点
    #用函数drawContours画出上下两个盖子，它的第一个参数是原始图像，第二个参数是轮廓，一个python列表，第三个参数是轮廓的索引（当设置为-1时绘制所有轮廓）
    img = cv2.drawContours(img, [imgpts[4:]],-1,(255,215,0),3)
    img = cv2.drawContours(img, [imgpts[:4]],-1,(255,215,0),3)
    return img
objp= np.zeros((w*h,3), np.float32)
objp[:,:2] = np.mgrid[0:w*len:len,0:h*len:len].T.reshape(-1,2)
axis = np.float32([[0,0,0], [0,2*len,0], [2*len,2*len,0], [2*len,0,0],
[0,0,-2*len],[0,2*len,-2*len],[2*len,2*len,-2*len],[2*len,0,-2*len] ])
axis2= np.float32([[3*len,0,0], [0,3*len,0], [0,0,-3*len]]).reshape(-1,3)
images = glob.glob('*.jpg')
i=1;
for fname in images:
    img = cv2.imread(fname)
    img = cv2.resize(img,None,fx=0.4, fy=0.4, interpolation = cv2.INTER_CUBIC)
    gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
    # 找到棋盘格角点
    #寻找角点，存入corners，ret是找到角点的flag
    ret, corners = cv2.findChessboardCorners(gray, (w, h), None)
    if ret is True:
        corners2= cv2.cornerSubPix(gray,corners,(11,11),(-1,-1),criteria)
        #求解物体位姿的需要
        _,rvecs, tvecs, inliers =cv2.solvePnPRansac(objp, corners2, mtx, dist)
        #projectPoints()根据所给的3D坐标和已知的几何变换来求解投影后的2D坐标。
        #imgpts是整体的8个顶点
        imgpts,_=cv2.projectPoints(axis, rvecs, tvecs, mtx, dist)
        #imgpts2是三个坐标轴的x,y,z划线终点
        imgpts2, _ =cv2.projectPoints(axis2, rvecs, tvecs, mtx, dist)
        img=draw(img,corners2,imgpts,imgpts2)
        cv2.imshow('世界坐标系与小盒子',img)
        cv2.imwrite(str(i)+'.png', img)
        cv2.waitKey(1000)
        i+=1;
cv2.destroyAllWindows()
print("完毕")