相机校准—外参矩阵

点击上方“小白学视觉”，选择加"星标"或“置顶”

重磅干货，第一时间送达

在本文中，我们将探讨摄影机的外参，并通过Python中的一个实践示例来加强我们的理解。

相机外参

摄像头可以位于世界任何地方，并且可以指向任何方向。我们想从摄像机的角度来观察世界上的物体，这种从世界坐标系到摄像机坐标系的转换被称为摄像机外参。

那么，我们怎样才能找到相机外参呢？一旦我们弄清楚相机是如何变换的，我们就可以找到从世界坐标系到相机坐标系的基变换的变化。我们将详细探讨这个想法。

具体来说，我们需要知道相机是如何定位的，以及它在世界空间中的位置，有两种转换可以帮助我们：

• 有助于确定摄影机方向的旋转变换。

• 有助于移动相机的平移变换。

让我们详细看看每一个。

旋转

通过旋转改变坐标

让我们看一下将点旋转一个角度的变换。让我们举一个在ℝ²的简单例子，对点𝑃逆时针旋转角度𝜃 得到点𝑃′， 如下图所示：

𝑃的坐标是(𝑥,𝑦) 以及𝑃′的坐标是(𝑥′,𝑦′). 我们需要找到(𝑥′,𝑦′).

从图来看，
      sinα = y/r , cosα = x/r − [1]
⟹    xsinα = ycosα − [2]
同样的, x′ = rcos(θ+α)
⟹    x′ = (x/cosα) ∗ cos(θ+α) (from [1])
但,  cos(θ+α) = cosθcosα − sinθsinα
⟹    x′ = (x/cosα) ∗ (cosθcosα − sinθsinα)
⟹    x′ = xcosθ − xsinα ∗ (sinθ / cosα)
⟹    x′ = xcosθ − ycosα ∗ (sinθ / cosα) (from [2])
⟹    x′ = xcosθ − ysinθ

同样地，
      y′ = rsin(θ+α)
⟹    y′ = (y/sinα) ∗ sin(θ+α) (from [1])
但,  sin(θ+α) = sinθcosα + cosθsinα
⟹    y′ = (y/sinα) ∗ (sinθcosα + cosθsinα)
⟹    y′ = ycosθ + ycosα ∗ (sinθ / sinα)
⟹    y′ = ycosθ + xsinα ∗ (sinθ / sinα) (from [2])
⟹    y′ = ycosθ + xsinθ
⟹    y′ = xsinθ + ycosθ

因此我们有，
      x′ = xcosθ − ysinθ
      y′ = xsinθ + ycosθ

旋转是一种线性运算，上述方程可以表示为矩阵乘法：

这个操作是一个线性变换。

扩展到R3

我们可以很容易地将旋转变换扩展到𝐑³. 旋转的变换矩阵𝐑关于标准X轴、Y轴和Z轴，如下所示：

绕Z轴旋转：

绕X轴旋转：

绕Y轴旋转：

内参旋转与外参旋转

上述变换围绕标准轴执行旋转。轴将随时固定。这就是所谓的外参旋转。还有另一种类型的旋转称为内参旋转，我们在每一步都围绕其相对轴旋转对象，如下所示：

内参旋转很难用欧几里德代数来实现，我们将坚持外参旋转。

基变换

在基变换中，我们的目标是在新的基上找到点的坐标。

在下面的示例中𝑋𝑌 轴已经旋转了一个角度𝜃 得到𝑋′𝑌′. 给定在原有XY轴下点𝑃的坐标 , 我们的目标是找到在新轴𝑋′𝑌′下点𝑃的坐标 .

XY下点P的坐标是(x,y) ，新轴X'Y'下是(𝑥′, 𝑦′). 我们的目标是找到(𝑥′, 𝑦′).

从这个图来看，
      sinα = y′/r , cosα = x′/r − [1]
⟹    x′sinα = y′cosα − [2]
同样, x = rcos(θ+α)
⟹    x = (x′/cosα) ∗ cos(θ+α) (from [1])
但,  cos(θ+α) = cosθcosα − sinθsinα
⟹    x = (x′ / cosα) ∗ (cosθcosα − sinθsinα)
⟹    x = x′cosθ − xsinα ∗ (sinθ / cosα)
⟹    x = x′cosθ − y′cosα ∗ (sinθ / cosα) (from [2])
⟹    x = x′cosθ − y′sinθ

同样地，
      y = rsin(θ+α)
⟹    y = (y′/sinα) ∗ sin(θ+α) (from [1])
但,  sin(θ+α) = sinθcosα + cosθsinα
⟹    y = (y′/sinα) ∗ (sinθcosα + cosθsinα)
⟹    y = y′cosθ + y′cosα ∗ (sinθ / sinα)
⟹    y = y′cosθ + x′sinα ∗ (sinθ / sinα) (from [2])
⟹    y = y′cosθ + x′sinθ
⟹    y = x′sinθ + y′cosθ

因此我们有，
      x = x′cosθ − y′sinθ
      y = x′sinθ + y′cosθ

上述方程式可以矩阵形式表示为：

我们的目标是找到(𝑥′,𝑦′). 因此，我们将矩阵移到另一侧，取其逆：

理解线性变换和基变换变化之间的区别非常重要。

接下来我们将看到这两种转换是如何关联的。

线性变换与基变换的关系

如果你观察，基矩阵的变化是线性变换矩阵的逆。这意味着，如果我们知道摄像机变换矩阵，即在世界上负责旋转和移动摄像机的矩阵，我们可以取其逆矩阵，这可以帮助我们找到摄像机上点的坐标。

平移

通过平移改变坐标

平移的想法很简单——只要有一点𝑃, 我们移动它一个偏移量来得到点𝑃′ 如下图所示：

在这里，我们移动点𝑃 坐标(𝑥, 𝑦) ，偏移量为(−𝑎, 𝑏) ，得到点𝑃′ (𝑥′, 𝑦′).。我们的目标是找到(𝑥′, 𝑦′).

从这个图来看，
      x′ = x - a
      y′ = y + b

我们不能将上述方程表示为矩阵乘法——至少不能用它们当前的表示形式。

诀窍是增加一个额外的维度，然后我们可以将平移表示为线性变换，如下所示：

用额外维度表示的坐标称为齐次坐标。

为了从齐次坐标中得到欧几里德坐标，我们只需除以最后一个元素，如下所示：

[x, y, 1] ≅ [x/1, y/1] = [x, y]

通常，我们在齐次空间中执行所有操作，因为这很容易处理，最后，当我们完成时，我们转换回欧几里德空间。稍后我们将看到更多细节。

通过平移改变基

就像我们之前看到的，在基变换的变化中，我们变换轴而不是点。在下面的示例中，我们移动𝑋𝑌 坐标轴偏移以获得𝑋′𝑌′. 我们的目标是找到𝑋′𝑌′下的点𝑃的坐标.

旧轴XY下P的坐标是(x,y)，新轴𝑋′𝑌′ 下P的坐标(𝑥′,𝑦′). 这里的偏移量是(𝑎, 𝑏).

从图上看,
      x′ = x - a
      y′ = y - b

同样，为了将上述方程表示为矩阵乘法，我们使用齐次坐标：

即使在平移中，线性变换和基变换的变化也是相反的。

摄像机外参矩阵

我们分别研究了旋转和平移；然而，我们可以使用如下所示的矩阵组合一次性执行这两个操作：

在这里𝑅 是旋转矩阵，形状是（3，3）和𝑂 是偏移量矩阵，形状是（3，1）。

通过求最终变换矩阵的逆，可以得到基矩阵的变化。我们称这个矩阵是摄像机外参矩阵E，形状是（4，4）

使用𝐸, 我们可以找到相机上任何一点的坐标。

自由度

相机外参矩阵自由度是6。三个旋转角度和沿X、Y、Z轴的三个偏移。

简化矩阵

我们可以看到相机外参矩阵的最后一行是0和1。它不会给转换增加任何价值，它的唯一目的是增加一个额外的维度——这意味着，正如我们将在下面的示例中看到的，我们可以删除最后一行。

实例

我们举一个实际操作的例子！

设置

包含所有代码的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绘图具有交互性。

实例直觉

在本例中，我们将首先创建旋转和平移的变换矩阵，将它们组合成一个矩阵，并使用它来变换相机。

然后，我们将通过对变换矩阵求逆来创建基矩阵的变化，并将其应用于一个点，并将其坐标从世界帧更改为相机帧。

下面是示例的笔记本，也可以在存储库中找到。

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

%matplotlib widget

import matplotlib.pyplot as plt
from utils import *

创建转换矩阵

# 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 parameters for the image plane
f = 2
img_size = (7, 7)

# 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)

R_, T_

(array([[ 0.70710678,  0.        , -0.70710678,  0.        ],
        [ 0.        ,  1.        ,  0.        ,  0.        ],
        [ 0.70710678,  0.        ,  0.70710678,  0.        ],
        [ 0.        ,  0.        ,  0.        ,  1.        ]]),
 array([[ 1.,  0.,  0.,  0.],
        [ 0.,  1.,  0., -8.],
        [ 0.,  0.,  1.,  0.],
        [ 0.,  0.,  0.,  1.]]))

转换并绘制

# 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=(6, 4))
ax = fig.add_subplot(111,projection='3d')

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

# plot the global basis and the transformed camera basis
ax = pr.plot_basis(ax)
ax = pr.plot_basis(ax, R, offset)

# plot the original and transformed image plane
ax.plot_surface(xx, yy, Z, alpha=0.75)
ax.plot_surface(xxt, yyt, Zt, alpha=0.75)

ax.set_title("camera transformation")
ax.set_xlabel("X-axis")
ax.set_ylabel("Y-axis")
ax.set_zlabel("Z-axis")

Text(0.5, 0, 'Z-axis')

创建基的变换矩阵

E = np.linalg.inv(R_ @ T_)

# remove last row of E
E = E[:-1, :]

进行坐标变换

cw = np.array([-1/np.sqrt(2), -8, 1/np.sqrt(2), 1]) # homogeneous coordinates of the point wrt the world
cc = E @ cw.reshape(4, 1) # coordinates of the point wrt the camera
cc = cc.flatten()

cc

array([0., 0., 1.])

让我们一步一步地分解：

• 首先，我们导入必要的库。utils.py文件包含所有必要的帮助函数。%matplotlib widget启用了ipympl后端，使我们能够使用绘图。

• 接下来，我们定义必要的参数，如角度、旋转顺序、平移偏移、焦距和图像平面的大小。焦距和图像平面仅用于演示目的，我们将在下一篇文章中详细讨论它们。

• 在这里，我们保持简单，关于标准Y轴旋转𝜋/4。然而，我们可以围绕任何轴进行任意数量的旋转。注意旋转的顺序。我们的平移偏移量是[0，-8，0]，沿Y轴8个单位。

• 使用这些参数，我们为旋转和平移创建变换矩阵。

• 接下来，我们使用变换矩阵转换最初位于原点的相机并绘制它。多亏了ipympl，图表是互动的。试着摆弄一下图表，用不同的视角来观看。

• 接下来，我们通过对变换矩阵求逆来创建基矩阵的变化，即相机外参矩阵。

• 最后，我们取一个世界坐标[-1/√2, -8, 1/√2]，然后应用基变换的变化，得到相机的坐标为[0, 0, 1]。这是有意义的，因为该点位于相机Z轴的正上方。

感谢阅读！

参考引用

计算机视觉导论——Udacity：https://classroom.udacity.com/courses/ud810

好消息！

小白学视觉知识星球

开始面向外开放啦👇👇👇

下载1：OpenCV-Contrib扩展模块中文版教程

在「小白学视觉」公众号后台回复：扩展模块中文教程，即可下载全网第一份OpenCV扩展模块教程中文版，涵盖扩展模块安装、SFM算法、立体视觉、目标跟踪、生物视觉、超分辨率处理等二十多章内容。


下载2：Python视觉实战项目52讲
在「小白学视觉」公众号后台回复：Python视觉实战项目，即可下载包括图像分割、口罩检测、车道线检测、车辆计数、添加眼线、车牌识别、字符识别、情绪检测、文本内容提取、面部识别等31个视觉实战项目，助力快速学校计算机视觉。


下载3：OpenCV实战项目20讲
在「小白学视觉」公众号后台回复：OpenCV实战项目20讲，即可下载含有20个基于OpenCV实现20个实战项目，实现OpenCV学习进阶。


交流群

欢迎加入公众号读者群一起和同行交流，目前有SLAM、三维视觉、传感器、自动驾驶、计算摄影、检测、分割、识别、医学影像、GAN、算法竞赛等微信群（以后会逐渐细分），请扫描下面微信号加群，备注：”昵称+学校/公司+研究方向“，例如：”张三 + 上海交大 + 视觉SLAM“。请按照格式备注，否则不予通过。添加成功后会根据研究方向邀请进入相关微信群。请勿在群内发送广告，否则会请出群，谢谢理解~