计算机视觉编程 第四章

前言

本专栏按《python计算机视觉编程 ——Jan Erik Solem》一书为参考，第四章介绍照相机的相关内容，包含照相机建模和投影建模、计算照相机参数、增强现实等技术

4.1 针孔照相机模型

针孔照相机模型是一种用于描述相机成像原理的简化模型，用它来解释相机成像最合适不过，它基于光线传播的基本原理，主要思想是，相机的镜头被简化为一个微小的光圈，光线从被摄物体通过一个微小的孔（针孔），也就是照相机中心C进入相机内部。这些光线穿过针孔后，在相机内部的感光材料（胶片或图像传感器）上投影成图像

针孔照相机中，三维点$\mathbf{X}$投影为图像点$\mathbf{x}$，如下所示$$\lambda \mathbf{x}=\mathbf{PX}$$其中$\mathbf{P}$为照相机矩阵，三维点坐标$\mathbf{X}=[X,Y,Z,W]$，标量$\lambda$是三维点的逆深度

4.1.1 照相机矩阵

照相机矩阵$\mathbf{P}$可分解为$$\mathbf{P}=K[\mathbf{R}|t]$$其中$\mathbf{R}$描述照相机方向的旋转矩阵，$t$描述照相机中心位置的三维平移向量，内标定矩阵$\mathbf{K}$描述照相机的投影性质，包含了相机的内在参数，主要包括焦距、主点（图像中心点）和畸变系数等。它可以表示为$$K=\left[\begin{array}{ccc}\alpha f& s& c_x\\ 0& f& c_y\\ 0& 0& 1\end{array}\right]$$$f$表示图像与照相机中心的焦距，$s$为倾斜参数，通常可以设置为0，当像素为非正方形时，使用纵横比例参数$\alpha$，一般情况为1，光心坐标$\mathbf{c}=[c_x,c_y]$表示光线坐标轴与图像的交点，通常为高和宽的一般。因此内标定矩阵一般情况下为$$K=\left[\begin{array}{ccc}f& 0& c_x\\ 0& f& c_y\\ 0& 0& 1\end{array}\right]$$以上参数中唯一未知的是焦距$f$

4.1.2 三维点的投影

创建照相机类的代码如下

from scipy import linalg
from pylab import  *
 
class Camera(object):
    """ Class for representing pin-hole cameras. """
    
    def __init__(self,P):
        """ Initialize P = K[R|t] camera model. """
        self.P = P
        self.K = None # calibration matrix
        self.R = None # rotation
        self.t = None # translation
        self.c = None # camera center
        
    
    def project(self,X):
        """    Project points in X (4*n array) and normalize coordinates. """
        
        x = dot(self.P,X)
        for i in range(3):
            x[i] /= x[2]    
        return x

由于书上数据集的来源网站无法获取，仅给出代码示例。主要步骤是先载入数据点，再使用一个投影矩阵创建对象，并通过以上的Camera函数进行照相机参数的设置，最后绘制投影。同时为研究照相机移动对投影效果的影响，可以使用下列代码

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

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

4.1.3 照相机矩阵的分解

照相机矩阵的分解允许我们从照相机矩阵$\mathbf{P}$中恢复出有用的信息，如相机的位置$t$、姿态$\mathbf{R}$以及内在参数。通常，照相机矩阵可以分解为视图矩阵和投影矩阵，同时也可以从投影矩阵中分解出内标定矩阵。矩阵分块操作称为因子分解，下面使用RQ因子分解

    def factor(self):
        """    Factorize the camera matrix into K,R,t as P = K[R|t]. """
        
        # factor first 3*3 part
        K,R = linalg.rq(self.P[:,:3])
        
        # make diagonal of K positive
        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 is its own inverse
        self.t = dot(linalg.inv(self.K),self.P[:,3])
        
        return self.K, self.R, self.t

RQ因子分解结果并不唯一，结果存在二义性，可以通过在得到的结果中加入变换$\mathbf{T}$改变符号

from numpy import *
from pylab import *
from PIL import Image

from 视觉编程 import camera

import numpy as np
from scipy.spatial.transform import Rotation as R

# 定义内标定矩阵 K
K = np.array([[1000, 0, 500],
              [0, 1000, 300],
              [0, 0, 1]])

# 生成一个绕 [0, 0, 1] 轴的旋转矩阵
rotation_matrix = R.from_euler('xyz', [0, 0, 1], degrees=True)
tmp = rotation_matrix.as_matrix()[:3, :3]

# 构建外参矩阵 Rt，其中包括旋转矩阵和平移向量
Rt = hstack((tmp, np.array([[50], [40], [30]])))

# 计算相机矩阵
camera_matrix = np.dot(K, Rt)

# 输出内标定矩阵 K 和外参矩阵 Rt
print(K)
print(Rt)

# 使用相机矩阵构建一个相机对象
cam = camera.Camera(np.dot(K, Rt))

# 输出相机内参和外参的分解结果
print(cam.factor())

第一个是内标定矩阵K直接介绍过了，第二个矩阵Rt是相机的外参矩阵，它包含了相机的旋转矩阵和平移向量：前三列 [0.54030231 -0.84147098 0.] 表示旋转矩阵，它描述了相机的姿态，如旋转角度。第四列 [50. 40. 30.] 是平移向量，表示相机在世界坐标系中的位置。最后第一个矩阵是内标定矩阵 K 的分解结果，包含了焦距和主点的坐标。第二个矩阵是外参矩阵 Rt 的分解结果，包含了旋转矩阵（姿态）和一个单位平移向量（由于最后一列是 [0, 0, 0, 1]）。第三个矩阵是平移向量，表示相机在世界坐标系中的位置。注意到两个t矩阵的第二个元素符号不相同

4.1.4 计算照相机中心

通过照相机投影矩阵$\mathbf{P}$能计算空间上照相机的位置，中心$\mathbf{C}$满足$\mathbf{PC}=0$，对于$\mathbf{P}=K[\mathbf{R}|t]$的照相机有$$\mathbf{K}[\mathbf{R}|\mathbf{t}]\mathbf{C}=\mathbf{K}\mathbf{RC}+\mathbf{Kt}=0$$$$\mathbf{C}=-{\mathbf{R}}^T\mathbf{t}$$照相机中心与内标定矩阵无关，下面是定义中心的函数代码

    def center(self):
        """    Compute and return the camera center. """
    
        if self.c is not None:
            return self.c
        else:
            # compute c by factoring
            self.factor()
            self.c = -dot(self.R.T,self.t)
            return self.c

4.2 照相机标定

照相机标定的目标是确定相机的内参数和外参数，以及可能的畸变参数，从而实现准确的图像到世界坐标的映射，或者反之。标准方法是用多幅平面棋盘模式的图像再进行处理

可以使用一种基于已知物体尺寸的简单标定方法。这种方法被称为尺寸标定，它适用于一些特定的情况，例如你拥有一个已知物体的尺寸，并且可以测量它在图像中的像素尺寸

1. 选择已知物体：选择一个在现实世界中有已知尺寸的物体，例如标定板、纸片、图钉等。确保你知道这个物体的真实尺寸（例如，物体的宽度、高度或直径）
2. 拍摄标定图像：使用相机拍摄包含已知物体的图像，确保图像清晰且物体完整可见
3. 测量像素尺寸：在标定图像中，使用图像处理技术测量已知物体在图像中的像素尺寸。例如，你可以测量物体的宽度和高度，或者物体中心到图像边界的距离
4. 计算焦距：使用物体的已知尺寸和测得的像素尺寸，可以通过简单的比例计算来估计相机的焦距。假设已知物体的真实尺寸为 S，在图像中的像素尺寸为 P，焦距为 f，则可以使用以下关系估计焦距：$f=(P\ast D)/S$ 其中，D 是相机到物体的距离。注意，这里假设焦点位于物体上方

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

下面使用一个例子来演示姿态估计。先使用下列代码提取两图的SIFT特征，然后使用之前的RANSAC算法估计单应性矩阵

sift.process_image('./filelist/cam1.jpeg', 'im0.sift')
l0, d0 = sift.read_features_from_file('im0.sift')

sift.process_image('./filelist/cam2.jpeg', 'im1.sift')
l1, d1 = sift.read_features_from_file('im1.sift')

# match features and estimate homography
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, inliers = homography.H_from_ransac(fp, tp, model)

下面是构建一个立方体函数的代码，之后使用它来检验单应性矩阵的正确性

def cube_points(c, wid):
    """ Creates a list of points for plotting
        a cube with plot. (the first 5 points are
        the bottom square, some sides repeated). """
    p = []
    # bottom
    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])  # same as first to close plot

    # top
    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])  # same as first to close plot

    # vertical sides
    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

但是在尝试了n多次的实验之后，换了多张照片后，不是得到了如下的报错就是产生了奇奇怪怪的结果，令我匪夷所思


但是，经过查找，发现使用 OpenCV 库计算图像中的特征点，并进行特征匹配得到的效果更好，其中主要修改的代码如下，代码参考’别来这个网址’的博客: Python机器学习实战第四章 照相机模型与增强现实

# 计算特征点
sift = cv2.SIFT_create()
l0, d0 = sift.detectAndCompute(im0, None)
l1, d1 = sift.detectAndCompute(im1, None)

# 特征匹配
bf = cv2.BFMatcher(cv2.NORM_L1, crossCheck=True)
matches = bf.match(d0, d1)
matches = sorted(matches, key=lambda x: x.distance)

# 计算单应性矩阵
src_pts = array([l0[m.queryIdx].pt for m in matches])
dst_pts = array([l1[m.trainIdx].pt for m in matches])
H, mask = cv2.findHomography(src_pts, dst_pts, cv2.RANSAC, 5.0)

4.4 增强现实

增强现实（Augmented Reality，简称AR）是一种技术，它将数字信息（例如图像、音频、视频和3D模型）与现实世界中的环境相结合，创造出一个综合的、增强的视觉体验。AR技术通过使用传感器、摄像头、显示屏和计算设备，将虚拟内容叠加到真实世界中，使用户能够在真实环境中看到并与虚拟元素进行交互。下面就是利用PyGame和PyOpenGL库进行虚拟物体摆放的示例

4.4.1 从照相机矩阵到OpenGL格式

从照相机矩阵到OpenGL格式使原来的一个内标定矩阵变成了照相机和场景的GL_PROJECTION和GL_MODELVIEW两个矩阵，前者等价于内标定矩阵，后者处理物体与照相机之间的三维变换关系，下面是假设已知标定矩阵实现OpenGL中的投影矩阵代码

def set_projection_from_camera(K, width, height):
    """
    Set the OpenGL projection matrix based on the camera intrinsic matrix.

    Args:
        K (numpy.ndarray): The camera intrinsic matrix.
        width (int): The width of the viewport.
        height (int): The height of the viewport.
    """
    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)

下面函数实现移除内标定矩阵后的照相机矩阵并创建模拟视图，实际上这段代码就是根据相机的姿态信息设置OpenGL模型视图矩阵，以正确渲染虚拟物体（如茶壶模型）并实现增强现实效果

def set_modelview_from_camera(Rt):
    """
    Set the modelview matrix from camera pose.

    Args:
        Rt (numpy.ndarray): The camera pose matrix (3x4).
    """
    glMatrixMode(GL_MODELVIEW)
    glLoadIdentity()

    # Rotate the teapot 90 degrees around the x-axis to make the z-axis point upward
    Rx = np.array([[1, 0, 0], [0, 0, -1], [0, 1, 0]])

    # Get the best approximation of the rotation
    R = Rt[:, :3]
    U, S, V = np.linalg.svd(R)
    R = np.dot(U, V)
    R[0, :] = -R[0, :]  # Change the sign of the x-axis

    # Get the translation vector
    t = Rt[:, 3]

    # Get a 4x4 modelview matrix
    M = np.eye(4)
    M[:3, :3] = np.dot(R, Rx)
    M[:3, 3] = t

    # Transpose and flatten to get column-major values
    M = M.T
    m = M.flatten()

    # Replace the modelview matrix with the new matrix
    glLoadMatrixf(m)

4.4.2 在图像中放置虚拟物体

直接给出完整代码

from pylab import *
from OpenGL.GL import *
from OpenGL.GLU import *
from OpenGL.GLUT import *
import pygame, pygame.image
from pygame.locals import *
import cv2

from 视觉编程 import camera, homography


def cube_points(c, wid):  # 绘制立方体的一各点列表
    """ Creates a list of points for plotting
        a cube with plot. (the first 5 points are
        the bottom square, some sides repeated). """
    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


def my_calibration(sz):
    row, col = sz
    fx = 2555 * col / 2592
    fy = 2586 * row / 1936
    K = diag([fx, fy, 1])
    K[0, 2] = 0.5 * col
    K[1, 2] = 0.5 * row
    return K


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)


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


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)  # 清除纹理


def draw_teapot(size):  # 红色茶壶
    glEnable(GL_LIGHTING)
    glEnable(GL_LIGHT0)
    glEnable(GL_DEPTH_TEST)
    glClear(GL_DEPTH_BUFFER_BIT)
    # 绘制红色茶壶
    glMaterialfv(GL_FRONT, GL_AMBIENT, [0, 0, 0, 0])
    glMaterialfv(GL_FRONT, GL_DIFFUSE, [0.5, 0.0, 0.0, 0.0])
    glMaterialfv(GL_FRONT, GL_SPECULAR, [0.7, 0.6, 0.6, 0.0])
    glMaterialf(GL_FRONT, GL_SHININESS, 0.25 * 128.0)
    glutSolidTeapot(size)


def drawFunc(size):  # 白色茶壶
    glRotatef(0.5, 5, 5, 0)  # (角度,x,y,z)
    glutWireTeapot(size)
    # 刷新显示
    glFlush()


width, height = 4032, 3024


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


imname1 = './filelist/can4.jpeg'
imname2 = './filelist/can3.jpeg'

# compute features
sift = cv2.SIFT_create()

# Calculate features for the first image
im0 = cv2.imread(imname1)
l0, d0 = sift.detectAndCompute(im0, None)

# Calculate features for the second image
im1 = cv2.imread(imname2)
l1, d1 = sift.detectAndCompute(im1, None)

# Feature matching and homography estimation
bf = cv2.BFMatcher(cv2.NORM_L1, crossCheck=True)
matches = bf.match(d0, d1)
matches = sorted(matches, key=lambda x: x.distance)

src_pts = array([l0[m.queryIdx].pt for m in matches])
dst_pts = array([l1[m.trainIdx].pt for m in matches])
H, mask = cv2.findHomography(src_pts, dst_pts, cv2.RANSAC, 5.0)

# 计算照相机标定矩阵
K = my_calibration((3024, 4032))
# 位于边长为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))

Rt = dot(linalg.inv(K), cam2.P)

setup()
draw_background("./filelist/can4.bmp")
set_projection_from_camera(K)
set_modelview_from_camera(Rt)

draw_teapot(0.05)  # 显示红色茶壶
# drawFunc(0.05)  # 显示白色空心茶壶
pygame.display.flip()
while True:
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            sys.exit()

可能由于图像的位置以及特征的匹配度，这个茶壶摆的并不是很端正，但基本能够实现