reconstruction_mesh.py代码阅读

三角形、平面法线、顶点法线

在Direct3D中，三角形是构成实体的基本单位，因为一个三角形正好是一个平面，以三角形面为单位进行渲染效率最高。

一个三角形由三个点构成，习惯上把这些点称为顶点（Vertex）。三角形平面有正反面之分，由顶点的排序决定：顶点按顺时针排列的表面是正面，如图。

其中与三角形平面垂直、且指向正面的矢量称为该平面的法线（Normal）。
在Direct3D中，为提高渲染效率，缺省条件下只有正面可见，
顶点法线（Vertex Normal）是过顶点的一个矢量（法线是一个向量），用于在高洛德着色（Gouraud Shading）中的计算光照和纹理效果。在生成曲面时，通常令顶点法线和相邻平面的法线保持等角，如图1，这样进行渲染时，会在平面接缝处产生一种平滑过渡的效果。如果是多边形，则令顶点法线等于该点所属平面（三角形）的法线，如图2，以便在接缝处产生突出的边缘。

在opengl中为了模拟光线或进行光照计算和阴影计算，我们往往需要首先计算法线。表面的光照强度（即反射光量）是和光线方向与法线方向的夹角成正比的，夹角越小表面就会看起来越亮
为了使表面看起来更加光滑，有两种方法，第一种我们可以采用计算8个相邻顶点的法线，具体做法与计算四个顶点类似。第二种方法是对法线进行插值，具体是在着色器中完成，由于现在大多使用可编程管线，我们可以采用Phong光照模型，把顶点法线传到顶点着色器中，然后在片段着色器中对法线进行插值，使之看起来更光滑，但这种方法做出来的效果并不适合所有模型，看起来太鲜亮，有的模型需要偏暗色彩。

计算旋转矩阵

使用三个欧拉角（应该是roll / yaw / pitch）来计算旋转矩阵计算方法参考这里

欧拉角–>旋转矩阵

旋转矩阵–>欧拉角

投影层（将3D人脸投影到图像平面上）

根据经验选择焦距和相机位置（一般是固定下来的，Position是原点，look at direction是-z方向（右手坐标系下）,Up_direction朝向y轴）

光照层

如何看多维矩阵

比如维度为（1，2，3，4）

先看最后两个数，代表的是3行4列，就可以写出来（假设值都为1）

[[1,1,1,1],
[1,1,1,1],
[1,1,1,1]]

写完最后两个维度继续看第-3个维度是2，那么说明有2个3*4维，写成矩阵形式就是

[[[1,1,1,1],[1,1,1,1],[1,1,1,1]], [[1,1,1,1],[1,1,1,1],[1,1,1,1]]]

这时候继续看第-4维，例子中是1，那么直接再加上一个框即可，矩阵形式为

[[[[1,1,1,1],[1,1,1,1],[1,1,1,1]],[[1,1,1,1],[1,1,1,1],[1,1,1,1]]]]

看到矩阵怎么判断维度
比如如下：

[[[[1], [2]],
[[3], [4]]]]

首先，数最前面（或者最后面）边框，例子中共4个框说明共4维，把4维先写出来（None，None，None，None）
然后找单框，看到单框中只有一个1，所以是一列。然后在一个框中有两个单框，所以是两行。因此可知是两行一列，那么最后两维就是（2,1），总维度就变成（None，None，2, 1）。
接着找双框，发现只有一个逗号将两对双框隔开，所以第-3维是2，总维度变成（None，2，2, 1）。
最后找三框，发现只有一对三框并且没有发现有逗号隔开的三框，所以第一个维度是1，总维度为（1，2，2, 1）

实践一下

tensor([[[1.0150e+03, 0.0000e+00, 1.1200e+02],[0.0000e+00, 1.0150e+03, 1.1200e+02],[0.0000e+00, 0.0000e+00, 1.0000e+00]],[[1.0150e+03, 0.0000e+00, 1.1200e+02],[0.0000e+00, 1.0150e+03, 1.1200e+02],[0.0000e+00, 0.0000e+00, 1.0000e+00]]])
p_matrix = p_matrix.permute(0, 2, 1)

这里的话可以看到是将行和列进行交换，有点类似二位矩阵中的转置

tensor([[[1.0150e+03, 0.0000e+00, 0.0000e+00],[0.0000e+00, 1.0150e+03, 0.0000e+00],[1.1200e+02, 1.1200e+02, 1.0000e+00]],[[1.0150e+03, 0.0000e+00, 0.0000e+00],[0.0000e+00, 1.0150e+03, 0.0000e+00],[1.1200e+02, 1.1200e+02, 1.0000e+00]]])

SH 球谐函数

BRDF （双向反射分布函数）

代码

import torch
import math
import numpy as np
from utils import LeastSquaresdef split_coeff(coeff):# input: coeff with shape [1,257]id_coeff = coeff[:, :80]  # identity(shape) coeff of dim 80ex_coeff = coeff[:, 80:144]  # expression coeff of dim 64tex_coeff = coeff[:, 144:224]  # texture(albedo) coeff of dim 80angles = coeff[:, 224:227]  # ruler angles(x,y,z) for rotation of dim 3# lighting coeff for 3 channel SH function of dim 27gamma = coeff[:, 227:254]translation = coeff[:, 254:]  # translation coeff of dim 3return id_coeff, ex_coeff, tex_coeff, angles, gamma, translationclass _need_const:a0 = np.pia1 = 2 * np.pi / np.sqrt(3.0)a2 = 2 * np.pi / np.sqrt(8.0)c0 = 1 / np.sqrt(4 * np.pi)c1 = np.sqrt(3.0) / np.sqrt(4 * np.pi)c2 = 3 * np.sqrt(5.0) / np.sqrt(12 * np.pi)d0 = 0.5 / np.sqrt(3.0)illu_consts = [a0, a1, a2, c0, c1, c2, d0]origin_size = 300target_size = 224camera_pos = 10.0#用的是3DMM的特征脸形成模型
def shape_formation(id_coeff, ex_coeff, facemodel):# compute face shape with identity and expression coeff, based on BFM model# input: id_coeff with shape [1,80]#         ex_coeff with shape [1,64]# output: face_shape with shape [1,N,3], N is number of vertices'''S = mean_shape + \alpha * B_id + \beta * B_exp'''n_b = id_coeff.size(0)face_shape = torch.einsum('ij,aj->ai', facemodel.idBase, id_coeff) + \torch.einsum('ij,aj->ai', facemodel.exBase, ex_coeff) + \facemodel.meanshapeface_shape = face_shape.view(n_b, -1, 3)# re-center face shapeface_shape = face_shape - \facemodel.meanshape.view(1, -1, 3).mean(dim=1, keepdim=True)return face_shapedef texture_formation(tex_coeff, facemodel):# compute vertex texture(albedo) with tex_coeff# input: tex_coeff with shape [1,N,3]# output: face_texture with shape [1,N,3], RGB order, range from 0-255'''T = mean_texture + \gamma * B_texture'''n_b = tex_coeff.size(0)face_texture = torch.einsum('ij,aj->ai', facemodel.texBase, tex_coeff) + facemodel.meantexface_texture = face_texture.view(n_b, -1, 3)return face_texturedef compute_norm(face_shape, facemodel):# compute vertex normal using one-ring neighborhood (8 points)# input: face_shape with shape [1,N,3]# output: v_norm with shape [1,N,3]# https://fredriksalomonsson.files.wordpress.com/2010/10/mesh-data-structuresv2.pdf# vertex index for each triangle face, with shape [F,3], F is number of facesface_id = facemodel.tri - 1  #这里减去1是因为坐标从0开始# adjacent face index for each vertex, with shape [N,8], N is number of vertexpoint_id = facemodel.point_buf - 1shape = face_shapev1 = shape[:, face_id[:, 0], :]v2 = shape[:, face_id[:, 1], :]v3 = shape[:, face_id[:, 2], :]e1 = v1 - v2e2 = v2 - v3face_norm = e1.cross(e2)  # compute normal for each face 可是感觉这里指计算了一个面的法线呀。按理说要计算8个面的法线然后相加最后得到顶点法线# normalized face_norm firstface_norm = torch.nn.functional.normalize(face_norm, p=2, dim=2)empty = torch.zeros((face_norm.size(0), 1, 3),dtype=face_norm.dtype, device=face_norm.device)# concat face_normal with a zero vector at the endface_norm = torch.cat((face_norm, empty), 1)# compute vertex normal using one-ring neighborhoodv_norm = face_norm[:, point_id, :].sum(dim=2)v_norm = torch.nn.functional.normalize(v_norm, p=2, dim=2)  # normalize normal vectorsreturn v_normdef compute_rotation_matrix(angles):# compute rotation matrix based on 3 ruler angles# input: angles with shape [1,3]# output: rotation matrix with shape [1,3,3]n_b = angles.size(0)# https://www.cnblogs.com/larry-xia/p/11926121.htmldevice = angles.device# compute rotation matrix for X-axis, Y-axis, Z-axis respectivelyrotation_X = torch.cat([torch.ones([n_b, 1]).to(device),torch.zeros([n_b, 3]).to(device),torch.reshape(torch.cos(angles[:, 0]), [n_b, 1]),- torch.reshape(torch.sin(angles[:, 0]), [n_b, 1]),torch.zeros([n_b, 1]).to(device),torch.reshape(torch.sin(angles[:, 0]), [n_b, 1]),torch.reshape(torch.cos(angles[:, 0]), [n_b, 1])],axis=1)rotation_Y = torch.cat([torch.reshape(torch.cos(angles[:, 1]), [n_b, 1]),torch.zeros([n_b, 1]).to(device),torch.reshape(torch.sin(angles[:, 1]), [n_b, 1]),torch.zeros([n_b, 1]).to(device),torch.ones([n_b, 1]).to(device),torch.zeros([n_b, 1]).to(device),- torch.reshape(torch.sin(angles[:, 1]), [n_b, 1]),torch.zeros([n_b, 1]).to(device),torch.reshape(torch.cos(angles[:, 1]), [n_b, 1]),],axis=1)rotation_Z = torch.cat([torch.reshape(torch.cos(angles[:, 2]), [n_b, 1]),- torch.reshape(torch.sin(angles[:, 2]), [n_b, 1]),torch.zeros([n_b, 1]).to(device),torch.reshape(torch.sin(angles[:, 2]), [n_b, 1]),torch.reshape(torch.cos(angles[:, 2]), [n_b, 1]),torch.zeros([n_b, 3]).to(device),torch.ones([n_b, 1]).to(device),],axis=1)rotation_X = rotation_X.reshape([n_b, 3, 3])rotation_Y = rotation_Y.reshape([n_b, 3, 3])rotation_Z = rotation_Z.reshape([n_b, 3, 3])# R = Rz*Ry*Rxrotation = rotation_Z.bmm(rotation_Y).bmm(rotation_X)# because our face shape is N*3, so compute the transpose of R, so that rotation shapes can be calculated as face_shape*Rrotation = rotation.permute(0, 2, 1)return rotationdef projection_layer(face_shape, fx=1015.0, fy=1015.0, px=112.0, py=112.0):# we choose the focal length and camera position empirically# project 3D face onto image plane# input: face_shape with shape [1,N,3]#          rotation with shape [1,3,3]#         translation with shape [1,3]# output: face_projection with shape [1,N,2]#           z_buffer with shape [1,N,1]cam_pos = 10p_matrix = np.concatenate([[fx], [0.0], [px], [0.0], [fy], [py], [0.0], [0.0], [1.0]],axis=0).astype(np.float32)  # projection matrixp_matrix = np.reshape(p_matrix, [1, 3, 3])p_matrix = torch.from_numpy(p_matrix)gpu_p_matrix = Nonen_b, nV, _ = face_shape.size()if face_shape.is_cuda:gpu_p_matrix = p_matrix.cuda()p_matrix = gpu_p_matrix.expand(n_b, 3, 3)else:p_matrix = p_matrix.expand(n_b, 3, 3)face_shape[:, :, 2] = cam_pos - face_shape[:, :, 2]aug_projection = face_shape.bmm(p_matrix.permute(0, 2, 1))face_projection = aug_projection[:, :, 0:2] / aug_projection[:, :, 2:]z_buffer = cam_pos - aug_projection[:, :, 2:]return face_projection, z_bufferdef illumination_layer(face_texture, norm, gamma):# CHJ: It's different from what I knew.# compute vertex color using face_texture and SH function lighting approximation# input: face_texture with shape [1,N,3]#          norm with shape [1,N,3]#         gamma with shape [1,27]# output: face_color with shape [1,N,3], RGB order, range from 0-255#          lighting with shape [1,N,3], color under uniform texturen_b, num_vertex, _ = face_texture.size()n_v_full = n_b * num_vertexgamma = gamma.view(-1, 3, 9).clone()gamma[:, :, 0] += 0.8gamma = gamma.permute(0, 2, 1)a0, a1, a2, c0, c1, c2, d0 = _need_const.illu_constsY0 = torch.ones(n_v_full).float() * a0*c0if gamma.is_cuda:Y0 = Y0.cuda()norm = norm.view(-1, 3)nx, ny, nz = norm[:, 0], norm[:, 1], norm[:, 2]arrH = []arrH.append(Y0)arrH.append(-a1*c1*ny)arrH.append(a1*c1*nz)arrH.append(-a1*c1*nx)arrH.append(a2*c2*nx*ny)arrH.append(-a2*c2*ny*nz)arrH.append(a2*c2*d0*(3*nz.pow(2)-1))arrH.append(-a2*c2*nx*nz)arrH.append(a2*c2*0.5*(nx.pow(2)-ny.pow(2)))H = torch.stack(arrH, 1)Y = H.view(n_b, num_vertex, 9)# Y shape:[batch,N,9].# shape:[batch,N,3]lighting = Y.bmm(gamma)face_color = face_texture * lightingreturn face_color, lightingdef rigid_transform(face_shape, rotation, translation):n_b = face_shape.shape[0]face_shape_r = face_shape.bmm(rotation)  # R has been transposedface_shape_t = face_shape_r + translation.view(n_b, 1, 3)return face_shape_tdef compute_landmarks(face_shape, facemodel):# compute 3D landmark postitions with pre-computed 3D face shapekeypoints_idx = facemodel.keypoints - 1face_landmarks = face_shape[:, keypoints_idx, :]return face_landmarksdef compute_3d_landmarks(face_shape, facemodel, angles, translation):rotation = compute_rotation_matrix(angles)face_shape_t = rigid_transform(face_shape, rotation, translation)landmarks_3d = compute_landmarks(face_shape_t, facemodel)return landmarks_3ddef transform_face_shape(face_shape, angles, translation):rotation = compute_rotation_matrix(angles)face_shape_t = rigid_transform(face_shape, rotation, translation)return face_shape_tdef render_img(face_shape, face_color, facemodel, image_size=224, fx=1015.0, fy=1015.0, px=112.0, py=112.0, device='cuda:0'):'''ref: https://github.com/facebookresearch/pytorch3d/issues/184The rendering function (just for test)Input:face_shape:  Tensor[1, 35709, 3]face_color: Tensor[1, 35709, 3] in [0, 1]facemodel: contains `tri` (triangles[70789, 3], index start from 1)'''from pytorch3d.structures import Meshesfrom pytorch3d.renderer.mesh.textures import TexturesVertexfrom pytorch3d.renderer import (PerspectiveCameras,PointLights,RasterizationSettings,MeshRenderer,MeshRasterizer,SoftPhongShader,BlendParams)face_color = TexturesVertex(verts_features=face_color.to(device))face_buf = torch.from_numpy(facemodel.tri - 1)  # index start from 1face_idx = face_buf.unsqueeze(0)mesh = Meshes(face_shape.to(device), face_idx.to(device), face_color)R = torch.eye(3).view(1, 3, 3).to(device)R[0, 0, 0] *= -1.0T = torch.zeros([1, 3]).to(device)half_size = (image_size - 1.0) / 2focal_length = torch.tensor([fx / half_size, fy / half_size], dtype=torch.float32).reshape(1, 2).to(device)principal_point = torch.tensor([(half_size - px) / half_size, (py - half_size) / half_size], dtype=torch.float32).reshape(1, 2).to(device)cameras = PerspectiveCameras(device=device,R=R,T=T,focal_length=focal_length,principal_point=principal_point)raster_settings = RasterizationSettings(image_size=image_size,blur_radius=0.0,faces_per_pixel=1)lights = PointLights(device=device,ambient_color=((1.0, 1.0, 1.0),),diffuse_color=((0.0, 0.0, 0.0),),specular_color=((0.0, 0.0, 0.0),),location=((0.0, 0.0, 1e5),))blend_params = BlendParams(background_color=(0.0, 0.0, 0.0))renderer = MeshRenderer(rasterizer=MeshRasterizer(cameras=cameras,raster_settings=raster_settings),shader=SoftPhongShader(device=device,cameras=cameras,lights=lights,blend_params=blend_params))images = renderer(mesh)images = torch.clamp(images, 0.0, 1.0)return imagesdef estimate_intrinsic(landmarks_2d, transform_params, z_buffer, face_shape, facemodel, angles, translation):# estimate intrinsic parametersdef re_convert(landmarks_2d, trans_params, origin_size=_need_const.origin_size, target_size=_need_const.target_size):# convert landmarks to un_cropped imagesw = (origin_size * trans_params[2]).astype(np.int32)h = (origin_size * trans_params[2]).astype(np.int32)landmarks_2d[:, :, 1] = target_size - 1 - landmarks_2d[:, :, 1]landmarks_2d[:, :, 0] = landmarks_2d[:, :, 0] + w / 2 - target_size / 2landmarks_2d[:, :, 1] = landmarks_2d[:, :, 1] + h / 2 - target_size / 2landmarks_2d = landmarks_2d / trans_params[2]landmarks_2d[:, :, 0] = landmarks_2d[:, :, 0] + trans_params[3] - origin_size / 2landmarks_2d[:, :, 1] = landmarks_2d[:, :, 1] + trans_params[4] - origin_size / 2landmarks_2d[:, :, 1] = origin_size - 1 - landmarks_2d[:, :, 1]return landmarks_2ddef POS(xp, x):# calculating least sqaures problem# ref https://github.com/pytorch/pytorch/issues/27036ls = LeastSquares()npts = xp.shape[1]A = torch.zeros([2*npts, 4]).to(x.device)A[0:2*npts-1:2, 0:2] = x[0, :, [0, 2]]A[1:2*npts:2, 2:4] = x[0, :, [1, 2]]b = torch.reshape(xp[0], [2*npts, 1])k = ls.lstq(A, b, 0.010)fx = k[0, 0]px = k[1, 0]fy = k[2, 0]py = k[3, 0]return fx, px, fy, py# convert landmarks to un_cropped imageslandmarks_2d = re_convert(landmarks_2d, transform_params)landmarks_2d[:, :, 1] = _need_const.origin_size - 1.0 - landmarks_2d[:, :, 1]landmarks_2d[:, :, :2] = landmarks_2d[:, :, :2] * (_need_const.camera_pos - z_buffer[:, :, :])# compute 3d landmarkslandmarks_3d = compute_3d_landmarks(face_shape, facemodel, angles, translation)# compute fx, fy, px, pylandmarks_3d_ = landmarks_3d.clone()landmarks_3d_[:, :, 2] = _need_const.camera_pos - landmarks_3d_[:, :, 2]fx, px, fy, py = POS(landmarks_2d, landmarks_3d_)return fx, px, fy, pydef reconstruction(coeff, facemodel):# The image size is 224 * 224# face reconstruction with coeff and BFM modelid_coeff, ex_coeff, tex_coeff, angles, gamma, translation = split_coeff(coeff)# compute face shapeface_shape = shape_formation(id_coeff, ex_coeff, facemodel)# compute vertex texture(albedo)face_texture = texture_formation(tex_coeff, facemodel)# vertex normalface_norm = compute_norm(face_shape, facemodel)# rotation matrixrotation = compute_rotation_matrix(angles)face_norm_r = face_norm.bmm(rotation)# print(face_norm_r[:, :3, :])# do rigid transformation for face shape using predicted rotation and translationface_shape_t = rigid_transform(face_shape, rotation, translation)# compute 2d landmark projectionface_landmark_t = compute_landmarks(face_shape_t, facemodel)# compute 68 landmark on image plane (with image sized 224*224)landmarks_2d, z_buffer = projection_layer(face_landmark_t)landmarks_2d[:, :, 1] = _need_const.target_size - 1.0 - landmarks_2d[:, :, 1]# compute vertex color using SH function lighting approximationface_color, lighting = illumination_layer(face_texture, face_norm_r, gamma)return face_shape, face_texture, face_color, landmarks_2d, z_buffer, angles, translation, gamma