Learn Opencv ---- 使用Opencv构建哈哈镜

阅读原文

前言

    在童年时期，有一个东西时常能给我带来欢乐，那就是哈哈镜。哈哈镜不是传统的平面镜，而是凹凸不平的镜面。当人们来到镜子面前时，镜子中的自己会发生奇怪的扭曲，带来欢笑声。
    这篇博文的主要动机是鼓励读者去学习相机的内外参数，相机的投影矩阵以及如何使用几何来表现图像。
    通过这篇文章，我希望读者能够了解到一个事实，那就是如果想要创建一个有趣的东西，必须要有相对应的明确的概念和理论。

形成图像的理论

    简单来说，一张图像其实是一个三维世界场景的二维投影。而要将无数个三维点进行投影，需要使用下述公式：



    其中，$P$为相机的投影矩阵，$X_w,Y_w,Z_w$为三维空间的点坐标，$u,v$为二维空间的像素坐标。

如何实现

    整个哈哈镜项目主要分为三个步骤：

1. 创建虚拟摄像机
2. 定义一个三维曲面（镜面），并使用合适的投影矩阵将该曲面投影到虚拟相机当中。
3. 将三维曲面投影点的图像坐标应用于图像网格的扭曲，以获得哈哈镜的效果。

    整个步骤如下图所示。

图1: 创建哈哈镜的步骤图。第一步，创建三维曲面（左图）；第二步，使用虚拟相机捕获平面，得到对应的二维点（中图）；第三步，将得到的二维点应用到网格变形，得到最终的效果（右图）。
    如果大伙暂时还未理解，别担心，后面会详细解释每一步。

创建虚拟相机

    从上述的介绍可以知道一个三维点是如何使用投影矩阵然后得到二维坐标点的。现在让我们了解一下虚拟相机的含义，以及如何使用这个虚拟相机来拍摄图像。
    虚拟相机的本质是投影矩阵$P$，因为它告诉我们3D世界坐标和响应的图像像素坐标之间的关系。我们知道，一个投影矩阵是由相机内参矩阵（$K$）和相机外参（$M_1$）矩阵组成。

import numpy as np

# Tx,Ty,Tz表示虚拟相机在世界坐标系中的位置
# 定义平移矩阵
T = np.array([[1,0,0,-Tx],[0,1,0,-Ty],[0,0,1,-Tz]])

# alpha,beta,gamma分别为虚拟相机的偏离方向
# 定义旋转矩阵
# x轴
Rx = np.array([[1, 0, 0], [0, math.cos(alpha), -math.sin(alpha)], [0, math.sin(alpha), math.cos(alpha)]])
# y轴
Ry = np.array([[math.cos(beta), 0, -math.sin(beta)],[0, 1, 0],[math.sin(beta),0,math.cos(beta)]])
# z轴
Rz = np.array([[math.cos(gamma), -math.sin(gamma), 0],[math.sin(gamma),math.cos(gamma), 0],[0, 0, 1]])
# 构建成最终的旋转矩阵
R = np.matmul(Rx, np.matmul(Ry, Rz))

# 计算相机外参
M1 = np.matmul(R,T)

# 构建相机内参矩阵
# sx和sy是x和y坐标的缩放
# ox和oy是相机光学中心的坐标
K = np.array([[-focus/sx,sh,ox],[0,focus/sy,oy],[0,0,1]])

P = np.matmul(K,RT)

有了投影矩阵$P$后，我们该如何得到图像呢？
    首先，我们假设原始图像或视频帧是三维平面。当然，我们知道实际的三维场景通常都不是平面，但这影响不大，所以可以假设场景是平面的。记住，我们的目标不是无限的模拟一个真实的哈哈镜，因此可以容忍误差。
    一旦我们将所面对的三维场景视为平面，那我们就可以简单地将世界坐标与投影矩阵相乘，得到像素坐标$(u,v)$。（由于我们不考虑渲染问题，因此对于投射的光线，没做其他额外的处理）。
    总的来说，我们需要做的是捕捉（投影），首先，先将原始图像（视频帧）表示为你虚拟相机中的三维平面，然后使用投影矩阵将该平面上的每个点投影到虚拟相机的成像平面上。

定义三维曲面

    对于平面镜来说，Z是一个固定的常数。而对于非平面镜来说，先构建一个关于X和Y的网格，然后在其中更新Z的值，使其发生起伏。

图2：简单的三维曲面

    下面，给出如何构建一个三维曲面，并投影到相机平面的代码。

# 获取图像的大小
H,W = image.shape[:2]
# 定义x和y的坐标。这里是让图像的中心为零点
x = np.linspace(-W/2, W/2, W)
y = np.linspace(-H/2, H/2, H)

# 使用x和y的坐标构建网格
xv,yv = np.meshgrid(x,y)

# 扩充原有的二维网格到三维
# 设置Z为常数1，即平面
X = xv.reshape(-1,1)
Y = yv.reshape(-1,1)
Z = X*0+1
# 构建成三维世界的坐标
pts3d = np.concatenate(([X],[Y],[Z],[X*0+1]))[:,:,0]
# 使用投影矩阵将三维坐标投影到二维
pts2d = np.matmul(P,pts3d)
# 计算得到像素坐标
u = pts2d[0,:]/(pts2d[2,:]+0.00000001)
v = pts2d[1,:]/(pts2d[2,:]+0.00000001)

使用Python的VCAM库

    有没有觉得上述的代码很麻烦，尤其是每次都要重复的去写。而且当我们想要动态调整一些参数时，也不好进行修改。放心好啦，大神们已经为你铺好了路。vcam库的目的在于简化创建此类三维曲面，定义虚拟摄像机，设置所有参数以及进行投影任务。
    话不多说，直接开搞。首先当然是pip先行

pip install vcam

装好vcam库之后，我们来看看他的使用。该库的原理与我们上述的内容相同，但更为但方便和简洁。

import cv2
import numpy as np
import math
from vcam import vcam,meshGen

# 创建一个虚拟相机，并设定输入图像的大小
c1 = vcam(H=H,W=W)

# 创建一个与输入图像大小相同的网格
plane = meshGen(H,W)

# 修改Z的值，默认为1，即平面
# 将每个3D点的Z坐标定义为Z = 10*sin(2*pi[x/w]*10)
plane.Z = 10*np.sin((plane.X/plane.W)*2*np.pi*10)

# 获取得到最终的三维曲面
pts3d = plane.getPlane()

# 将三维曲面投影到二维图像坐标
pts2d = c1.project(pts3d)

现在得到的投影后的二维点可以用于基于网格的形变(重新映射)

图像重映射

    图像重映射说得是，将输入图像的每个像素经过重映射函数的运算后，来到新的位置，进而得到一张新的图像。这一过程可以使用数学公式来表现：

这样的方法被称为前向重映射或者前向变形。其中$ma{p}_x$和$ma{p}_y$为每个像素点$(x,y)$提供了新的位置。
    如果通过$ma{p}_x$和$ma{p}_y$映射和得不到有效的整数值呢？我们可以根据最接近的整数值将$(x,y)$处的像素强度值扩散到相邻的像素。但是这样做，会在新得到的图像中出现孔洞。该如何避开这些孔洞呢？
    答案是，使用反向重映射（反向变形）。意思就是，$ma{p}_x$和$ma{p}_y$将为我们提供旧图像的像素位置。它的数学公式如下：

    好了，我们现在只需要获取得到两个映射函数，即可完成有趣的图像变形啦。可是这两个函数该如何得到呢？别担心，Opencv为我们解决了一切！

# 使用投影得到的二维点集构建映射函数
# 这里的二维点集使用三维曲面投影得到
map_x,map_y = c1.getMaps(pts2d)

# 将两个映射函数作用与图像，得到最终图像
output = cv2.remap(img,map_x,map_y,interpolation=cv2.INTER_LINEAR)

cv2.imshow("Funny mirror",output)
cv2.waitKey(0)

图3：左图为输入图像，右图为经过变化后的图像。

总结

到这里基本上就结束啦，接下啦只要改变不同的Z的方式，就可以得到不同形变的图像了。更多的带注释的代码，可以去这里免费下载。