详解相机模型：矩阵变化+代码解析

0. 写在前面

在3dmm中，重要的一步是对3d模型进行拍照。

这里引出问题：怎么给3d模型拍照？

下面解决这个拍照问题。

1. 相机模型

本节目标：

1. 理解针孔相机的模型，内参与径向畸变参数。
2. 理解一个空间点是如何投影到相机成像。

我们知道：
一张照片（二维）：由多个像素组成，每个像素记录了色彩或亮度的信息。
一个物体（三维世界）：物体反射或发出的光线，通过相机光心后，投影到相机的成像平面。

相机将三维世界中的坐标点（单位米）映射到二维图像平面（单位像素）的过程能够用一个几何模型进行描述。

这个模型有很多种，其中最简单的为针孔模型。事实上，真实的相机镜头是透镜，会使得光线投影到成像平面的过程会产生畸变。

那么总结一下，所谓相机模型实际是是：针孔相机模型+畸变模型。

在3dmm中我们暂不考虑畸变模型。

1.1 针孔相机模型 (Pinhole camera)

在小孔成像过程中，小孔模型将三维世界中的蜡烛投影到一个二维成像平面。

首先我们要先认识一下4种坐标：

1. 世界坐标 (World reference system)
2. 相机坐标 (camera coordinate)
3. 归一化相机坐标 (normalized coordinate)
4. 像素坐标 (pixel coordinate)

我们来品一下这张图><。

1. 世界坐标 是我们客观存在的世界，它有自己的固有坐标。在这里我们定义三维空间的三个方向分别为:$$x_w,y_w,z_w$$
   例子：
   例子1: 长城。它就存在在那里，有它自己的坐标。
   例子2: 3d模型。3d模型在计算机中是以点的坐标来存储的，这个坐标代表点在$x_w,y_w,z_w$三个方向上的大小，比如一个点$B$，它的存储形式是$(12,15,18)$。这样的点有$n$个，在3dmm中实际上有53215个这样的三维点，这些点组成了基本的人脸3d模型。

2. 相机坐标 是从相机的角度去看世界，相机本身是这个坐标系的原点。
   在这里我们定义在这个角度的三维空间的三个方向分别为:$$x_c,y_c,z_c$$
   例子：
   例子1: 小时候我总觉得门前的山特别高，后来长大之后回到老地方，发现这个山也没那么高嘛。山变了吗？山没变，是我看山的角度变了。
   例子2: 正所谓【横看成岭侧成峰】，说的也是这回事。

3. 归一化相机坐标 物理成像平面。在这里我们可以定义三维空间的三个方向分别为:$$x_n,y_n,z_n$$
   原点$\acute{O}$是相机坐标的$z_c$轴与物理成像平面的交点。一般取$z_n$的方向与$z_c$方向相同，将点$P$在针孔相机模型到物理成像平面的距离统一为$z_n=f$，其中$f$为焦距，这就相当于对相机坐标进行归一化。

4. 像素坐标 水平方向是U，垂直方向是V，通过这个平面的，二维的UV坐标系。我们可以定位图象上的任意一个像素。
   

注意 这里的像素坐标的单位是像素，而上述3个坐标系的单位为米。
在这里我们定义二维空间的二个方向分别为:$$u,v$$

1.2 计算过程

================== 计算过程（1） 从世界坐标到相机坐标 ==================

已知我们有一个点$P$，它的世界坐标（我们常说的坐标）为:
$$P_w=\left[\begin{array}{c}{x}_w\\ y_w\\ z_w\end{array}\right]$$

这个点$P$，它在相机视角下的坐标为：
$$P_c=\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\end{array}\right]$$

我们知道$P_w$和$P_c$是同一点在不同坐标系下的表达方式，那么$P_w$和$P_c$之间是什么关系呢？其实是：

$P_w$ $\underrightarrow{旋转变换+平移变换}$ $P_c$

这里需要引入一下齐次坐标的概念。

齐次坐标就是将一个原本是n维的向量用一个n+1维向量来表示，是指一个用于投影几何里的坐标系统，如同用于欧氏几何里的笛卡儿坐标一般。

比如从欧式坐标变换到齐次坐标时，即$欧式坐标\to 齐次坐标$：
$$(x,y)=\left[\begin{array}{c}x\\ y\\ z\end{array}\right]$$ $$(x,y,z)=\left[\begin{array}{c}x\\ y\\ z\\ 1\end{array}\right]$$
从齐次坐标变换到欧式坐标时，即$齐次坐标\to 欧式坐标$：

$$\left[\begin{array}{c}x\\ y\\ w\end{array}\right]=(\frac{x}{w},\frac{y}{w})$$ $$\left[\begin{array}{c}x\\ y\\ z\\ w\end{array}\right]=(\frac{x}{w},\frac{y}{w},\frac{z}{w})$$

接下来回到$P_w$ $\underrightarrow{旋转变换+平移变换}$ $P_c$。

(1)旋转变换

$R_x(\alpha )=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \alpha & -\sin \alpha \\ 0& \sin \alpha & \cos \alpha \end{array}\right]$
$R_y(\beta )=\left[\begin{array}{ccc}\cos \beta & 0& \sin \beta \\ 0& 1& 0\\ -\sin \beta & 0& \cos \beta \end{array}\right]$
$R_x(\gamma )=\left[\begin{array}{ccc}\cos \gamma & -\sin \gamma & 0\\ \sin \gamma & \cos \gamma & 0\\ 0& 0& 1\end{array}\right]$

此时：
$$R=R_x(\alpha )R_y(\beta )R_x(\gamma )$$

(2)平移变换
$$T=\left[\begin{array}{c}{T}_x\\ T_y\\ T_z\end{array}\right]$$

结合（1）（2），那么对于点$P$（下角标${}_h$表示齐次坐标）:
$$P_{hc}=\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\\ 1\end{array}\right]={\left[\begin{array}{cc}R& T\\ 0& 1\end{array}\right]}_{4\times 4}{\left[\begin{array}{c}{x}_w\\ y_w\\ z_w\\ 1\end{array}\right]}_{4\times 1}={\left[\begin{array}{cc}R& T\\ 0& 1\end{array}\right]}_{4\times 4}{P}_{hw}$$

================== 计算过程（2） 从相机坐标到像素坐标 ==================

由图可知：

$△O\acute{O}\acute{P}$与$△O\acute{P}\acute{\stackrel{ˊ}{P}}$相似，那么对于点$P$在相机坐标下的表示$P_c=\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\end{array}\right]$ 及在归一化相机坐标系里的表现$P_n=\left[\begin{array}{c}{x}_n\\ y_n\\ z_n\end{array}\right]$, （$z_n=f$）：
$$x_n=f\frac{x_c}{z_c}$$ $$y_n=f\frac{y_c}{z_c}$$

接下来要从归一化相机坐标系继续变换到像素坐标系，归一化相机坐标系与像素坐标系之间，相差了一个缩放和原点的平移。

我们设像素坐标在$u$轴上缩放了$\alpha$倍，在$v$轴上缩放了$\beta$倍。同时，原点平移了 $\left[\begin{array}{c}{c}_x\\ c_y\end{array}\right]$ ，那么像素坐标系下的$P_{uv}$与归一化相机坐标系下的$P_n$的关系为：
$$u=\alpha x_n+c_x=\alpha f\frac{x_c}{z_c}+c_x$$ $$v=\beta y_n+c_y=\beta f\frac{y_c}{z_c}+c_y$$
我们设$\alpha f=f_x$, $\beta f=f_y$,得到：
$$u=f_x\frac{x_c}{z_c}+c_x$$ $$v=f_y\frac{y_c}{z_c}+c_y$$

上述内容总结一下，用齐次坐标来表示：
$$P_{huv}=\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]=\left[\begin{array}{c}{f}_x\frac{x_c}{z_c}+c_x\\ f_y\frac{y_c}{z_c}+c_y\\ 1\end{array}\right]$$
因为齐次坐标乘以非零常数后表达相同含义：

$$\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]=\left[\begin{array}{c}{f}_x\frac{x_c}{z_c}+c_x\\ f_y\frac{y_c}{z_c}+c_y\\ 1\end{array}\right]=\left[\begin{array}{c}{f}_x{x}_c+c_x{z}_c\\ f_y{y}_c+c_y{z}_c\\ z_c\end{array}\right]$$ $$=\left[\begin{array}{ccc}{f}_x& 0& c_x\\ 0& f_y& c_y\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\end{array}\right]$$

回忆一下的（1）部分的结果：
$$P_{hc}=\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\\ 1\end{array}\right]={\left[\begin{array}{cc}R& T\\ 0& 1\end{array}\right]}_{4\times 4}{\left[\begin{array}{c}{x}_w\\ y_w\\ z_w\\ 1\end{array}\right]}_{4\times 1}={\left[\begin{array}{cc}R& T\\ 0& 1\end{array}\right]}_{4\times 4}{P}_{hw}$$

这里的$P_{ch}$是$P_c$齐次坐标，为了将这一部分带入$\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]$,
$\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]=\left[\begin{array}{c}{f}_x\frac{x_c}{z_c}+c_x\\ f_y\frac{y_c}{z_c}+c_y\\ 1\end{array}\right]$
$=\left[\begin{array}{c}{f}_x{x}_c+c_x{z}_c\\ f_y{y}_c+c_y{z}_c\\ z_c\end{array}\right]$
$=\left[\begin{array}{ccc}{f}_x& 0& c_x\\ 0& f_y& c_y\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\end{array}\right]$
$\uparrow$
我们称这个矩阵为内部参数$K$，

$={\left[\begin{array}{cccc}{f}_x& 0& c_x& 0\\ 0& f_y& c_y& 0\\ 0& 0& 1& 0\end{array}\right]}_{3\times 4}{\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\\ 1\end{array}\right]}_{4\times 1}$
$\uparrow \uparrow$
这可以写作$K\left[\begin{array}{cc}I& 0\end{array}\right]$ &ThinSpace;&ThinSpace;&ThinSpace;&ThinSpace;&ThinSpace;&ThinSpace;&ThinSpace;&ThinSpace;&ThinSpace;&ThinSpace;&ThinSpace;&ThinSpace;&ThinSpace;&ThinSpace;&ThinSpace; \,\,\,\,\,\, \,\,\,\,\,\,\,\,\, 是不是很眼熟，这不是$P_c$的齐次坐标嘛！

$={\left[\begin{array}{cccc}{f}_x& 0& c_x& 0\\ 0& f_y& c_y& 0\\ 0& 0& 1& 0\end{array}\right]}_{3\times 4}{\left[\begin{array}{cc}R& T\\ 0& 1\end{array}\right]}_{4\times 4}{\left[\begin{array}{c}{x}_w\\ y_w\\ z_w\\ 1\end{array}\right]}_{4\times 1}$

$=K\left[\begin{array}{cc}I& 0\end{array}\right]\left[\begin{array}{cc}R& T\\ 0& 1\end{array}\right]P_{hw}$

$=K\left[\begin{array}{cc}R& T\end{array}\right]P_{hw}$

============================== 总结 ==============================

$$\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]={\left[\begin{array}{cccc}{f}_x& 0& c_x& 0\\ 0& f_y& c_y& 0\\ 0& 0& 1& 0\end{array}\right]}_{3\times 4}{\left[\begin{array}{cc}R& T\\ 0& 1\end{array}\right]}_{4\times 4}{\left[\begin{array}{c}{x}_w\\ y_w\\ z_w\\ 1\end{array}\right]}_{4\times 1}$$

即：
$$P_{huv}=K\left[\begin{array}{cc}R& T\end{array}\right]P_{hw}$$

当我们有一个3d模型的n个点的三维坐标，可以通过这个变换得到它的照片。

2. 例子 + 代码

比如以我们提到的点$B$为例，它的世界坐标是$(12,15,18)$。

================== 计算过程（1） 从世界坐标到相机坐标 ==================

(1)旋转变换

$R_x(\alpha )=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \alpha & -\sin \alpha \\ 0& \sin \alpha & \cos \alpha \end{array}\right]$
$R_y(\beta )=\left[\begin{array}{ccc}\cos \beta & 0& \sin \beta \\ 0& 1& 0\\ -\sin \beta & 0& \cos \beta \end{array}\right]$
$R_x(\gamma )=\left[\begin{array}{ccc}\cos \gamma & -\sin \gamma & 0\\ \sin \gamma & \cos \gamma & 0\\ 0& 0& 1\end{array}\right]$

此时：
$$R=R_x(\alpha )R_y(\beta )R_x(\gamma )$$

def angle2matrix(angles):
    
    ''' 
    根据右手系三个旋转角，
    得到三个旋转矩阵。
    Args:
        angles: [3,]. x, y, z angles
        x: pitch. positive for looking down.
        y: yaw. positive for looking left. 
        z: roll. positive for tilting head right. 
    Returns:
        R: [3, 3]. rotation matrix.
    
    '''
    
    #np.deg2rad将角度变为弧度pi
    x, y, z = np.deg2rad(angles[0]), np.deg2rad(angles[1]), np.deg2rad(angles[2])
    # x
    Rx=np.array([[1,      0,       0],
                 [0, cos(x),  -sin(x)],
                 [0, sin(x),   cos(x)]])
    # y
    Ry=np.array([[ cos(y), 0, sin(y)],
                 [      0, 1,      0],
                 [-sin(y), 0, cos(y)]])
    # z
    Rz=np.array([[cos(z), -sin(z), 0],
                 [sin(z),  cos(z), 0],
                 [     0,       0, 1]])
    
    R=Rz.dot(Ry.dot(Rx))
    return R.astype(np.float32)

(2)平移变换
$$T=\left[\begin{array}{c}{T}_x\\ T_y\\ T_z\end{array}\right]$$

结合（1）（2），那么对于点$P$:
$$P_{hc}=\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\\ 1\end{array}\right]={\left[\begin{array}{cc}R& T\\ 0& 1\end{array}\right]}_{4\times 4}{\left[\begin{array}{c}{x}_w\\ y_w\\ z_w\\ 1\end{array}\right]}_{4\times 1}={\left[\begin{array}{cc}R& T\\ 0& 1\end{array}\right]}_{4\times 4}{P}_{hw}$$

def similarity_transform(vertices, s, R, t3d):
    ''' similarity transform. dof = 7.
    3D: s*R.dot(X) + t
    Homo: M = [[sR, t],[0^T, 1]].  M.dot(X)
    Args:(float32)
        vertices: [nver, 3]. 
        s: [1,]. scale factor.
        R: [3,3]. rotation matrix.
        t3d: [3,]. 3d translation vector.
    Returns:
        transformed vertices: [nver, 3]
    '''
    t3d = np.squeeze(np.array(t3d, dtype = np.float32))
    transformed_vertices = s * vertices.dot(R.T) + t3d[np.newaxis, :]

    return transformed_vertices

值得注意的是：函数angle2matrix(angles)中从世界坐标到相机坐标的计算并没有用到齐次坐标。只在欧式坐标下计算，
$$P_c=R_{3\times 3}\left[\begin{array}{c}{x}_w\\ y_w\\ z_w\end{array}\right]+\left[\begin{array}{c}{T}_x\\ T_y\\ T_z\end{array}\right]$$这并没有关系，因为齐次坐标和非齐次坐标本质上并没有区别:
$$P_{hc}=\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\\ 1\end{array}\right]P_c=\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\end{array}\right]$$只要在用的时候区分开就好。

我们之所以引入齐次坐标是为了直接计算$P_{huv}=K\left[\begin{array}{cc}R& T\end{array}\right]P_{hw}$。

================== 计算过程（2） 从相机坐标到像素坐标 ==================

我们已经知道：

$$\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]=\left[\begin{array}{ccc}{f}_x& 0& c_x\\ 0& f_y& c_y\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\end{array}\right]$$

还原一下这个过程，相机坐标$\to$归一化相机坐标$\to$像素坐标。

（2.1） 相机坐标$\to$归一化相机坐标
$$z_n=f$$ $$x_n=f\frac{x_c}{z_c}$$ $$y_n=f\frac{y_c}{z_c}$$

我们可以假设物品的深度，远远小于物体与相机间的距离，比如两个点$i,j$，两个点的$z_c$方向坐标$z_i\approx z_j$，此时我们可以直接删去$z_c$方向坐标，且归一化时我们取$f=1$:
$$z_n=f=1$$ $$x_n=x_c$$ $$y_n=y_c$$
这样的映射也称为平行映射(orthographic project).

def orthographic_project(vertices):
    return vertices.copy()

这里的代码虽然保留$z_c$，计算上其实我们只取$x_c,y_c$。

(2.2) 归一化相机坐标$\to$像素坐标

设归一化相机坐标$\to$像素坐标时无缩放， 归一化相机坐标的圆心$\acute{O}$落在像素为$(w,h)$的2d照片的中心，此时$\acute{O}$在像素坐标系下的坐标为 $$(w/2,h/2)$$

那么对于点$B$:

$$u=f_x\frac{x_c}{z_c}+c_x=x_c+w/2$$ $$v=f_y\frac{y_c}{z_c}+c_y=-y_c+h/2-1$$

def to_image(vertices, h, w):
    ''' 
    Args:
        vertices: [nver, 3]
        h: height of the rendering
        w : width of the rendering
    Returns:
        projected_vertices: [nver, 3]  
    '''
    image_vertices = vertices.copy()

    # move to center of image
    image_vertices[:,0] = image_vertices[:,0] + w/2
    image_vertices[:,1] = image_vertices[:,1] + h/2
    # flip vertices along y-axis.
    image_vertices[:,1] = h - image_vertices[:,1] - 1
    return image_vertices

至此，得到三维点$B$在照片里的像素坐标。