用Python进行相机校准

汉无为 2022-05-02 发布于湖北

展开全文

介绍

相机校准的基本思想是给定世界上已知的一组点及其在图像中的相应投影，我们必须找到负责投影变换的矩阵。然后我们可以使用这个矩阵将世界上的任何点投影到图像上。我们将在本文中了解如何执行相机校准。

在此之前，让我们回顾一下。

回顾

相机外参矩阵

在本系列的第2部分中，我们看到相机外参矩阵是基础矩阵的一种变化，它将一个点的坐标从世界坐标系转换为相机坐标系。

https:///camera-extrinsic-matrix-with-example-in-python-cfe80acab8dd

它让我们从相机的角度来观察世界。我们还看到，它是旋转矩阵和平移矩阵的组合——旋转矩阵定向相机，平移矩阵移动相机。

因此，给定世界上一个点的坐标，我们可以应用相机外参矩阵来改变它在相机上的坐标，如下所示：

此处的符号如下所示：

在这里，这些点被表示为齐次坐标，这实质上是在原始坐标上增加一个额外的维度。

相机外参矩阵的最后一行只是0和1，它们不会给变换增加任何值，因此我们可以安全地删除最后一行并重写如下等式：

注意这里输出的形状是3×1，而之前是4×1。这意味着点在这里用欧几里德形式表示，这是一件好事，因为我们不需要额外的步骤来从齐次坐标转换回来。

此外，摄像机的内参矩阵，即管道中的下一个矩阵，接受3×1的输入形状，我们将在下面看到。

总之，给定世界坐标系中某个点的坐标，相机外参矩阵将坐标转换为相机坐标系。这种转变可以表示为：

摄像机内参矩阵

现在，我们已经使用相机外参矩阵得到了相机上的一个点的坐标，下一步是将该点投影到相机的图像平面上，并形成一个图像。这是相机内参矩阵的工作。

本系列的第3部分已经深入讨论了相机内参矩阵：

https:///camera-intrinsic-matrix-with-example-in-python-d79bf2478c12

但总而言之，相机内参矩阵将相机给定坐标的点投影到相机的图像平面上。从本质上说，我们可以使用相机内参矩阵来获得图像中点的像素位置。如下所示：

这些符号是：

输出将是同质形式的投影。要将其转换为欧几里德形式，我们只需除以最后一个坐标，并丢弃最后一个维度，如下所示：

这里（u，v）是图像中点P的像素位置。摄像机内参矩阵由以下表达式表示：𝜅, 这种转变可以表示为：

结合相机的外参和内参矩阵

结合相机的外在和内参矩阵，我们可以建立一个管道，以世界坐标系中的一个点为输入，计算其在摄像机形成的图像上的投影。它可以表示为：

此外，使用矩阵组合，我们可以将两个矩阵组合成一个矩阵M，如下所示：

形状为3×4的矩阵M具有所有必需的信息。

整合

使用矩阵M，我们可以将整个管道表示为：

如果你观察，矩阵M有12个元素，但我们在之前的文章中看到，内参矩阵和外参矩阵有11个自由度（6个来自外参矩阵，5个来自内参矩阵）。所以M中有一个元素依赖于其他元素。

假设最后一个元素M（3，4）是依赖元素。这意味着如果我们用这个元素划分矩阵M，M中的信息不会受到影响，因为它是一个依赖元素。

在上面的等式中，我们将矩阵M与最后一个元素M（3，4）分开，得到一个新的矩阵。这里需要注意的重要一点是，尽管进行了操作，但两种情况下的输出或投影都是相同的。矩阵M的这个性质被称为尺度不变性，这意味着我们可以用任何因子对矩阵M进行缩放，而不会影响输出。

通常，在实际应用中，矩阵M是未知的，相机校准的目标是使用一组已知点及其投影来找到它。

几何相机校准

我们可以从上述矩阵乘法中找到图像坐标（ui，vi），如下所示：

我们可以将这个等式改写为：

我们可以进一步简化为：

如果你观察，上面的方程看起来像一个齐次线性方程，以 m⃗ 为单位，通常有如下形式：

现在，齐次线性方程组可以用矩阵形式表示为A𝑥⃗ = 0 ⃗，其中A是m×n矩阵，𝑥⃗ 是包含n个条目的列向量，0 ⃗ 是包含m个条目的零向量。

现在我们可以用矩阵形式来表示我们的方程Am⃗ = 0 ⃗

上述方程式代表Am⃗ = 0 ⃗。m⃗ 是向量形式的平坦矩阵M。记住，我们的目标是求矩阵M的系数，这与求解这个齐次系统是一样的。

现在，每个齐次系统都至少有一个解，称为零解，它是通过给每个元素赋值得到的。但这不是我们正在寻找的解。那么，我们怎样才能找到这个齐次系统的非零解呢？

我们可以找到一个近似解，而不是找到一个精确解，我们甚至不确定它是否存在。

从数学上讲，这意味着我们不能找到Am⃗*=0* ⃗的精确解, 但是我们可以找到m⃗ ，以至于|Am⃗| 最小。本质上，我们试图最小化代数误差。

此外，我们可以将m⃗标准化，所以它变成了单位向量。

所以，现在我们的问题是求*|Am|*的最小值，且以|m⃗|=1。

在本系列的第4部分中（https:///find-the-minimum-stretching-direction-of-positive-definite-matrices-79c2a3b397fc），我们知道了如果单位向量𝑥⃗ 沿着𝐴⊺𝐴的最小特征向量的方向，那么|𝐴𝑥⃗| 能被最小化

让我们先看看矩阵A：

我们的想法是，我们将找到一些由（Xi，Yi，Zi）表示的点，并在由（ui，vi）表示的图像中找到它们相应的投影，然后我们可以计算矩阵A。

这个标记过程是手动完成的。现在，如果我们标记n个点，矩阵A的形状将是2n×12，而m⃗ 的形状将保持固定在12×1。零向量的形状是2n×1。

这里有一个问题要问你：我们需要标注多少点？m⃗的大小是12，这意味着有12个未知数，所以我们需要解12个独立的方程来找到这12个未知数。所以我们要求n至少为6，这意味着我们应该标记至少6个独立点来求解m⃗ 。我们还可以标记更多的点，一般来说，点越多越好；但6是最低要求。

整合

好的，计划是标记至少6个点，然后计算矩阵A，然后我们可以计算𝐴⊺𝐴 的特征向量还有m⃗ 。最后，我们可以重塑向量m⃗ 得到3×4矩阵M，这将是我们校准的摄像机矩阵。使用M，我们可以找到图像中任何世界点的投影。

让我们用一个例子来说明这一点。

实例

设置

包含所有代码的GitHub存储库可以在这里找到。

https://github.com/wingedrasengan927/Image-formation-and-camera-calibration

假设之前没有设置环境，现在可以通过运行以下命令来完成：

# create a virtual environment in anaconda
conda create -n camera-calibration-python python=3.6 anaconda
conda activate camera-calibration-python

# clone the repository and install dependencies
git clone https://github.com/wingedrasengan927/Image-formation-and-camera-calibration.git
cd Image-formation-and-camera-calibration
pip install -r requirements.txt

注意：这假设你已经安装了anaconda。

我们将使用两个主要的库：

pytransform3d：这个库对于三维空间中的可视化和转换具有强大的功能。
ipympl：它使matplotlib绘图具有交互性，允许我们在笔记本电脑中实时执行平移、缩放和旋转，这在使用3D绘图时非常有用。

例子

为了进行摄像机校准，我们首先需要准备真实标签，它本质上是世界上的一组点及其在图像上的相应投影。

在现实世界中，我们手动测量点与相机的距离，并在图像中找到相应的像素。然而，在计算机上，我们可以模拟这个过程——我们可以创建一个相机外参矩阵和内参矩阵，并构建一个管道来计算世界点的投影。

我们已经在本系列的第2部分和第3部分中看到了如何创建矩阵。一旦我们得到了真实标签，我们就可以用它来构造代数矩阵A，然后找到它的转置并进行计算𝐴⊺𝐴, 最后，我们可以取𝐴⊺𝐴 用最小特征值对其进行重塑，得到矩阵M。

以下是代码实现：

%matplotlib widget

import matplotlib.pyplot as plt
from utils import *
np.random.seed(42)

创建矩阵

定义参数

# define extrinsic parameters
# -------------------------------

# rotate an angle of pi/4 along the standard Y axis
angles = [np.pi/4]
order = 'y'

# transalte by the given offset
offset = np.array([0, -8, 0])

# define intrinsic parameters
# -------------------------------

f = 2
s = 0
a = 1
cx = 0
cy = 0
img_size = (10, 10)

创建外参和内参矩阵

# create extrinsic matrix
# --------------------------

# create rotation transformation matrix
R = create_rotation_transformation_matrix(angles, order)
R_ = np.identity(4)
R_[:3, :3] = R

# create translation transformation matrix
T_ = create_translation_matrix(offset)

E = np.linalg.inv(R_ @ T_)
E = E[:-1, :]

# create intrinsic matrix
# ---------------------------
K = compute_intrinsic_parameter_matrix(f, s, a, cx, cy)

生成随机点

# choose the lower limit of the points such they're always beyond the image plane

n_points = 12
rand_points = generate_random_points(n_points, (-10, 0), (-10, 10), (f, 10))

绘图设置

# create an image grid
xx, yy, Z = create_image_grid(f, img_size)

# convert the image grid to homogeneous coordinates
pt_h = convert_grid_to_homogeneous(xx, yy, Z, img_size)

# transform the homogeneous coordinates
pt_h_transformed = R_ @ T_ @ pt_h

# convert the transformed homogeneous coordinates back to the image grid
xxt, yyt, Zt = convert_homogeneous_to_grid(pt_h_transformed, img_size)

# define axis and figure
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111,projection='3d')

# set limits
ax.set(xlim=(-10, 5), ylim=(-15, 5), zlim=(0, 10))

# plot the camera in the world
ax = pr.plot_basis(ax, R, offset)
ax.plot_surface(xxt, yyt, Zt, alpha=0.75)

# plot the generated random points
c = 0
for i in range(n_points):
    point = rand_points[:, c]
    ax.scatter(*point, color='orange')
    ax.plot(*make_line(offset, point), color='purple', alpha=0.25)
    c += 1

ax.set_title('The Setup')
ax.set_xlabel('X-axis')
ax.set_ylabel('Y-axis')
ax.set_zlabel('Z-axis')

计算这些点的投影并形成图像

rand_points_camera = compute_coordniates_wrt_camera(rand_points, E, is_homogeneous=False)

projections = compute_image_projection(rand_points_camera, K)

fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111)

for i in range(n_points):
ax.scatter(*projections.reshape(-1, 2)[i], color='orange')

ax.set_title('projection of points in the image')

直接线性校正

建立代数矩阵A并求出m

# compute the algebraic matrix A
A = create_algebraic_matrix(rand_points, projections)

# compute At x A
A_ = np.matmul(A.T, A)
# compute its eigenvectors and eigenvalues
eigenvalues, eigenvectors = np.linalg.eig(A_)
# find the eigenvector with the minimum eigenvalue
# (numpy already returns sorted eigenvectors wrt their eigenvalues)
m = eigenvectors[:, 11]

# reshape m back to a matrix
M = m.reshape(3, 4)

从校准的矩阵M计算预测

predictions = compute_world2img_projection(rand_points, M, is_homogeneous=False)

画出预测和真实值

fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111)

for i in range(n_points):
    if i == 0:
        o_label = 'groundtruth'
        g_label = 'predictions'
    else:
        o_label = ''
        g_label = ''
        
    ax.scatter(*projections.reshape(-1, 2)[i], color='orange', alpha=0.75, label=o_label)
    ax.scatter(*predictions.reshape(-1, 2)[i], color='green', alpha=0.75, label=g_label)
    
ax.set_title('groundtruth vs predictions - direct linear calibration')
ax.legend()

优化wrt几何误差

from scipy.optimize import minimize

result = minimize(geometric_error, m, args=(rand_points, projections))

M_ = result.x.reshape(3, 4)

predictions_v2 = compute_world2img_projection(rand_points, M_, is_homogeneous=False)

fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111)

for i in range(n_points):
    if i == 0:
        o_label = 'groundtruth'
        g_label = 'predictions'
    else:
        o_label = ''
        g_label = ''

    ax.scatter(*projections.reshape(-1, 2)[i], color='orange', alpha=0.5, label=o_label)
    ax.scatter(*predictions_v2.reshape(-1, 2)[i], color='green', alpha=0.5, label=g_label)

ax.set_title('groundtruth vs predictions - optimization wrt geometric error')
ax.legend()

fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(8, 6))

for i in range(n_points):
        
    axes[0].scatter(*projections.reshape(-1, 2)[i], color='orange', label=o_label)
    axes[1].scatter(*predictions_v2.reshape(-1, 2)[i], color='green', label=g_label)
    
axes[0].set_title('groundtruth')
axes[1].set_title('predictions')

plt.tight_layout()

让我们一步一步地浏览，了解发生了什么：

首先，我们定义必要的参数，并创建相机外参矩阵和内参矩阵。这些都是建造管道和准备标签必需的。然后我们在世界上生成n个随机点（这里我选择了n=12）。这是一个用相机和随机点绘制的空间图。

我们可以计算这些点的投影并形成一个图像。首先，我们必须应用外参矩阵来表示摄像机坐标系中的点，然后我们可以应用内参矩阵来获得投影，最后，我们可以在图像中绘制这些投影，如下所示。

现在我们已经准备好了真实值，就是世界上的点（Xi，Yi，Zi）以及它们在图像中的对应投影（ui，vi）。使用它们我们可以计算代数矩阵A。如果你记得的话，这就是矩阵A：

在计算矩阵A之后，我们可以计算𝐴⊺𝐴 并用最小的特征值选取其特征向量。这将给我们大小为12的向量m⃗。最后，我们可以重塑m⃗为 3×4得到我们要找的校准矩阵M。
我们可以直接使用矩阵M来计算世界点在图像上的投影。它有所有必要的信息。让我们将矩阵M应用于我们生成的n个随机点；通过这种方式，我们可以将校准矩阵M的投影与真实标签进行比较。下图显示了比较结果。

嗯……我们校准矩阵M的预测并不十分准确。有些点接近实际情况，有些则相去甚远，但总的来说并不好。在下一节中，我们将讨论可能发生了什么，以及我们可以做些什么来改进。

几何误差

我们上面讨论的方法叫做直接线性校准。它不是很好的一个原因是因为我们试图最小化代数误差——我们基本上是在试图找到m⃗，它很好地符合代数方程Am⃗ = 0 ⃗。

这种线性优化的问题是，当非线性误差和不确定性（如径向畸变）潜入相机时（这在现实世界中经常发生），算法将严重失败。我们应该针对正确的错误类型进行优化。

我们应该看到的误差叫做几何误差。几何误差基本上给出了预测与真实标签之间的距离估计。它通过测量点的预测投影与其真实情况的投影之间的距离来实现。

当我们最小化几何误差时，我们实际上是在最小化预测和真实标签之间的距离。

几何误差由以下公式给出：

这里d是距离度量，我们通常使用欧几里德距离。在我们的例子中，预测𝑥′𝑖 由MXi给出，其中M是我们标定的摄像机矩阵，Xi是世界坐标系中的一个点。我们可以把方程改写为：

好的，我们的想法是执行某种非线性优化，并更新矩阵M的权重，以最小化几何误差。

这种方法在直觉上类似于机器学习，我们使用梯度下降算法更新模型以最小化损失。幸运的是，我们不需要自己编写优化算法。scipy在scipy中提供了几十种优化算法。优化模块。

让我们看看如何在代码中实现：

from scipy.optimize import minimize

def geometric_error(m, world_points, projections):
    '''
    compute the geometric error wrt the
    prediction projections and the groundtruth projections

    Parameters
    ------------
    m - np.ndarray, shape - (12)
        an 12-dim vector which is to be updated

    world_points - np.ndarray, shape - (3, n)
                   points in the world coordinate system

    projections - np.ndarray(2, n)
                  projections of the points in the image

    Returns
    --------
    error - float
            the geometric error
    '''

    assert world_points.shape[1] == projections.shape[1]

    error = 0
    n_points = world_points.shape[1]

    for i in range(n_points):

        X, Y, Z = world_points[:, i]
        u, v = projections[:, i]

        u_ = m[0] * X + m[1] * Y + m[2] * Z + m[3]
        v_ = m[4] * X + m[5] * Y + m[6] * Z + m[7]
        d = m[8] * X + m[9] * Y + m[10] * Z + m[11]

        u_ = u_/d
        v_ = v_/d

        error += np.sqrt(np.square(u - u_) + np.square(v - v_))

    return error

result = minimize(geometric_error, m, args=(rand_points, projections))
M_ = result.x.reshape(3, 4)

predictions_v2 = compute_world2img_projection(rand_points, M_, is_homogeneous=False)

在上面的代码中，我们首先定义几何误差函数。此函数接受向量m⃗作为参数以及真实标签，这是一组世界点及其在图像中的相应投影。它使用我们上面讨论的公式计算预测：

最后，我们计算所有点的预测投影和实际投影之间的欧氏距离，并计算几何误差。

从模块scipy.optimize导入的函数minimize包含两个重要参数——误差函数和初始权重。我们可以将上面定义的几何误差函数作为第一个参数传递，作为第二个参数，我们可以将任何12维向量作为初始状态传递——然而，因为我们已经从直接线性校准方法计算了m⃗，我们可以通过它作为初始状态，而不是一些随机向量。作为第三个参数，我们可以传递一个包含error函数参数的元组。

执行最小化函数后，它运行并完成优化过程，并返回一个包含结果的对象。我们可以通过x属性访问该对象的最终更新权重。然后，我们可以将这个12维的权重向量重塑为一个3×4矩阵，并计算预测。下图是我们从该方法得到的预测与真实投影的比较：