计算机视觉--python+opencv相机标定

前言

张正友标定法简介

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

一、相机标定原理

1.1 相机模型

1.1.1 各个坐标系

    确定空间某点的三维几何位置与其在图像中对应点之间的相互关系，必须建立相机成像的几何模型（各个坐标系），这些坐标系之间的转换参数就是相机参数，求解参数的过程叫做相机标定（摄像机标定）。建立立体视觉系统所需要的各个坐标系，包括世界坐标系、相机坐标系、以及图像坐标系（物理和像素坐标系）。

坐标系名称	坐标系描述
世界坐标系（3D）	描述目标在真实世界中的位置引入的参考坐标系(Xw，Yw，Zw)
相机坐标系（3D）	联系世界坐标系与图像坐标系的桥梁，一般取摄像机的光学轴为ｚ轴(Xc，Yc，Zc)
图像物理坐标系（2D）	根据投影关系引入，方便进一步得到像素坐标，单位为毫米，坐标原点为摄像机光轴与图像物理坐标系的交点位置(x，y)
图像像素坐标系（2D）	真正从相机内读到的信息，图像物理坐标的离散化，以像素为单位，坐标原点在左上角(u,v)

1.1.2 相机畸变模型

一般只考虑径向畸变k和切向畸变p
$$x_{distorted}=x\ast (1+k_1{r}^2+k_2{r}^4+k_3{r}^6)+2P_{1}xy+p_2({\gamma }^2+2x^2)$$$$y_{distorted}=y\ast (1+k_1{r}^2+k_2{r}^4+k_3{r}^6)+2P_{1}xy+p_2({\gamma }^2+2y^2)$$
其中${\gamma }^2=x^2+y^2$ 。
一般选择2或3个k值，这个是经验证过可以获得较好的结果，如果k取得再多影响不会很大可以忽略，甚至可能反而导致效果不好。
畸变模型：枕型畸变(k>0)和桶型畸变(k<0)

当k>0时，r越大（点离中心越远），畸变量越大，r越小，畸变量越小，呈枕型。
当k<0时，r越大（点离中心越远），畸变量越小，r越小，畸变量越大，呈桶型。

1.1.3 相机标定参数

内参：
（单位长度的像素个数）$f_x、f_y$
（主点坐标）$C_x、C_y$
（畸变系数）$k_1、k_2、k_3、p_1、p_2$
外参：R、T（旋转和平移矩阵）

1.2 张正友标定法

1.2.1 计算单应性矩阵H

单应性：在计算机视觉中被定义为一个平面到另一个平面的投影映射。首先确定，图像平面与标定物棋盘格平面的单应性。

• 不妨设棋盘格位于Z = 0
• 定义旋转矩阵R的第i列为r_i 则有：
  $$s\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]=K\left[\begin{array}{cccc}{r}_1& r_2& r_3& t\end{array}\right]\left[\begin{array}{c}X\\ Y\\ 0\\ 1\end{array}\right]=K\left[\begin{array}{ccc}{r}_1& r_2& t\end{array}\right]\left[\begin{array}{c}X\\ Y\\ 1\end{array}\right]$$
• 于是空间到图像的映射可改为:
  $$s\tilde{m}=H\tilde{M},H=K\left[\begin{array}{ccc}{r}_1& r_2& t\end{array}\right]$$
• 其中H是描述Homographic矩阵,可通过最小二乘 ，从角点世界坐标到图像坐标的关系求解。

1.2.2 计算内参数矩阵

• 令 H 为$H=\left[\begin{array}{ccc}{h}_1& h_2& h_3\end{array}\right]$，则$\left[\begin{array}{ccc}{h}_1& h_2& h_3\end{array}\right]=\lambda K\left[\begin{array}{ccc}{r}_1& r_2& t\end{array}\right]$
• Homography 有 8 个自由度
• 通过上述等式的矩阵运算，根据r1和r2正交，以及 归一化的约束可以得到如下等式：
  
• 即每个单应性矩阵能提供两个方程，而内参数矩阵包含5个参数，要求解，至少需要3个单应性矩阵。为了得到三个不同的单应性矩阵，我们使用至少三幅棋盘格平面的图片进行标定。通过改变相机与标定板之间的相对位置来得到三个不同的图片。为了方便计算，我们定义：
  
• B 是对称阵，其未知量可表示为6D 向量 b
  
• 设H中的第i列为 h_i
• 根据b的定义，可以推导出如下公式
  
• 最后推导出
  
  通过上式，我们可知当观测平面 n ≥ 3 时，即至少3幅棋盘格图像，可以得到b的唯一解，求得相机内参数矩阵K。

1.2.3 计算外参数矩阵

外部参数可通过Homography求解，由 H = [h1 h2 h3] = λK[r1 r2 t]，可推出

1.2.4 最大似然估计

上述的推导结果是基于理想情况下而言，但由于可能存在一些其他干扰，所以使用最大似然估计进行优化。假设拍摄了n张棋盘格图像，每张图像有m个角点。最终获得的最大似然估计公式为

二、棋盘格标定实验

2.1 实验步骤

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

2.2 数据集及环境介绍

数据集介绍：

• 打印棋盘格并贴在一个平板上，针对棋盘格拍摄若干张图片（共拍摄15张）
• 棋盘格的规格为9（行）*7（列），每小方格边长28mm.
  

环境介绍：
win10，64位，python3.7，opencv4.2.0，手机型号：华为荣耀9（STF-AL10）

2.3 代码实现

# -*- coding: utf-8 -*-
import os
import numpy as np
import cv2
import glob
import re

def calib(inter_corner_shape, size_per_grid, img_dir, img_type):
    # criteria: only for subpix calibration, which is not used here.
    # criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 30, 0.001)
    w, h = inter_corner_shape
    # cp_int: corner point in int form, save the coordinate of corner points in world sapce in 'int' form
    # like (0,0,0), (1,0,0), (2,0,0) ....,(10,7,0).
    cp_int = np.zeros((w * h, 3), np.float32)
    cp_int[:, :2] = np.mgrid[0:w, 0:h].T.reshape(-1, 2)
    # cp_world: corner point in world space, save the coordinate of corner points in world space.
    cp_world = cp_int * size_per_grid

    obj_points = []  # the points in world space
    img_points = []  # the points in image space (relevant to obj_points)
    images = glob.glob(img_dir + os.sep + '**.' + img_type)
    for fname in images:
        # print(re.sub("\D", "", fname))
        img = cv2.imread(fname)
        gray_img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
        # find the corners, cp_img: corner points in pixel space.
        ret, cp_img = cv2.findChessboardCorners(gray_img, (w, h), None)
        # if ret is True, save.
        if ret == True:
            # cv2.cornerSubPix(gray_img,cp_img,(11,11),(-1,-1),criteria)
            obj_points.append(cp_world)
            img_points.append(cp_img)
            # view the corners
            cv2.drawChessboardCorners(img, (w, h), cp_img, ret)
            cv2.imwrite('D:\\Alike\\python_dw\\Code_a\\admin_code\\data\save\\p' + re.sub("\D", "", fname) + '.jpg', img)
            # cv2.waitKey(1000)
    # cv2.destroyAllWindows()
    # calibrate the camera
    ret, mat_inter, coff_dis, v_rot, v_trans = cv2.calibrateCamera(obj_points, img_points, gray_img.shape[::-1], None,
                                                                   None)
    print(("ret:"), ret)
    print(("internal matrix:\n"), mat_inter)
    # in the form of (k_1,k_2,p_1,p_2,k_3)
    print(("distortion cofficients:\n"), coff_dis)
    print(("rotation vectors:\n"), v_rot)
    print(("translation vectors:\n"), v_trans)
    # calculate the error of reproject
    total_error = 0
    for i in range(len(obj_points)):
        img_points_repro, _ = cv2.projectPoints(obj_points[i], v_rot[i], v_trans[i], mat_inter, coff_dis)
        error = cv2.norm(img_points[i], img_points_repro, cv2.NORM_L2) / len(img_points_repro)
        total_error += error
    print(("Average Error of Reproject: "), total_error / len(obj_points))

    return mat_inter, coff_dis


def dedistortion(inter_corner_shape, img_dir, img_type, save_dir, mat_inter, coff_dis):
    w, h = inter_corner_shape
    images = glob.glob(img_dir + os.sep + '**.' + img_type)
    for fname in images:
        img_name = fname.split(os.sep)[-1]
        img = cv2.imread(fname)
        newcameramtx, roi = cv2.getOptimalNewCameraMatrix(mat_inter, coff_dis, (w, h), 0, (w, h))  # 自由比例参数
        dst = cv2.undistort(img, mat_inter, coff_dis, None, newcameramtx)
        # clip the image
        # x,y,w,h = roi
        # dst = dst[y:y+h, x:x+w]
        cv2.imwrite(save_dir + os.sep + img_name, dst)
    print('Dedistorted images have been saved to: %s successfully.' % save_dir)


if __name__ == '__main__':
    inter_corner_shape = (8, 6)
    size_per_grid = 0.028

    img_dir = "D:\\Alike\\python_dw\\Code_a\\admin_code\\data\data"
    img_type = "jpg"
    # calibrate the camera
    mat_inter, coff_dis = calib(inter_corner_shape, size_per_grid, img_dir, img_type)
    # dedistort and save the dedistortion result.
    save_dir = "D:\\Alike\\python_dw\\Code_a\\admin_code\\data\pic\\save_dedistortion"
    if (not os.path.exists(save_dir)):
        os.makedirs(save_dir)
    dedistortion(inter_corner_shape, img_dir, img_type, save_dir, mat_inter, coff_dis)

2.4 实验结果

2.4.1 角点检测结果

角点检测结果：
单张——

所有图像检测——

2.4.2 内部参数计算结果

内部参数：
以下为控制台输出的内部参数部分，分别表示：内参数矩阵、畸变系数。

2.4.3 外部参数计算结果

外部参数：
以下为控制台输出的外部参数部分，分别表示：旋转向量、平移向量。

2.4.4 外部参数可视化结果

外部参数可视化
①.以棋盘格为中心

②.以相机为中心

2.4.5 平均重投影误差

以下为控制台输出的平均重投影误差，大约为0.3744.

2.5 实验分析

• 在准备实验数据时（即棋盘格图像），最好将打印出的棋盘格A4纸贴在一块平板上再进行拍摄。起初我没有将其固定于平板上，由于A4纸本身有褶皱，会导致实验结果不准确。并且拍摄时最好选择光线稍微暗的地方，不要出现反光现象。
• 从实验结果可以得出我手机的相机参数，internal matrix是内参数矩阵。归一化焦距fx=320.90727491，fy=320.52453665，像主点（光心）的坐标Cx=143.6485789，Cy=191.96528843。
• 见以上平均重投影误差结果，发现为0.3744，该结果可以用于评估标定效果的好坏。为了更方便观察，除了使用上述的python语言完成实验，我又使用matlab对图片集进行标定实验，得到结果如下：
  
  可以看到，使用matlab进行实验，得到的平均重投影误差为0.40，与先前的实验结果大致相同，柱形图有助于观察每幅图像的误差，可以看到编号为2的图像误差最大，编号1的误差最小，将其调出，如下图：
  
      比较图1和图2，可以看到图2存在径向畸变现象，产生径向畸变的主要原因是由于透镜的制造原理使光线在远离透镜中心的地方折射效果更加明显，越远离中心点畸变的效果越显著。由于拍摄角度原因，看到图2中对于空间中一条直线，由于径向畸变在成像平面上会变成一条曲线，整体误差大。而图片1的棋盘格几乎处在图像正中间，且保持平坦，所以重投影误差很小。
      就总体来说，平均重投影误差小于0.5，即可认为标定效果较好。

三、实验小结

注意点：
实验时需要关注自己打印出来的棋盘格行列数以及每个小方格的边长，代码中与之相关的参数要注意保持一致。

遇到问题：
①opencv输出图片时会一闪而过。
解决方案：
经检查，发现代码中原本写的是:

cv2.waitKey(1)

将其改为

cv2.waitKey(0)

可以实现图片输出。
或者也可以将图片写入文件夹，如下：

 cv2.imwrite('D:\\Alike\\python_dw\\Code_a\\admin_code\\data\save\\p' + re.sub("\D", "", fname) + '.jpg', img)

②用matlab进行标定实验时报错：

解决方案：
将所有图片像素值改成一致，可以适当调小像素，运行速度会更快。