Face3D学习笔记（2）pipeline示例源码解析

写在前面

第一部分的Pipeline程序总体流程为将标准平均模型放到相机空间得到一张图片，利用这张图片来对当前模型进行变换。此程序目的为让读者了解Face3D的基本使用，为以后的重建程序编写打好基础。

• 为了保证整个示例项目更加直观，方便理解，在展示一些函数的源码时会使用numpy版本进行展示，而在示例程序中并未使用numpy版本的库，在Cython版本与numpy版本出现差异的原码前会有标注，希望读者留意。
• pipeline实例程序的jupyter版本在这里https://download.csdn.net/download/qq_45912037/85031147
  ==完全免费==，欢迎大家下载

正文

话不多说，直接开始

0.把该引的库引一引

import os, sys
import numpy as np
import scipy.io as sio
from skimage import io
from time import time
import matplotlib.pyplot as plt
sys.path.append('..')
import face3d
from face3d import mesh

1.载入网格数据

网格数据包括: 顶点, 三角网格数据, 颜色（可选）,纹理(可选)
这里用颜色来表示面部纹理

C = sio.loadmat('Data/example1.mat')
vertices = C['vertices']; colors = C['colors']; triangles = C['triangles']
colors = colors/np.max(colors)

这里的网格数据为.mat文件，分别取得其中的vertices、colors、triangles数据，
且colors数据进行了归一化。
数据规格如下：

2.改变顶点位置

改变网格对象在世界坐标中的位置

s = 180/(np.max(vertices[:,1]) - np.min(vertices[:,1]))

R = mesh.transform.angle2matrix([0, 30, 0]) 
t = [0, 0, 0]
transformed_vertices = mesh.transform.similarity_transform(vertices, s, R, t)

其中mesh.transform.angle2matrix源码如下

def angle2matrix(angles):
    ''' get rotation matrix from three rotation angles(degree). right-handed.
    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.
    '''
    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)

挺好理解的，是根据输入的角度值生成对应的旋转矩阵

另外的mesh.transform.similarity_transform源码如下

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

输入为网格顶点vertices、缩放比例s、旋转矩阵R和平移向量t3d
其中s = 180/(np.max(vertices[:,1]) - np.min(vertices[:,1]))
就是要缩放到纵向距离为180
执行空间坐标变换s*R.dot(X) + t后输出变换后的顶点位置

3. 修改颜色/纹理（添加光线）

加入点光源，光线位置在世界坐标系定义。

light_positions = np.array([[-128, -128, 300]])
light_intensities = np.array([[1, 1, 1]])
lit_colors = mesh.light.add_light(transformed_vertices, triangles, colors, light_positions, light_intensities)

其中mesh.light.add_light的源码如下

def add_light(vertices, triangles, colors, light_positions = 0, light_intensities = 0):
    ''' Gouraud shading. add point lights.
    In 3d face, usually assume:
    1. The surface of face is Lambertian(reflect only the low frequencies of lighting)
    2. Lighting can be an arbitrary combination of point sources
    3. No specular (unless skin is oil, 23333)

    Ref: https://cs184.eecs.berkeley.edu/lecture/pipeline    
    Args:
        vertices: [nver, 3]
        triangles: [ntri, 3]
        light_positions: [nlight, 3] 
        light_intensities: [nlight, 3]
    Returns:
        lit_colors: [nver, 3]
    '''
    nver = vertices.shape[0]
    normals = get_normal(vertices, triangles) # [nver, 3]

    # ambient
    # La = ka*Ia

    # diffuse
    # Ld = kd*(I/r^2)max(0, nxl)
    direction_to_lights = vertices[np.newaxis, :, :] - light_positions[:, np.newaxis, :] # [nlight, nver, 3]
    direction_to_lights_n = np.sqrt(np.sum(direction_to_lights**2, axis = 2)) # [nlight, nver]
    direction_to_lights = direction_to_lights/direction_to_lights_n[:, :, np.newaxis]
    normals_dot_lights = normals[np.newaxis, :, :]*direction_to_lights # [nlight, nver, 3]
    normals_dot_lights = np.sum(normals_dot_lights, axis = 2) # [nlight, nver]
    diffuse_output = colors[np.newaxis, :, :]*normals_dot_lights[:, :, np.newaxis]*light_intensities[:, np.newaxis, :]
    diffuse_output = np.sum(diffuse_output, axis = 0) # [nver, 3]
    
    # specular
    # h = (v + l)/(|v + l|) bisector
    # Ls = ks*(I/r^2)max(0, nxh)^p
    # increasing p narrows the reflectionlob

    lit_colors = diffuse_output # only diffuse part here.
    lit_colors = np.minimum(np.maximum(lit_colors, 0), 1)
    return lit_colors

其中get_normal函数在源码中另有定义，这里不再赘述。

在3d face中，通常假设：
1.人脸表面是Lambertian（朗伯表面），只会反射低频率的光
2.光照可以是点光源的任意组合
3.无镜面反射
这些参考了https://cs184.eecs.berkeley.edu/lecture/pipeline
但是这个网站好像挂掉了。

输入的参数分别有顶点坐标、三角网格数据、光源位置、光线强度。
经过运算后输出加入光源的颜色数据。
（这部分运算并不是太了解，以后可能会专门出一下光线这方面的解读）

4. 修改顶点位置

将对象从世界空间坐标转换到相机空间，即观察者视角。
（如果使用标准相机则忽略此步）

camera_vertices = mesh.transform.lookat_camera(transformed_vertices, eye = [0, 0, 200], at = np.array([0, 0, 0]), up = None)
# -- project object from 3d world space into 2d image plane. orthographic or perspective projection
projected_vertices = mesh.transform.orthographic_project(camera_vertices)

其中mesh.transform.lookat_camera的源码如下：

def normalize(x):
    epsilon = 1e-12
    norm = np.sqrt(np.sum(x**2, axis = 0))
    norm = np.maximum(norm, epsilon)
    return x/norm
def lookat_camera(vertices, eye, at = None, up = None):
    """ 'look at' transformation: from world space to camera space
    standard camera space: 
        camera located at the origin. 
        looking down negative z-axis. 
        vertical vector is y-axis.
    Xcam = R(X - C)
    Homo: [[R, -RC], [0, 1]]
    Args:
      vertices: [nver, 3] 
      eye: [3,] the XYZ world space position of the camera.5
      at: [3,] a position along the center of the camera's gaze.
      up: [3,] up direction 
    Returns:
      transformed_vertices: [nver, 3]
    """
    if at is None:
      at = np.array([0, 0, 0], np.float32)
    if up is None:
      up = np.array([0, 1, 0], np.float32)

    eye = np.array(eye).astype(np.float32)
    at = np.array(at).astype(np.float32)
    z_aixs = -normalize(at - eye) # look forward
    x_aixs = normalize(np.cross(up, z_aixs)) # look right
    y_axis = np.cross(z_aixs, x_aixs) # look up

    R = np.stack((x_aixs, y_axis, z_aixs))#, axis = 0) # 3 x 3
    transformed_vertices = vertices - eye # translation
    transformed_vertices = transformed_vertices.dot(R.T) # rotation
    return transformed_vertices

标准相机空间为：
相机位于原点;向下看是负z轴;垂直向量为y轴。

输入分别为顶点坐标、摄像机的XYZ世界空间位置、沿着相机视线中心的位置[默认为(0，0，0)]、向上方向[默认为(0，1，0)]
根据输入计算出旋转矩阵R并通过Xcam=R(X-C)计算出新的顶点位置。

5. 转换为二维图像

设置图像的长和宽为256

h = w = 256
# change to image coords for rendering
image_vertices = mesh.transform.to_image(projected_vertices, h, w)
# render 
rendering =  mesh.render.render_colors(image_vertices, triangles, lit_colors, h, w)

mesh.transform.to_image部分的源码如下：

def to_image(vertices, h, w, is_perspective = False):
    ''' change vertices to image coord system
    3d system: XYZ, center(0, 0, 0)
    2d image: x(u), y(v). center(w/2, h/2), flip y-axis. 
    Args:
        vertices: [nver, 3]
        h: height of the rendering
        w : width of the rendering
    Returns:
        projected_vertices: [nver, 3]  
    '''
    image_vertices = vertices.copy()
    if is_perspective:
        # if perspective, the projected vertices are normalized to [-1, 1]. so change it to image size first.
        image_vertices[:,0] = image_vertices[:,0]*w/2
        image_vertices[:,1] = image_vertices[:,1]*h/2
    # 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

输入为顶点坐标、图片的长宽以及透视选项is_perspective（默认为FALSE）
通过计算得到二维的顶点坐标。

mesh.render.render_colors的源码如下：
==注意，此为numpy版本==

def render_colors(vertices, triangles, colors, h, w, c = 3):
    ''' render mesh with colors
    Args:
        vertices: [nver, 3]
        triangles: [ntri, 3] 
        colors: [nver, 3]
        h: height
        w: width    
    Returns:
        image: [h, w, c]. 
    '''
    assert vertices.shape[0] == colors.shape[0]
    
    # initial 
    image = np.zeros((h, w, c))
    depth_buffer = np.zeros([h, w]) - 999999.

    for i in range(triangles.shape[0]):
        tri = triangles[i, :] # 3 vertex indices

        # the inner bounding box
        umin = max(int(np.ceil(np.min(vertices[tri, 0]))), 0)
        umax = min(int(np.floor(np.max(vertices[tri, 0]))), w-1)

        vmin = max(int(np.ceil(np.min(vertices[tri, 1]))), 0)
        vmax = min(int(np.floor(np.max(vertices[tri, 1]))), h-1)

        if umax<umin or vmax<vmin:
            continue

        for u in range(umin, umax+1):
            for v in range(vmin, vmax+1):
                if not isPointInTri([u,v], vertices[tri, :2]): 
                    continue
                w0, w1, w2 = get_point_weight([u, v], vertices[tri, :2])
                point_depth = w0*vertices[tri[0], 2] + w1*vertices[tri[1], 2] + w2*vertices[tri[2], 2]

                if point_depth > depth_buffer[v, u]:
                    depth_buffer[v, u] = point_depth
                    image[v, u, :] = w0*colors[tri[0], :] + w1*colors[tri[1], :] + w2*colors[tri[2], :]
    return image

输入为顶点坐标、三角网格数据、网格颜色数据以及目标图片的长宽。
输出为目标的带纹理的二维图像数据。

6.图片保存

save_folder = 'results/pipeline'
if not os.path.exists(save_folder):
    os.mkdir(save_folder)
io.imsave('{}/rendering.jpg'.format(save_folder), rendering)

Lossy conversion from float32 to uint8. Range [0, 1]. Convert image to uint8 prior to saving to suppress this warning.

7.结果展示

生成的二维图片展示

plt.imshow(rendering)
plt.show()

相机空间的Mesh也可以用里面的库函数mesh.vis.plot_mesh展示

# ---- show mesh
mesh.vis.plot_mesh(camera_vertices, triangles)
plt.show()