Python计算机视觉编程——第四章 照相机模型与增强现实

1 针孔照相机模型

针孔照相机模型是计算机视觉中广泛使用的照相机模型，该模型简单且具有足够的精确度。该照相机从一个小孔采集射到暗箱内部的光线。在针孔照相机模型中，在光线投影到图像平面之前，从唯一一个点经过，也就是照相机中心C。

在针孔照相机中，三维点X投影为图像点x，如下所示：
$$\lambda \mathbf{x}=P\mathbf{X}$$

这里，3×4 的矩阵$P$为照相机矩阵(或投影矩阵)。注意，在齐次坐标系中，三维点 X 的坐标由 4 个元素组成，$\mathbf{X}=[X,Y,Z,W]$。

1.1 照相机矩阵

照相机矩阵可以分解为：$$P=K[R|t]$$

其中，$R$是描述照相机方向的旋转矩阵，$t$是描述照相机中心位置的三维平移向量，内标定矩阵$K$描述照相机的投影性质。标定矩阵通常情况下可以写为：$\mathbf{K}=\left[\begin{array}{ccc}\alpha f& s& c_x\\ 0& f& c_y\\ 0& 0& 1\end{array}\right]$

图像平面和照相机中心间的距离为焦距$f$。在大多数情况下，$s$可以设置成 0。也就是说：
$$K=\left[\begin{array}{ccc}{f}_x& 0& c_x\\ 0& f_y& c_y\\ 0& 0& 1\end{array}\right]$$

使用了另外的记号$f_x$ 和$f_y$, 两者关系为$f_x=\alpha f_y$。

通常情况下，默认设置 $\alpha =1$。经过这些假设，标定矩阵变为.
$$\mathbf{K}=\left[\begin{array}{ccc}f& 0& c_x\\ 0& f& c_y\\ 0& 0& 1\end{array}\right]$$

1.2 三维点的投影

• 创建照相机类（处理对照相机和投影建模所需操作）
• 使用Model Housing数据集将这些三维点投影到图像平面执行绘制操作，并且进行增量旋转。

from scipy import linalg
import  numpy as np
import camera
from matplotlib import pyplot as plt
class Camera(object):
  
  def __init__(self,P):
    self.P = P
    self.K = None # 标定矩阵
    self.R = None # 旋转
    self.t = None # 平移
    self.c = None # 照相机中心
  def project(self,X):
    x = np.dot(self.P,X)
    for i in range(3):
     x[i] /= x[2]
    return x

def rotation_matrix(a):
    R = np.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 = np.loadtxt('house.p3d').T
points = np.vstack((points, np.ones(points.shape[1])))
P = np.hstack((np.eye(3), np.array([[0], [0], [-10]])))
cam = Camera(P)
x = cam.project(points)

r = 0.05 * np.random.rand(3)
rot = rotation_matrix(r)
plt.figure()
plt.subplot(1, 2, 1)
plt.plot(x[0], x[1], 'k.')

plt.subplot(1, 2, 2)
for t in range(20):
  cam.P = np.dot(cam.P, rot)
  x = cam.project(points)
  plt.plot(x[0], x[1], 'k.')

plt.show()

结果如下图

1.3 照相机矩阵的分解

矩阵分块操作称为因子分解，这里介绍一种矩阵因子分解的方法——RQ因子分解。
将下面代码添加到Camera类中：

def factor(self):

  K,R = linalg.rq(self.P[:,:3])
  T = diag(sign(diag(K)))
  if linalg.det(T) < 0:
    T[1,1] *= -1
  self.K = dot(K, T)
  self.R = dot(T, R)  
  self.t = dot(linalg.inv(self.K), self.P[:, 3])
  return self.K, self.R, self.t

使用如下代码

K = array([[1000,0,500],[0,1000,300],[0,0,1]])
tmp = rotation_matrix([0,0,1])[:3,:3]
Rt = hstack((tmp,array([[50],[40],[30]])))
cam = Camera(dot(K,Rt))
print (K,Rt)
print (cam.factor)

得到结果：

RQ因子分解的结果并不是唯一的。在该因子分解中，分解的结果存在符号二义性。如果需要，可以在求解到的结果中加入变换T来改变符号。

1.4 计算照相机中心

给定照相机投影矩阵$P$,我们可以计算出空间上照相机的所在位置。照相机的中心
$\mathbb{C}$,是一个三维点，满足约束$P\mathbb{C}=0$。对于投影矩阵为$P=K[R|t]$的照相机，有：
$$K[R|t]\mathbf{C}=KR\mathbf{C}+Kt=0$$
照相机的中心可以由下述式子来计算：
$$\mathbf{C}=-R^{T}t$$
可使用下面代码来计算照相机中心：

    def center(self):
        if self.c is not None:
            return self.c
        else:
            self.factor()
            self.c = -dot(self.R.T,self.t)
            return self.c

2 照相机标定

标定照相机是指计算出该照相机的内参数，指计算矩阵K。标定照相机的标准方法为，拍摄多幅平面棋盘模式的图像，然后进行处理计算。

2.1 一个简单的标定方法

大多数参数可以使用基本的假设来设定，比较难处理的是获得正确的焦距。需要准备一个平面矩形的标定物体(一个书本即可)、用于测量的卷尺和直尺，以及一个平面。
具体操作步骤：

• 测量你选定矩形标定物体的边长 d$X$和 d$Y;$
• 将照相机和标定物体放置在平面上，使得照相机的背面和标定物体平行，同时物体位于照相机图像视图的中心，你可能需要调整照相机或者物体来获得良好的对齐效果，
• 测量标定物体到照相机的距离 dZ,
• 拍摄一副图像来检验该设置是否正确，即标定物体的边要和图像的行和列对齐；
• 使用像素数来测量标定物体图像的宽度和高度 dx 和 dy。

使用下面的相似三角形关系得到焦距：
$f_x=\frac{\mathrm{d}x}{\mathrm{d}X}\mathrm{d}Z,f_y=\frac{\mathrm{d}y}{\mathrm{d}Y}\mathrm{d}Z$

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

如果图像中包含平面状的标记物体，并且已经对照相机进行了标定，可以计算出照相机的姿态。

先提取两幅图像的SIFT特征，然后使用RANSAC算法估计单应性矩阵：
通过下面代码获得单应性矩阵，该矩阵将一幅图像中标记物上的点映射到另一幅图像中的对应点。

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):
  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))
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]))
box_trans = homography.normalize(dot(H,box_cam1))
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))
point = array([1,1,0,1]).T
print (homography.normalize(dot(dot(H,cam1.P),point)))
print (cam2.project(point))

4 增强现实

增强现实是将物体和相应信息放置在图像数据上的一系列操作的总称。将会使用两个工具包：PyGame和PyOpenGL。

4.1 PyGame和PyOpenGL

PyGame是流行的游戏开发工具包，用于简单地处理显示窗口，输入设备，事件以及其他内容。
PyOpenGL 是 OpenGL 图形编程的 Python 绑定接口。OpenGL 可以安装在几乎所有的系统上，并且具有很好的图形性能。OpenGL具有跨平台性，能够在不同的操作系统之间工作。
为了使用PyGame和PyOpenGL，需要载入下面命令：

from OpenGL.GL import *
from OpenGL.GLU import *
import pygame, pygame.image
from pygame.locals import *

4.2 从照相机矩阵到OpenGL格式

OpenGL使用4*4矩阵来表示变换。照相机与场景的变换分成了两个矩阵， GL_PROJECTION 矩阵和 GL_MODELVIEW 矩阵。GL_PROJECTION 矩阵处理图像成像的性质，等价于我们的内标定矩阵$K$。GL_MODELVIEW 矩阵处理物体和照相机之间的三维变换关系，对应于我们照相机矩阵中的$R$和$t$部分。
下面的函数将照相机参数转换为OpenGL中的投影矩阵：

def set_projection_from_camera(K): 
	glMatrixMode(GL_PROJECTION) 
	glLoadIdentity()
	fx = K[0,0] 
	fy = K[1,1] 
	fovy = 2*math.atan(0.5*height/fy)*180/math.pi 
	aspect = (width*fy)/(height*fx)
 
	near = 0.1 
	far = 100.0
 
	gluPerspective(fovy,aspect,near,far) 
	glViewport(0,0,width,height)

下面函数实现如何获得移除标定矩阵后的3*4针孔照相机矩阵，并创建一个模拟试图：

def set_modelview_from_camera(Rt):
    glMatrixMode(GL_MODELVIEW)
    glLoadIdentity()
    Rx = np.array([[1,0,0],[0,0,-1],[0,1,0]])
    R = Rt[:,:3]
    U,S,V = linalg.svd(R)
    R = np.dot(U,V)
    R[0,:] = -R[0,:]
    t = Rt[:,3]
    M = np.eye(4)
    M[:3,:3] = np.dot(R,Rx)
    M[:3,3] = t
    M = M.T
    m = M.flatten()
    glLoadMatrixf(m)

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

4.3 在图像中放置虚拟物体

第一件事是将图像作为背景添加进来。在OpenGL中，通过创建一个四边形的方式来完成，该四边形为整个视觉。
下面函数载入一幅图像，将其转换成一个OpenGL纹理，将该纹理放置在四边形上：

def draw_background(imname):
  
  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)
  glBindTexture(GL_TEXTURE_2D,glGenTextures(1))
  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(1)