计算机视觉----相机标定模型

简介

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

---

前言

摄像机标定简单来说是从世界坐标系转换为相机坐标系，再由相机坐标系转换为图像坐标系的过程，也就是求最终的投影矩阵P的过程
相机标定的目的和意义在于我们所处的世界是三维的，而照片是二维的，这样我们可以把相机认为是一个函数，输入量是一个场景，输出量是一幅灰度图。这个从三维到二维的过程的函数是不可逆的。而相机标定的目标是我们找一个合适的数学模型，求出这个模型的参数，这样我们能够近似这个三维到二维的过程，使这个三维到二维的过程的函数找到反函数。

张正友相机标定方法介于传统标定法和自标定法之间，但克服了传统标定法需要的高精度标定物的缺点，而仅需使用一个打印出来的棋盘格就可以。同时也相对于自标定而言，提高了精度，便于操作。因此张氏标定法被广泛应用于计算机视觉方面。

---

提示：以下是本篇文章正文内容，下面案例可供参考

一、相机标定简介

相机标定指建立相机图像像素位置与场景点位置之间的关系，根据相机成像模型，由特征点在图像中坐标与世界坐标的对应关系，求解相机模型的参数。相机需要标定的模型参数包括内部参数和外部参数。

针孔相机成像原理其实就是利用投影将真实的三维世界坐标转换到二维的相机坐标上去，其模型示意图如下图所示：

从图中我们可以看出，在世界坐标中的一条直线上的点在相机上只呈现出了一个点，其中发生了非常大的变化，同时也损失和很多重要的信息，这正是我们3D重建、目标检测与识别领域的重点和难点。实际中，镜头并非理想的透视成像，带有不同程度的畸变。理论上镜头的畸变包括径向畸变和切向畸变，切向畸变影响较小，通常只考虑径向畸变。

径向畸变：径向畸变主要由镜头径向曲率产生（光线在远离透镜中心的地方比靠近中心的地方更加弯曲）。导致真实成像点向内或向外偏离理想成像点。

• 枕形畸变：畸变像点相对于理想像点沿径向向外偏移，远离中心；
• 桶状畸变：径向畸点相对于理想点沿径向向中心靠拢。
  
• 数学公式
  x ~ K[R|t]X
  
  x为相机中的坐标；
  X为真实世界坐标；
  [R|t]为外参矩阵,
  K为内参矩阵，是相机内部参数组成的一个3*3的矩阵，
  其中，f代表焦距；s为畸变参数；(x0,y0)为中心点坐标；a为纵横比例参数，我们可以默认设为1，所以x0=ay0。
  [R|t]为外参矩阵，R是描述照相机方向的旋转矩阵，t是描述照相机中心位置的三维平移向量。

二、张友正黑白棋盘标定

1.思想

张正友标定法利用黑白棋盘格标定板，在得到一张标定板的图像之后，用Harris角点检测算法得到每一个角点的像素坐标 (u,v) 。张正友标定法将世界坐标系固定于棋盘格上，则棋盘格上任一点的物理坐标W=0，由于标定板的世界坐标系是人为事先定义好的，标定板上每一个格子的大小是已知的，我们可以计算得到每一个角点在世界坐标系下的物理坐标(U,V,W=0) 。我们将利用这些信息：每一个角点的像素坐标 、每一个角点在世界坐标系下的物理坐标，来进行相机的标定，获得相机的内外参矩阵、畸变参数。

2.原理

3.模型求解

• 计算单应性矩阵
  单应性：在计算机视觉中被定义为一个平面到另一个平面的投影映射。首先确定，图像平面与标定物棋盘格平面的单应性。
  设棋盘格位于Z=0,定义旋转矩阵R的第i列为ri, 则有：
  
  则空间到图像的映射可变为：
  

其中H 是描述Homographic矩阵,可通过最小二乘，从角点世界坐标到图像坐标的关系求解。
计算内参数矩阵，令H= [h1 h2 h3 ]

Homography有8个自由度，通过上述等式的矩阵运算，根据r1和r2正交，以及归一化的约束可以得到如下等式：

定义：

其中：

B是对称阵，其未知量可表示为6D向量b：

设H为第i列中的hi:

根据b的定义，可以推出如下公式：

即可推导出：

内部参数可通过如下公式计算(cholesky分解)：

3.外部参数求解

4.极大似然估计
给定 n张棋盘格图像，每张图像有 m 个角点，最小化下述公式等同于极大似然估计：

上述非线性优化问题可以利用Levenberg-Marquardt 算法求解，需要初值K,Ri,ti|i=1…n

5.径向畸变估计
张氏标定法只关注了影响最大的径向畸变。则数学表达式为：

其中，(u,v)是理想无畸变的像素坐标，(u,v)(u,v)是实际畸变后的像素坐标。(u0,v0)代表主点，(x,y)是理想无畸变的连续图像坐标，(x,y)(x,y)是实际畸变后的连续图像坐标。k1和k2为前两阶的畸变参数。

化作矩阵形式：

记做：Dk=d，则可得：

计算得到畸变系数k。使用最大似然的思想优化得到的结果，即像上一步一样，LM法计算下列函数值最小的参数值(K为相机内参矩阵A)：

三、实验内容及过程

3.1 实验要求

1.打印一张棋盘格A4纸张（黑白间距已知），并贴在一个平板上
2.针对棋盘格拍摄若干张图片（一般10-20张）
3.在图片中检测特征点（Harris角点）
4.根据角点位置信息及图像中的坐标，求解Homographic矩阵
5.利用解析解估算方法计算出5个内部参数，以及6个外部参数
6.根据极大似然估计策略，设计优化目标并实现参数的refinement

3.2 实验数据及环境

1.实验数据

手机（华为荣耀V20 后置：4800万像素）
拍摄各个角度下打印的黑白棋盘格图片13张：

拍摄相同角度下打印的黑白棋盘格图片13张：

棋盘格规格为10*7，每个小格边长为0.038m

2.实验环境

win 10 ，64位，python3.6.5 ，OpenCV 4.2.0

3.3 实现代码

# -*- coding:utf-8 -*-
import cv2
import glob
import numpy as np
from ipython_genutils.py3compat import xrange
import matplotlib.pyplot as plt
from pylab import *

from Calibration.work1.Harris import plot_harris_points
# from refine_all_param import refinall_all_param

cbraw = 9  # 有6行角点
cbcol = 6  # 有4列角点

objp = np.zeros((cbraw * cbcol, 3), np.float32)
'''
设定世界坐标下点的坐标值，因为用的是棋盘可以直接按网格取；
假定棋盘正好在x-y平面上，这样z值直接取0，简化初始化步骤。
mgrid把列向量[0:cbraw]复制了cbcol列，把行向量[0:cbcol]复制了cbraw行。
转置reshape后，每行都是4×6网格中的某个点的坐标。
'''
objp[:, :2] = np.mgrid[0:cbraw, 0:cbcol].T.reshape(-1, 2)

objpoints = []  # 3d point in real world space
imgpoints = []  # 2d points in image plane.

# glob是个文件名管理工具
images = glob.glob("images1/*.jpg")
for fname in images:
	# 对每张图片，识别出角点，记录世界物体坐标和图像坐标
	img = cv2.imread(fname)  # source image

	# 我用的图片太大，缩小了一半
	# img = cv2.resize(img, None, fx=0.5, fy=0.5, interpolation=cv2.INTER_CUBIC)
	gray1 = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)  # 转灰度


	# 寻找角点，存入corners，ret是找到角点的flag
	ret, corners = cv2.findChessboardCorners(gray1, (9, 6), None)
	img = cv2.drawChessboardCorners(gray1, (9,6), corners,ret)
	cv2.waitKey(1000)
	# criteria:角点精准化迭代过程的终止条件
	criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 30, 0.001)
	# 执行亚像素级角点检测
	corners2 = cv2.cornerSubPix(gray1, corners, (11, 11), (-1, -1), criteria)
	objpoints.append(objp)
	imgpoints.append(corners2)
	# 在棋盘上绘制角点,只是可视化工具
	img = cv2.drawChessboardCorners(gray1, (9, 6), corners2, ret)
	# cv2.imshow('img', img)
# cv2.waitKey(1000)

'''
传入所有图片各自角点的三维、二维坐标，相机标定。
每张图片都有自己的旋转和平移矩阵，但是相机内参和畸变系数只有一组。
mtx，相机内参；dist，畸变系数；外参数：revcs，旋转矩阵；tvecs，平移矩阵。
'''
ret, mtx, dist, rvecs, tvecs = cv2.calibrateCamera(objpoints, imgpoints, gray1.shape[::-1], None, None)
img = cv2.imread('images1/empire1.jpg')
# extrinsics_param\
rvecs = np.array(rvecs)
# print("rvecs="+str(rvecs))
# rvecsnp, _ = cv2.Rodrigues(rvecs)

# print("rvecsnp="+str(rvecsnp))
# print("tvecs="+str(tvecs))
# 注意这里跟循环开头读取图片一样，如果图片太大要同比例缩放，不然后面优化相机内参肯定是错的。
# img = cv2.resize(img, None, fx=0.5, fy=0.5, interpolation=cv2.INTER_CUBIC)
h, w = img.shape[:2]

'''
优化相机内参（camera matrix），这一步可选。
参数1表示保留所有像素点，同时可能引入黑色像素，
设为0表示尽可能裁剪不想要的像素，这是个scale，0-1都可以取。
'''
newcameramtx, roi = cv2.getOptimalNewCameraMatrix(mtx, dist, (w, h), 1, (w, h))#显示更大范围的图片（正常重映射之后会删掉一部分图像）

# 纠正畸变
dst = cv2.undistort(img, mtx, dist, None, newcameramtx)

# 这步只是输出纠正畸变以后的图片
x, y, w, h = roi
dst = dst[y:y + h, x:x + w]
cv2.imwrite('calibresult2.png', dst)

# 打印我们要求的两个矩阵参数
"""参数newcameramtx为内参数矩阵，
dist为畸变系数，total error为总体误差。
newcameramtx:
 [[1.25045215e+03 0.00000000e+00 1.48972123e+03]
 [0.00000000e+00 7.76162048e+02 1.31742293e+03]
 [0.00000000e+00 0.00000000e+00 1.00000000e+00]]
dist:
 [[ 2.49509032e-02  3.66095655e-01 -1.75899580e-02 -5.07301401e-04
  -8.18362820e-01]]
total error:  0.10671872187032518
从上面的运行结果可以看出，
归一化焦距fx=1.25045215e+03，fy=7.76162048e+02，
像主点（光心）的坐标cx=1.48972123e+03，cy=1.31742293e+03。
所有图像投影坐标和亚像素角点坐标之间的总体的平均误差大概为0.10，误差较小，标定结果不错"""

print("newcameramtx:\n", newcameramtx)
print("dist:\n", dist)
print("ret:\n",ret)
print("mtx:\n",mtx)
# 计算误差
tot_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: ", tot_error / len(objpoints))

绘画角点的代码（可用opencv写好的函数）

# 引入所有函数，绘图板块
from PIL import Image
from imutils.feature import harris
from pylab import *
from numpy import *
# n维图像库，我们这里使用filters函数
from scipy.ndimage import filters


# 使用高斯滤波器进行卷积，标准差为3
# 这个函数得到Harris响应函数值的图像
def compute_harris_response(im, sigma=3):
	""" Compute the Harris corner detector response function
		for each pixel in a graylevel image. """

	# derivatives
	imx = zeros(im.shape)
	# x方向上的高斯导数
	filters.gaussian_filter(im, (sigma, sigma), (0, 1), imx)
	imy = zeros(im.shape)
	# y方向上的高斯导数
	filters.gaussian_filter(im, (sigma, sigma), (1, 0), imy)

	# compute components of the Harris matrix
	# 导数运算的高斯模糊值，第一个参数表示input，通过这样的方式得到矩阵的分量
	Wxx = filters.gaussian_filter(imx * imx, sigma)
	Wxy = filters.gaussian_filter(imx * imy, sigma)
	Wyy = filters.gaussian_filter(imy * imy, sigma)

	# 计算特征值与迹
	Wdet = Wxx * Wyy - Wxy ** 2
	Wtr = Wxx + Wyy

	return Wdet / Wtr


# 从响应图像中获得harris角点检测结果，min_dist表示最少数目的分割角点，threshold是阈值
def get_harris_points(harrisim, min_dist=10, threshold=0.6):
	""" Return corners from a Harris response image
		min_dist is the minimum number of pixels separating
		corners and image boundary. """

	# find top corner candidates above a threshold
	# 寻找高于阈值的候选角点，.max是numpy的函数·
	corner_threshold = harrisim.max() * threshold
	harrisim_t = (harrisim > corner_threshold) * 1

	# get coordinates of candidates
	# nonzeros(a)返回数组a中值不为零的元素的下标，它的返回值是一个长度为a.ndim(数组a的轴数)的元组,.T是矩阵转置操作
	coords = array(harrisim_t.nonzero()).T

	# ...and their values，harrisim的响应值
	candidate_values = [harrisim[c[0], c[1]] for c in coords]

	# 从小到大输出，注意输出的是下标
	index = argsort(candidate_values)

	# 标记分割角点范围的坐标
	allowed_locations = zeros(harrisim.shape)
	allowed_locations[min_dist:-min_dist, min_dist:-min_dist] = 1

	# select the best points taking min_distance into account,选择harris点
	filtered_coords = []
	for i in index:
		if allowed_locations[coords[i, 0], coords[i, 1]] == 1:
			# 添加坐标
			filtered_coords.append(coords[i])
			# 删除过近区域
			allowed_locations[(coords[i, 0] - min_dist):(coords[i, 0] + min_dist),
			(coords[i, 1] - min_dist):(coords[i, 1] + min_dist)] = 0

	return filtered_coords


# 绘制检测角点，灰度图 ,上个函数的角点需要满足的条件，1：角点范围之内；2：高于harrisim响应值的阈值
def plot_harris_points(image, filtered_coords):
	""" Plots corners found in image. """

	figure()
	gray()
	imshow(image)
	# 将坐标以*标出
	plot([p[1] for p in filtered_coords],
		 [p[0] for p in filtered_coords], 'ro')
	axis('off')
	show()


# 读入图像
for i in range(1,14):
	print("第"+str(i)+"张图")
	im = array(Image.open('images1/empire'+str(i)+'.jpg').convert('L'))
	# 检测harris角点
	harrisim = compute_harris_response(im)
	# figure()
	# imshow(harrisim)
	# axis('off')
	# threshold = [0.1]
	# for j, thres in enumerate(threshold):
	filtered_coords = get_harris_points(harrisim, 6)
	plot_harris_points(im, filtered_coords)
# show()

3.4 实验结果

1、不同角度拍摄的角点

内参矩阵：

畸变系数：

旋转向量：

rvecs: [[[-0.38331118]
  [ 0.48723626]
  [ 1.71109437]]

 [[-0.45411493]
  [-0.03231818]
  [-0.01473678]]

 [[-0.04345208]
  [-0.78809805]
  [-3.0115678 ]]

 [[-0.01145075]
  [ 0.40085887]
  [ 0.33333145]]

 [[-0.54167611]
  [ 0.4791829 ]
  [ 1.48189766]]

 [[-0.44199928]
  [ 0.25740441]
  [ 0.91299673]]

 [[-0.0463189 ]
  [-0.63133509]
  [-1.62688687]]

 [[-0.55588811]
  [-0.18347327]
  [-0.87782993]]

 [[-0.05706544]
  [-0.01193983]
  [-1.47530636]]

 [[-0.43935382]
  [ 0.38784614]
  [ 1.51912972]]

 [[-0.4176884 ]
  [-0.48684335]
  [-1.47771514]]

 [[-0.29551665]
  [ 0.12034187]
  [-1.08950363]]

 [[ 0.05496156]
  [-0.49431   ]
  [-1.75260836]]]

平移向量：

tvecs: [array([[ 2.7148961 ],
       [-4.16079681],
       [16.68959021]]), array([[-3.64450399],
       [-2.18762748],
       [15.65285177]]), array([[ 3.81963211],
       [ 1.40552132],
       [13.25435986]]), array([[-3.73334425],
       [-3.12204986],
       [15.13478662]]), array([[ 2.76260801],
       [-4.17157634],
       [16.87731711]]), array([[-1.0071051 ],
       [-5.13343633],
       [17.25474235]]), array([[-0.08193984],
       [ 2.87492183],
       [12.42362365]]), array([[-4.71138451],
       [ 0.09700443],
       [14.22382484]]), array([[-3.36180182],
       [ 4.22473216],
       [13.81892064]]), array([[ 2.35659495],
       [-4.60907617],
       [17.35177962]]), array([[-2.51863381],
       [ 1.8754995 ],
       [11.18001005]]), array([[-4.63085747],
       [ 2.23809406],
       [15.74990388]]), array([[-0.08244932],
       [ 4.17750851],
       [13.37648316]])]

总误差：

2、相同角度拍摄的角点

内参数矩阵：

畸变系数：

旋转向量：

rvecs: [[[-0.08251195]
  [-0.00774483]
  [ 1.58220765]]

 [[-0.0980028 ]
  [ 0.0286664 ]
  [ 1.58236595]]

 [[-0.15378806]
  [ 0.06036396]
  [ 1.55295778]]

 [[-0.1474118 ]
  [ 0.0454044 ]
  [ 1.55619334]]

 [[-0.16825274]
  [ 0.12941535]
  [ 1.5465853 ]]

 [[-0.05994238]
  [-0.01253248]
  [ 1.57987408]]

 [[-0.1113285 ]
  [ 0.0430649 ]
  [ 1.58810691]]

 [[-0.10699552]
  [ 0.04157332]
  [ 1.58185485]]

 [[-0.10110703]
  [ 0.03530747]
  [ 1.58014635]]

 [[-0.10359539]
  [ 0.03666896]
  [ 1.58478425]]

 [[-0.10140888]
  [ 0.03361181]
  [ 1.58035485]]

 [[-0.10091694]
  [ 0.03114847]
  [ 1.58719659]]

 [[-0.10052098]
  [ 0.03258916]
  [ 1.58223094]]]

平移向量：

[array([[ 2.22253207],
       [-4.26875302],
       [18.17558029]]), array([[ 2.02899116],
       [-4.47806431],
       [19.06540396]]), array([[ 2.03943056],
       [-4.80539584],
       [19.4636413 ]]), array([[ 2.03003722],
       [-5.01055769],
       [19.39913314]]), array([[ 1.83214688],
       [-5.0281785 ],
       [19.93167163]]), array([[ 1.94609342],
       [-4.62603081],
       [18.55547549]]), array([[ 2.09854521],
       [-4.27777057],
       [19.31213953]]), array([[ 1.99358885],
       [-4.33545851],
       [19.23956555]]), array([[ 1.98220393],
       [-4.3656387 ],
       [19.15736314]]), array([[ 2.07126341],
       [-4.3729312 ],
       [19.1971351 ]]), array([[ 2.00080433],
       [-4.4686167 ],
       [19.12629221]]), array([[ 2.09121347],
       [-4.43094166],
       [19.15081471]]), array([[ 2.00591686],
       [-4.42331201],
       [19.11459241]])]

总误差：total error: 0.13716845175700373

四、总结

• 代码中的行列数据是自己自作的棋盘格的行列各-1
• 注意尺寸的单位是m
• 当图片像素过高时，运行代码会报错
• 标定效果不好与棋盘格是否标准和平整有关，表面褶皱，弯曲的棋盘格，实验效果会比标准的差很多
• 拍照尽量不要有反光的现象，因为反光的地方成像偏白色，影响标定效果

坐标系转化学习链接