Python计算机视觉第4章-照相机模型与增强现实

4.1　针孔照相机模型

针孔照相机模型（有时称为射影照相机模型）是计算机视觉中广泛使用的照相机模型。对于大多数应用来说，针孔照相机模型简单，并且具有足够的精确度。这个名字来源于一种类似暗箱机的照相机。该照相机从一个小孔采集射到暗箱内部的光线。在针孔照相机模型中，在光线投影到图像平面之前，从唯一一个点经过，也就是照相机中心C。图4-1为从照相机中心前画出图像平面的图解。事实上，在真实的照相机中，图像平面位于照相机中心之后，但是照相机的模型和图4-1的模型是一样的。
在这里插入图片描述

图4-1 针孔照相机模型。图像点x是由图像平面与连接三维点X和照相机中心C的直线相交而成的。虚线表示该照相机的光学坐标轴

由图像坐标轴和三维坐标系中的x轴和y轴对齐平行的假设，可以得出针孔照相机的投影性质。照相机的光学坐标轴和 $z$ 轴一致，该投影几何可以简化成相似三角形。在投影之前通过旋转和平移变换，对该坐标系加入三维点，会出现完整的投影变换。

4.1.1　照相机矩阵

照相机矩阵可以分解为：
$P = K [R ∣ t]$
其中， $R$ 是描述照相机方向的旋转矩阵， $t$ 是描述照相机中心位置的三维平移向量，内标定矩阵 $K$ 描述照相机的投影性质。

在计算机视觉和摄影测量中，照相机矩阵（Camera Matrix）也被称为内参矩阵（Intrinsic Matrix），是描述照相机内在参数的一个3x3矩阵。它用于将三维空间中的点映射到相机的二维图像平面上。

照相机矩阵的形式如下：

$\begin{bmatrix} f_x & 0 & c_x \\ 0 & f_y & c_y \\ 0 & 0 & 1 \end{bmatrix}$

其中：

$f_x$ 和 $f_y$ 是图像在 x 和 y 方向上的焦距，通常以像素为单位。
$c_x$ 和 $c_y$ 是光轴在图像平面上的主点（Principal Point）坐标，即图像的中心位置。
1 是一个尺度因子。

这个矩阵 $K$ 用于描述相机的内在几何和光学特性，如焦距、光轴位置等。它是相机标定过程的核心部分，通过标定可以求得这些参数，从而准确地进行三维重建、图像矫正、测量等任务。

4.1.2　三维点的投影

下面来创建照相机类，用来处理对照相机和投影建模所需要的全部操作：

from scipy import linalg

class Camera(object):
    """表示针孔照相机的类"""

    def __init__(self, P):
        """初始化 P = K[R|t] 照相机模型"""
        self.P = P
        self.K = None  # 标定矩阵
        self.R = None  # 旋转
        self.t = None  # 平移
        self.c = None  # 照相机中心

    def project(self, X):
        """X（4×n的数组）的投影点，并且进行坐标归一化"""
        x = dot(self.P, X)
        for i in range(3):
            x[i] /= x[2]
        return x

下面的例子展示如何将三维中的点投影到图像视图中。在这个例子中，将使用牛津多视图数据集中的“Model Housing”数据集。可以从http://www.robots.ox.ac.uk/~vgg/data/data-mview.html 下载。

import camera

# 载入点
points = loadtxt('house.p3d').T
points = vstack((points, ones(points.shape[1])))

# 设置照相机参数
P = hstack((eye(3), array([[0], [0], [-10]])))
cam = camera.Camera(P)
x = cam.project(points)

# 绘制投影
figure()
plot(x[0], x[1], 'k.')
show()

首先，使用齐次坐标来表示这些点。然后使用一个投影矩阵来创建Camera对象将这些三维点投影到图像平面并执行绘制操作，输出结果如图4-2中间图像所示。

为了研究照相机的移动会如何改变投影的效果，可以使用下面的代码。该代码围绕一个随机的三维向量，进行增量旋转的投影。

# 创建变换
r = 0.05 * random.rand(3)
rot = camera.rotation_matrix(r)

# 旋转矩阵和投影
figure()
for t in range(20):
    cam.P = dot(cam.P, rot)
    x = cam.project(points)
    plot(x[0], x[1], 'k.')

show()

在上面的代码中，使用了rotation_matrix()函数，该函数能够创建围绕一个向量进行三维旋转的旋转矩阵

def rotation_matrix(a):
    """创建一个用于围绕向量 a 轴旋转的三维旋转矩阵"""
    R = eye(4)
    R[:3, :3] = linalg.expm([[0, -a[2], a[1]], 
                             [a[2], 0, -a[0]], 
                             [-a[1], a[0], 0]])
    return R

结合以上代码段，完整代码如下：

from scipy import linalg
import random
from numpy import dot, eye, array, loadtxt, vstack, ones
import matplotlib.pyplot as plt

class Camera(object):
    """表示针孔照相机的类"""

    def __init__(self, P):
        """初始化 P = K[R|t] 照相机模型"""
        self.P = P
        self.K = None  # 标定矩阵
        self.R = None  # 旋转
        self.t = None  # 平移
        self.c = None  # 照相机中心

    def project(self, X):
        """X（4×n的数组）的投影点，并且进行坐标归一化"""
        x = dot(self.P, X)
        for i in range(3):
            x[i] /= x[2]
        return x

def rotation_matrix(a):
    """创建一个用于围绕向量 a 轴旋转的三维旋转矩阵"""
    R = eye(4)
    R[:3, :3] = linalg.expm([[0, -a[2], a[1]], 
                             [a[2], 0, -a[0]], 
                             [-a[1], a[0], 0]])
    return R

# 载入点
points = loadtxt('house.p3d').T
points = vstack((points, ones(points.shape[1])))

# 设置照相机参数
P = hstack((eye(3), array([[0], [0], [-10]])))
cam = Camera(P)
x = cam.project(points)

# 绘制投影
plt.figure()
plt.plot(x[0], x[1], 'k.')
plt.show()

# 创建变换
r = 0.05 * random.rand(3)
rot = rotation_matrix(r)

# 旋转矩阵和投影
plt.figure()
for t in range(20):
    cam.P = dot(cam.P, rot)
    x = cam.project(points)
    plt.plot(x[0], x[1], 'k.')

plt.show()

运行结果如下：

在这里插入图片描述

实验图1 实验结果

4.1.3　照相机矩阵的分解

将使用一种矩阵因子分解的方法，称为RQ因子分解。

将下面的方法添加到Camera类中：

def factor(self):
    """将照相机矩阵分解为 K、R、t，其中，P = K[R|t]"""
    # 分解前3×3的部分
    K, R = linalg.rq(self.P[:, :3])
    
    # 将 K 的对角线元素设为正值
    T = diag(sign(diag(K)))
    if linalg.det(T) < 0:
        T[1, 1] *= -1
    
    self.K = dot(K, T)
    self.R = dot(T, R)  # T 的逆矩阵为其自身
    self.t = dot(linalg.inv(self.K), self.P[:, 3])
    
    return self.K, self.R, self.t

RQ因子分解的结果并不是唯一的。在该因子分解中，分解的结果存在符号二义性。由于需要限制旋转矩阵R为正定的（否则，旋转坐标轴即可），所以如果需要，可以在求解到的结果中加入变换T来改变符号。

以下为照相机矩阵分解的代码：

import camera
from numpy import array, hstack, dot

K = array([[1000, 0, 500],
           [0, 1000, 300],
           [0, 0, 1]])

tmp = camera.rotation_matrix([0, 0, 1])[:3, :3]
Rt = hstack((tmp, array([[50], [40], [30]])))
cam = camera.Camera(dot(K, Rt))

print(K, Rt)
print(cam.factor())

4.2　照相机标定

标定照相机是指计算出该照相机的内参数。例子中，是指计算矩阵K。如果你的应用要求高精度，那么可以扩展该照相机模型,使其包含径向畸变和其他条件。

标定照相机的标准方法是，拍摄多幅平面棋盘模式的图像，然后进行处理计算。

4.2.1　一个简单的标定方法

这里将要介绍一个简单的照相机标定方法。大多数参数可以使用基本的假设来设定（正方形垂直的像素，光心位于图像中心），比较难处理的是获得正确的焦距。对于这种标定方法，你需要准备一个平面矩形的标定物体（一个书本即可）、用于测量的卷尺和直尺，以及一个平面。下面是具体操作步骤：
1、测量你选定矩形标定物体的边长 $d X$ 和 $d Y$ ；
2、将照相机和标定物体放置在平面上，使得照相机的背面和标定物体平行，同时物体位于照相机图像视图的中心，你可能需要调整照相机或者物体来获得良好的对齐效果；
3、测量标定物体到照相机的距离 $d Z$ ；
4、拍摄一副图像来检验该设置是否正确，即标定物体的边要和图像的行和列对齐；
5、使用像素数来测量标定物体图像的宽度和高度 $d x$ 和 $d y$

4.3　以平面和标记物进行姿态估计

如果图像中包含平面状的标记物体，并且已经对照相机进行了标定，那么可以计算出照相机的姿态（旋转和平移）。这里的标记物体可以为对任何平坦的物体。

使用下面的代码来提取两幅图像的SIFT特征，然后使用RANSAC算法稳健地估计单应性矩阵：

import homography
import camera
import sift

# 计算特征
sift.process_image('book_frontal.JPG', 'im0.sift')
l0, d0 = sift.read_features_from_file('im0.sift')

sift.process_image('book_perspective.JPG', 'im1.sift')
l1, d1 = sift.read_features_from_file('im1.sift')

# 匹配特征，并计算单应性矩阵
matches = sift.match_twosided(d0, d1)
ndx = matches.nonzero()[0]
fp = homography.make_homog(l0[ndx, :2].T)

ndx2 = [int(matches[i]) for i in ndx]
tp = homography.make_homog(l1[ndx2, :2].T)

model = homography.RansacModel()
H = homography.H_from_ransac(fp, tp, model)

这样就得到了单应性矩阵，为了检验单应性矩阵结果的正确性，需要将一些简单的三维物体放置在标记物上，这里使用一个立方体。可以使用下面的函数来产生立方体上的点：

def cube_points(c, wid):
    """创建用于绘制立方体的一个点列表（前5个点是底部的正方形，一些边重合了）"""
    p = []
    
    # 底部
    p.append([c[0] - wid, c[1] - wid, c[2] - wid])
    p.append([c[0] - wid, c[1] + wid, c[2] - wid])
    p.append([c[0] + wid, c[1] + wid, c[2] - wid])
    p.append([c[0] + wid, c[1] - wid, c[2] - wid])
    p.append([c[0] - wid, c[1] - wid, c[2] - wid])  # 为了绘制闭合图像，和第一个相同

    # 顶部
    p.append([c[0] - wid, c[1] - wid, c[2] + wid])
    p.append([c[0] - wid, c[1] + wid, c[2] + wid])
    p.append([c[0] + wid, c[1] + wid, c[2] + wid])
    p.append([c[0] + wid, c[1] - wid, c[2] + wid])
    p.append([c[0] - wid, c[1] - wid, c[2] + wid])  # 为了绘制闭合图像，和第一个相同

    # 竖直边
    p.append([c[0] - wid, c[1] - wid, c[2] + wid])
    p.append([c[0] - wid, c[1] + wid, c[2] + wid])
    p.append([c[0] - wid, c[1] + wid, c[2] - wid])
    p.append([c[0] + wid, c[1] + wid, c[2] - wid])
    p.append([c[0] + wid, c[1] + wid, c[2] + wid])
    p.append([c[0] + wid, c[1] - wid, c[2] + wid])
    p.append([c[0] + wid, c[1] - wid, c[2] - wid])
    
    return array(p).T

有了单应性矩阵和照相机的标定矩阵，可以使用以下代码得出两个视图间的相对变换：

# 计算照相机标定矩阵
K = my_calibration((747, 1000))

# 位于边长为0.2，z=0平面上的三维点
box = cube_points([0, 0, 0.1], 0.1)

# 投影第一幅图像上底部的正方形
cam1 = camera.Camera(hstack((K, dot(K, array([[0], [0], [-1]])))))

# 底部正方形上的点
box_cam1 = cam1.project(homography.make_homog(box[:, :5]))

# 使用 H 将点变换到第二幅图像中
box_trans = homography.normalize(dot(H, box_cam1))

# 从 cam1 和 H 中计算第二个照相机矩阵
cam2 = camera.Camera(dot(H, cam1.P))
A = dot(linalg.inv(K), cam2.P[:, :3])
A = array([A[:, 0], A[:, 1], cross(A[:, 0], A[:, 1])]).T
cam2.P[:, :3] = dot(K, A)

# 使用第二个照相机矩阵投影
box_cam2 = cam2.project(homography.make_homog(box))

# 测试：将点投影在 z=0 上，应该能够得到相同的点
point = array([1, 1, 0, 1]).T
print(homography.normalize(dot(dot(H, cam1.P), point)))
print(cam2.project(point))

4.4　增强现实

增强现实（Augmented Reality，AR）是将物体和相应信息放置在图像数据上的一系列操作的总称。最经典的例子是放置一个三维计算机图形学模型，使其看起来属于该场景；如果在视频中，该模型会随着照相机的运动很自然地移动。

4.4.1　PyGame和PyOpenGL

PyGame 是非常流行的游戏开发工具包，它可以非常简单地处理显示窗口、输入设备、事件，以及其他内容。PyGame是开源的，可以从http://www.pygame.org/下载。事实上，它是一个Python绑定的SDL游戏引擎。

PyOpenGL 是 OpenGL 图形编程的Python绑定接口。OpenGL可以安装在几乎所有的系统上，并且具有很好的图形性能。OpenGL具有跨平台性，能够在不同的操作系统之间工作。关于OpenGL的更多信息，参见http://www.opengl.org/。

4.4.2　从照相机矩阵到OpenGL格式

OpenGL 使用4×4的矩阵来表示变换（包括三维变换和投影）。这和使用的3×4照相机矩阵略有差别。但是，照相机与场景的变换分成了两个矩阵， GL_PROJECTION 矩阵和GL_MODELVIEW矩阵。GL_PROJECTION矩阵处理图像成像的性质，等价于的内标定矩阵K。GL_MODELVIEW矩阵处理物体和照相机之间的三维变换关系，对应于照相机矩阵中的R和t部分。一个不同之处是，假设照相机为坐标系的中心，GL_MODELVIEW矩阵实际上包含了将物体放置在照相机前面的变换。

假设已经获得了标定好的照相机，即已知标定矩阵K，下面的函数可以将照相机参数转换为OpenGL中的投影矩阵：

def set_projection_from_camera(K):
    """从照相机标定矩阵中获得视图"""
    glMatrixMode(GL_PROJECTION)
    glLoadIdentity()
    
    fx = K[0, 0]
    fy = K[1, 1]
    
    fovy = 2 * arctan(0.5 * height / fy) * 180 / pi
    aspect = (width * fy) / (height * fx)
    
    # 定义近的和远的剪裁平面
    near = 0.1
    far = 100.0
    
    # 设定透视
    gluPerspective(fovy, aspect, near, far)
    glViewport(0, 0, width, height)

下面的函数实现如何获得移除标定矩阵后的3×4针孔照相机矩阵（将P和K-1相乘），并创建一个模拟视图：

def set_modelview_from_camera(Rt):
    """从照相机姿态中获得模拟视图矩阵"""
    glMatrixMode(GL_MODELVIEW)
    glLoadIdentity()
    
    # 围绕 x 轴将茶壶旋转 90 度，使 z 轴向上
    Rx = array([[1, 0, 0],
                [0, 0, -1],
                [0, 1, 0]])
    
    # 获得旋转的最佳逼近
    R = Rt[:, :3]
    U, S, V = linalg.svd(R)
    R = dot(U, V)
    R[0, :] = -R[0, :]  # 改变 x 轴的符号
    
    # 获得平移量
    t = Rt[:, 3]
    
    # 获得 4×4 的模拟视图矩阵
    M = eye(4)
    M[:3, :3] = dot(R, Rx)
    M[:3, 3] = t
    
    # 转置并压平以获取列序数值
    M = M.T
    m = M.flatten()
    
    # 将模拟视图矩阵替换为新的矩阵
    glLoadMatrixf(m)

该操作使用SVD分解方法，旋转矩阵的最佳逼近可以通过 $R=UV^T$ 来获得

4.4.3　在图像中放置虚拟物体

我们需要做的第一件事是将图像（打算放置虚拟物体的图像）作为背景添加进来。在OpenGL中，该操作可以通过创建一个四边形的方式来完成，该四边形为整个视图。完成该操作最简单的方式是绘制出四边形，同时将投影和模拟试图矩阵重置，使得每一维的坐标范围在-1到1之间。

下面的函数可以载入一幅图像，然后将其转换成一个OpenGL纹理，并将该纹理放置在四边形上：

def draw_background(imname):
    """使用四边形绘制背景图像"""
    # 载入背景图像（应该是 .bmp 格式），转换为 OpenGL 纹理
    bg_image = pygame.image.load(imname).convert()
    bg_data = pygame.image.tostring(bg_image, "RGBX", 1)
    
    glMatrixMode(GL_MODELVIEW)
    glLoadIdentity()
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    
    # 绑定纹理
    glEnable(GL_TEXTURE_2D)
    texture_id = glGenTextures(1)
    glBindTexture(GL_TEXTURE_2D, texture_id)
    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, width, height, 0, GL_RGBA, GL_UNSIGNED_BYTE, bg_data)
    glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST)
    glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST)
    
    # 创建四方形填充整个窗口
    glBegin(GL_QUADS)
    glTexCoord2f(0.0, 0.0); glVertex3f(-1.0, -1.0, -1.0)
    glTexCoord2f(1.0, 0.0); glVertex3f( 1.0, -1.0, -1.0)
    glTexCoord2f(1.0, 1.0); glVertex3f( 1.0,  1.0, -1.0)
    glTexCoord2f(0.0, 1.0); glVertex3f(-1.0,  1.0, -1.0)
    glEnd()
    
    # 清除纹理
    glDeleteTextures(texture_id)

该函数首先使用PyGame中的一些函数来载入一幅图像，将其序列化为能够在PyOpenGL 中使用的原始字符串表示。然后，重置模拟视图，清除颜色和深度缓存。接下来，绑定这个纹理，使其能够在四边形和指定插值中使用它。四边形是在每一维分别为-1和1的点上定义的。注意，纹理图像的坐标是从0到1。最后，清除该纹理，避免其干扰之后准备绘制的图像。

4.4.4　综合集成

下面的完整脚本可以生成如图4-5所示的图像：

from OpenGL.GL import *
from OpenGL.GLU import *
from OpenGL.GLUT import *
import pygame
import pygame.image
from pygame.locals import *
import pickle

width, height = 1000, 747

def setup():
    """设置窗口和 pygame 环境"""
    pygame.init()
    pygame.display.set_mode((width, height), OPENGL | DOUBLEBUF)
    pygame.display.set_caption('OpenGL AR demo')

# 载入照相机数据
with open('ar_camera.pkl', 'rb') as f:
    K = pickle.load(f)
    Rt = pickle.load(f)

setup()
draw_background('book_perspective.bmp')
set_projection_from_camera(K)
set_modelview_from_camera(Rt)
draw_teapot(0.02)

while True:
    event = pygame.event.poll()
    if event.type in (QUIT, KEYDOWN):
        break

    pygame.display.flip()