目录
简介
在图像测量过程以及机器视觉应用中，为确定空间物体表面某点的三维几何位置与其在图像中对应点之间的相互关系，必须建立相机成像的几何模型，这些几何模型参数就是相机参数。在大多数条件下这些参数必须通过实验与计算才能得到，这个求解参数的过程就称之为相机标定（或摄像机标定）。无论是在图像测量或者机器视觉应用中，相机参数的标定都是非常关键的环节，其标定结果的精度及算法的稳定性直接影响相机工作产生结果的准确性。因此，做好相机标定是做好后续工作的前提，提高标定精度是科研工作的重点所在。
一、相机模型
1.坐标系
- 世界坐标系：代表物体在真实世界里的三维坐标，坐标系用Xw，Yw， Zw表示。
- 相机坐标系：代表以相机光学中心为原点的坐标系，光轴与z轴重合，坐标系用Xc，Yc，Zc表示。
- 图像坐标系：代表相机拍摄图像的坐标系，原点为相机光轴与成像平面的交点，是图像的中心点，坐标系用X，Y表示。
- 像素坐标系：由于图像的基本单位是像素，所以该坐标系是图像上点在图像储存矩阵中的像素位置，坐标原点在左上角，坐标系用u，v表示。前三个坐标系的单位是毫米，而最后一个坐标系的单位是像素。
2.坐标系变化
其中Cx、Cy，是图像坐标系原点在像素坐标系下的坐标。
3.相机畸变模型
（我们一般只考虑径向畸变k和切向畸变p）
其中r2=x2+y2
- 枕形畸变：畸变像点相对于理想像点沿景象向外偏移，远离中心。
- 桶状畸变：径向畸点相对于理想点沿径向向中心靠拢。
二、相机标定原理
针对针孔相机模型， 相机标定，标定的是指内参矩阵和外参矩阵，就可以确定为一的相机模型。
- 相机的内参： 主要包含 焦距，像主点坐标，以及畸变参数即 fx fy cx cy k1 k2 k3 p1 p2
- 相机的外参： 包含 旋转矩阵和平移矩阵 即 R t
三、张正友黑白棋盘格标定
2.1.算法思想
张正友标定法利用黑白棋盘格标定板，在得到一张标定板的图像之后，用Harris角点检测算法得到每一个角点的像素坐标（u，v）。
张正友标定法将世界坐标系固定于棋盘格上，则棋盘格上任一点的物理坐标W=0，由于标定板的世界坐标系是人为事先定义好的，标定板上每一个格子的大小是已知的，我们可以计算得到每一个角点在世界坐标系下的物理坐标（U，V，W=0）。
我们将利用这些信息：每一个角点的像素坐标、每一个角点在世界坐标系下的物理坐标，来进行相机的标定，获得相机的内外参矩阵、畸变参数。
2.2.求解内参和外参的积
将世界坐标系固定于棋盘格上，则棋盘格上任一点的物理坐标W=0。因此，原单点无畸变的成像模型可以化为下式：
其中，R1，R2为旋转矩阵R的前两列。为了简便，将内参矩阵记为A。对于不同的图片，内参矩阵A为定值；对于同一张图片，内参矩阵A，外参矩阵（R1 R2 T）为定值；对于同一张图片上的单点，内参矩阵A，外参矩阵（R1 R2 T），尺度因子Z为定值。我们将A（R1 R2 T）记为矩阵H，H即为内参矩阵和外参矩阵的积，记矩阵H的三列为（H1 H2 H3），则有：
利用上式，消除尺度因子Z，可得：
此时，尺度因子Z已经被消去，因此上式对于同一张图片上所有的角点均成立。（u，v）是像素坐标系下的标定板角点的坐标，（U，V）是世界坐标系下的标定板角点的坐标。通过图像识别算法，我们可以得到标定板角点的像素坐标（u，v），又由于标定板的世界坐标系是人为定义好的，标定板上每一个格子的大小是已知的，我们可以计算得到世界坐标系下的（U，V）.
由这里的H是齐次矩阵，有8个独立未知元素。每一个标定板角点可以提供两个约束方程（u，U，V的对应关系，v，U，V的对应关系提供了两个约束方程）。因此，当一张图片上的标定板角点数量等于4时，即可求得该图片对应的矩阵H。当一张图片上的标定板角点数量大于4时，利用最小二乘法回归最佳的矩阵 H
2.3.求解内参矩阵
我们已知了矩阵H=A（R1 R2 T），接下来需要求解相机的内参矩阵A。
我们利用R1，R2作为旋转矩阵R的两列，存在单位正交的关系，即：
则由H和R1，R2的关系，可知：
代入可得：
另外，我们发现，上述两个约束方程中均存在矩阵A−TA−1。因此，我们记A-TA-1=B，则B为对称阵。我们视图先求解出矩阵B，通过矩阵B再求解相机的内参矩阵A。
同时，为了简便，我们记相机内参矩阵A为：
则：
则用矩阵A表示矩阵B得：
注意：由于B为对称阵，上式出现了两次B12，B13，B23。
这里，我们可以使用B=A-TA-1将前面通过R1，R2单位正交得到的约束方程化为：
因此，为了求解矩阵B ，我们必须要计算HjTBHj。则：
我们可以记：
则上述方程转化为：
此时，通过R1，R2单位正交得到的约束方程可化为：
即：
由于矩阵H已知，矩阵v又全部由矩阵H的元素构成，因此矩阵v已知。
此时，我们只要求解出向量b，即可得到矩阵B。每张标定板图片可以提供一个ub =0的约束关系，该约束关系含有两个约束方程。但是，向量6有6个未知元素。因此，单张图片提供的两个约束方程是不足以解出来向量b。因此，我们只要取3张标定板照片，得到3个ub=0的约束关系，即6个方程，即可求解向量b。当标定板图片的个数大于3时（事实上一般需要15到20张标定板图片），可采用最小二成拟合最佳的向量6，并得到矩阵B。
根据矩阵B的元素和相机内参a，ß，γ，uo，vo的对应关系（如上式），可得到：
2.4.求解外参矩阵
对于同一个相机，相机的内参矩阵取决于相机的内部参数，无论标定板和相机的位置关系是怎么样的，相机的内参矩阵不变。这也正是在第2部分”求解内参矩阵”中，我们可以利用不同的图片（标定板和相机位置关系不同）获取的矩阵H，共同求解相机内参矩阵A的原因。
但是，外参矩阵反映的是标定板和相机的位置关系。对于不同的图片，标定板和相机的位置关系已经改变，此时每一张图片对应的外参矩阵都是不同的。
在关系：A（R1 R2 T）=H中，我们已经求解得到了矩阵H（对于同一张图片相同，对于不同的图片不同）、矩阵A（对于不同的图片都相同）。通过公式：（R1 R2 T）= A-1H，即可求得每一张图片对应的外参矩阵（R1 R2 T）。
注意，这里值得指出，完整的外参矩阵为
但是，由于张正友标定板将世界坐标系的原点选取在棋盘格上，则棋盘格上任一点的物理坐标W =0，将旋转矩阵的R的第三列R3消掉。因此，R3在坐标转化中并没有作用。但是R3要使得R满足旋转矩阵的性质，即列与列之间单位正交。因此可以通过向量R1，R2的叉乘，即R3 =R1 x R2，计算得到R3。
此时相机的内外参均已经得到
2.5.得到相机畸变矫正参数
张正友标定法仅仅考虑了畸变模型中影响较大的径向畸变。
径向畸变公式（2阶）如下：
其中，（x，y），（â，y）分别为理想的无畸变的归一化的图像坐标、畸变后的归一化图像坐标，r为图像像素点到图像中心点的距离，即r2 = x2 + y2.
图像坐标和像素坐标的转化关系为：
其中，（u，v）为理想的无畸变的像素坐标。由于0接近于90°，则上式近似为：
同理可得畸变后的像素坐标（ü，v）的表达式为：
每一个角点，只要知道畸变后的像素坐标（ü，v）、理想的无畸变的像素坐标（u，v），就可以构造两个上述等式。那么，有m幅图像，每幅图像上有n个标定板角点，则将得到的所有等式组合起来，可以得到mn个未知数为k=[k1，k2]T的约束方程，将约束方程系数矩阵记为D，等式右端非齐次项记为d，可将其记为矩阵形式：
利用最小二乘法：
此时，相机畸变矫正参数已经标定好。
那么，如何获得畸变后的像素坐标（ü，ü）和理想的无畸变的像素坐标（u，v）呢？
（น，v）可以通过识别标定板的角点获得，（u，v）可以通过如下方法近似求得。世界坐标系下的每一个角点的坐标（U，V）是可以计算得到的，我们利用已经求得的外参矩阵（R1 R2 T）和内参矩阵A进行反投影。
利用上式，消去尺度因子Z，可得：
即可得到理想的、无畸变的像素坐标（u，v）。当然，由于外参矩阵（R1 R2 T）和内参矩阵A是在有畸变的情况下获得的，这里得到的像素坐标（u，v）并不是完全理想的、无畸变的。我们的总逻辑是，在进行内参矩阵和外参矩阵的求解的时候，我们假设不存在畸变；在进行畸变系数的求解的时候，我们假设求得的内参矩阵和外参矩阵都是无误差的。最后，我们再通过L-M算法对于参数进行选代优化。
2.6.L-M算法参数优化
从上述推导过程就可以看出，张正友标定法是有很多近似的，所以仅仅利用上述的过程进行标定误差肯定是很大的。所以张正友标定法还利用L-M（Levenberg-Marquardt）算法对参数进行了优化。
L-M算法计算以下函数值最小的参数（k为相机内参矩阵A）
三、实验
3.1 实验要求
1.打印一张棋盘格A4纸张（黑白间距已知），并贴在一个平板上
2.针对棋盘格拍摄若干张图片（一般10-20张）
3.在图片中检测特征点（Harris角点）
4.根据角点位置信息及图像中的坐标，求解Homographic矩阵
5.利用解析解估算方法计算出5个内部参数，以及6个外部参数6.根据极大似然估计策略，设计优化目标并实现参数的refinement
3.2 实验数据环境
华为荣耀V20后置： 4800万像素
拍摄不同角度黑白棋盘格 14张
拍摄相同角度 13张
棋盘格规格为10*7，每个小格边长为0.038m
实验环境：win 10，64位，python3.7，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))
绘画角点的代码：
引入所有函数，绘图板块
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.08590331e+00]
[-7.07566231e-02]
[ 1.73450336e+00]]
[[-1.62797426e-01]
[ 1.50836683e-01]
[ 9.18464323e-04]]
[[-5.34186221e-01]
[ 3.73091917e-03]
[ 1.94074283e-02]]
[[-4.93429591e-01]
[ 5.67761446e-01]
[ 2.95925386e+00]]
[[-5.80071923e-01]
[-7.22436444e-01]
[-1.57453636e+00]]
[[-9.48455655e-01]
[ 3.03024879e-01]
[ 1.29828347e+00]]
[[ 2.94832335e-01]
[ 3.47331204e-02]
[ 2.64544901e-02]]
[[-5.31299069e-02]
[ 3.71927268e-03]
[ 9.28221086e-03]]
[[-6.12810996e-01]
[-4.26476299e-01]
[-1.03092565e+00]]
[[-6.39400599e-01]
[-3.00311309e-01]
[-3.33890353e-01]]
[[-4.36546146e-01]
[ 7.72766764e-01]
[ 1.26422735e+00]]
[[-8.17639859e-01]
[ 3.92249552e-01]
[ 1.46334276e+00]]
[[-5.36634067e-01]
[-1.01048332e+00]
[-2.66320014e+00]]
[[-5.21928699e-01]
[ 1.60164447e-01]
[ 9.02542727e-02]]]
平移向量：
(array([[-1.98022321],
[ 3.90678311],
[21.60898903]]), array([[-4.85885981],
[-0.6289291 ],
[22.06139288]]), array([[-3.06791902],
[-1.41980029],
[17.73790133]]), array([[ 5.44462226],
[ 3.05527622],
[18.99937008]]), array([[-7.10288561],
[ 7.38181483],
[16.07917721]]), array([[-4.38042851],
[-0.1470518 ],
[23.31177057]]), array([[-3.28338334],
[-0.41788501],
[18.69628827]]), array([[-3.32515074],
[-1.11411687],
[18.98961375]]), array([[-3.51751419],
[ 2.3812171 ],
[20.44982348]]), array([[-2.20115011],
[-0.66322809],
[19.38790959]]), array([[-3.47406785],
[ 1.6223005 ],
[21.05495619]]), array([[-2.90648048],
[ 1.91986687],
[21.78508256]]), array([[-3.97603321],
[ 8.24589791],
[19.4315128 ]]), array([[-4.45135664],
[-1.37711117],
[20.59510491]]))
- 相同角度拍摄的角点
内参矩阵：
畸变矩阵：
旋转变量：
[[[ 0.24891805]
[-0.26431711]
[ 1.56918475]]
[[ 0.25431227]
[-0.26637111]
[ 1.5643228 ]]
[[ 0.25283782]
[-0.2733382 ]
[ 1.5657245 ]]
[[ 0.25292215]
[-0.27137721]
[ 1.57478432]]
[[ 0.25054907]
[-0.2749431 ]
[ 1.57254583]]
[[ 0.24767626]
[-0.27061521]
[ 1.57629155]]
[[ 0.27330193]
[-0.27878684]
[ 1.4781248 ]]
[[ 0.24911727]
[-0.27417874]
[ 1.57415463]]
[[ 0.25097148]
[-0.27088346]
[ 1.57368064]]
[[ 0.24943128]
[-0.27150587]
[ 1.58753057]]
[[ 0.24931312]
[-0.25788696]
[ 1.56927848]]
[[ 0.25003405]
[-0.26845351]
[ 1.57669215]]
[[ 0.25437311]
[-0.26630249]
[ 1.57098008]]]
平移变量：
(array([[ 2.55270454],
[ -5.36701444],
[144.98007943]]), array([[ 2.44283884],
[ -5.42338893],
[144.5119281 ]]), array([[ 2.49552813],
[ -5.37504836],
[144.18066185]]), array([[ 2.67937867],
[ -5.37984686],
[144.3546585 ]]), array([[ 2.62616415],
[ -5.41405278],
[144.0528159 ]]), array([[ 2.63307209],
[ -5.44047493],
[143.97896365]]), array([[ 1.58925241],
[ -5.56047437],
[139.15058038]]), array([[ 2.66180544],
[ -5.38680787],
[144.15643877]]), array([[ 2.66953543],
[ -5.40095104],
[143.91674766]]), array([[ 2.79989396],
[ -5.44070812],
[144.26010846]]), array([[ 2.52335486],
[ -5.32486579],
[145.17592865]]), array([[ 2.70933161],
[ -5.44882022],
[143.75022183]]), array([[ 2.55187977],
[ -5.34339341],
[144.54043873]]))
四、实验总结
- 特征点提取误差对标定精度的影响因为相机标定是在已知特征点的世界坐标系坐标和特征点的像面坐标的情况下求解内外参的，若是特征点提取存在误差的话势必会影响标定的精度，研究表明，标定误差与特征点提取误差的均值成正比。
- 标定图片数量对标定精度的影响
和标定点的数量一样，标定图片的数量也对标定精度有影响，实验表明，标定图片时，误差较大。随着标定图片增多，误差逐渐减小。超过一定量后，误差大小变化变小。 - 标定的效果与棋盘格标准以及表面是否平整，表面褶皱也有关系
- 同一相机用同一组数据 得到的内参相同，畸变系数不变，同一相机不同数据 相机内参不同因为相机的小孔成像并非标准。只是近似模型，不同距离相机光心不同，内参不同
Original: https://blog.csdn.net/weixin_46462313/article/details/124917748
Author: 小韩丸子
Title: 计算机视觉—— 相机标定
原创文章受到原创版权保护。转载请注明出处：https://www.johngo689.com/648067/
转载文章受原作者版权保护。转载请注明原作者出处！