寫在前面

• 為了保證整個示例項目更加直觀，方便理解，在展示一些函數的源碼時會使用numpy版本進行展示，而在示例程序中并未使用numpy版本的庫，在Cython版本與numpy版本出現差異的原碼前會有標注，希望讀者留意。
• 3DMM實例程序的jupyter版本后續會更新，完全免費，歡迎大家下載

在Face3d中的求解過程可以概述如下：
(1)初始化$\alpha ,\beta$為0；
(2)利用黃金標準算法得到一個仿射矩陣$P_A$，分解得到$s,R,t_{2d}$；
(3)將(2)中求出的$s,R,t_{2d}$帶入能量方程，解得$\beta$；
(4)將(2)和(3)中求出的$\alpha$代入能量方程，解得$\alpha$；
(5)更新$\alpha ,\beta$的值，重復(2)-(4)進行迭代更新。

代碼解析

(3).將(2)中求出的$s,R,t_{2d}$帶入能量方程，解得$\beta$

在上一篇文章我們通過黃金標準算法求解出了仿射矩陣$P_A$并將它分解的到了$s,R,t_{2d}$，這部分將繼續按照求解步驟進行源碼分析:

for i in range(max_iter):X = shapeMU + shapePC.dot(sp) + expPC.dot(ep)X = np.reshape(X, [int(len(X)/3), 3]).T#----- estimate poseP = mesh.transform.estimate_affine_matrix_3d22d(X.T, x.T)s, R, t = mesh.transform.P2sRt(P)rx, ry, rz = mesh.transform.matrix2angle(R)# print('Iter:{}; estimated pose: s {}, rx {}, ry {}, rz {}, t1 {}, t2 {}'.format(i, s, rx, ry, rz, t[0], t[1]))#----- estimate shape# expressionshape = shapePC.dot(sp)shape = np.reshape(shape, [int(len(shape)/3), 3]).Tep = estimate_expression(x, shapeMU, expPC, model['expEV'][:n_ep,:], shape, s, R, t[:2], lamb = 0.002)# shapeexpression = expPC.dot(ep)expression = np.reshape(expression, [int(len(expression)/3), 3]).Tsp = estimate_shape(x, shapeMU, shapePC, model['shapeEV'][:n_sp,:], expression, s, R, t[:2], lamb = 0.004)return sp, ep, s, R, t

對于公式

其中的形狀部分為$\sum _{i=1}^m{a}_i{S}_i$，通過
shape = shapePC.dot(sp)
定義shape的值為$\sum _{i=1}^m{a}_i{S}_i$。

shape的格式為（159645，1），再通過
shape = np.reshape(shape, [int(len(shape)/3), 3]).T
將shape的XYZ坐標分開，轉為（53215，3）格式。

下段代碼
ep = estimate_expression(x, shapeMU, expPC, model['expEV'][:n_ep,:], shape, s, R, t[:2], lamb = 0.002)
其中ep即是$\beta$, estimate_expression的源碼如下：

def estimate_expression(x, shapeMU, expPC, expEV, shape, s, R, t2d, lamb = 2000):'''Args:x: (2, n). image points (to be fitted)shapeMU: (3n, 1)expPC: (3n, n_ep)expEV: (n_ep, 1)shape: (3, n)s: scaleR: (3, 3). rotation matrixt2d: (2,). 2d translationlambda: regulation coefficientReturns:exp_para: (n_ep, 1) shape parameters(coefficients)'''x = x.copy()assert(shapeMU.shape[0] == expPC.shape[0])assert(shapeMU.shape[0] == x.shape[1]*3)dof = expPC.shape[1]n = x.shape[1]sigma = expEVt2d = np.array(t2d)P = np.array([[1, 0, 0], [0, 1, 0]], dtype = np.float32)A = s*P.dot(R) #(2,3)# --- calc pcpc_3d = np.resize(expPC.T, [dof, n, 3]) pc_3d = np.reshape(pc_3d, [dof*n, 3]) # (29n,3)pc_2d = pc_3d.dot(A.T) #(29n,2)pc = np.reshape(pc_2d, [dof, -1]).T # 2n x 29# --- calc b# shapeMUmu_3d = np.resize(shapeMU, [n, 3]).T # 3 x n# expressionshape_3d = shape# b = A.dot(mu_3d + shape_3d) + np.tile(t2d[:, np.newaxis], [1, n]) # 2 x nb = np.reshape(b.T, [-1, 1]) # 2n x 1# --- solveequation_left = np.dot(pc.T, pc) + lamb * np.diagflat(1/sigma**2)x = np.reshape(x.T, [-1, 1])equation_right = np.dot(pc.T, x - b)exp_para = np.dot(np.linalg.inv(equation_left), equation_right)return exp_para

數據處理

x = x.copy()assert(shapeMU.shape[0] == expPC.shape[0])assert(shapeMU.shape[0] == x.shape[1]*3)dof = expPC.shape[1]n = x.shape[1]sigma = expEVt2d = np.array(t2d)P = np.array([[1, 0, 0], [0, 1, 0]], dtype = np.float32)

首先是確認輸入的格式正確：
assert(shapeMU.shape[0] == expPC.shape[0])
assert(shapeMU.shape[0] == x.shape[1]*3)
然后此時輸入的表情主成分expPC格式為(159645,29)
令dof=29
dof = expPC.shape[1]
令n=68
n = x.shape[1]
另sigma=expEV 即表情主成分方差$\sigma$
sigma = expEV
$t_{2d}$轉化為array數組
t2d = np.array(t2d)
P即正交投影矩陣$P_{orth}=\left[\begin{array}{lcc}1& 0& 0\\ 0& 1& 0\end{array}\right]$
P = np.array([[1, 0, 0], [0, 1, 0]], dtype = np.float32)

A = s*P.dot(R) #(2,3)# --- calc pcpc_3d = np.resize(expPC.T, [dof, n, 3]) pc_3d = np.reshape(pc_3d, [dof*n, 3]) # (29n,3)pc_2d = pc_3d.dot(A.T) #(29n,2)pc = np.reshape(pc_2d, [dof, -1]).T # 2n x 29# --- calc b# shapeMUmu_3d = np.resize(shapeMU, [n, 3]).T # 3 x n# expressionshape_3d = shape# b = A.dot(mu_3d + shape_3d) + np.tile(t2d[:, np.newaxis], [1, n]) # 2 x nb = np.reshape(b.T, [-1, 1]) # 2n x 1

已知公式

定義和計算A、pc和b：

• 定義A:

A = s*P.dot(R)

• 計算pc：

這里的pc計算相當于下式

將表情主成分expPC轉換為（29，68，3）的新矩陣pc_3d
pc_3d = np.resize(expPC.T, [dof, n, 3])

注意：這里的expPC經過了
expPC = model['expPC'][valid_ind, :n_ep]
運算，只包含特征點的表情主成分，格式為(68*3,29)

將pc_3d轉換為(29*68,3)的格式：
pc_3d = np.reshape(pc_3d, [dof*n, 3])

計算出$p{c}_{2d}$=$p{c}_{3d}?A^T$，pc_2d格式為（29*68，2）：
pc_2d = pc_3d.dot(A.T)

得出將pc_2d展開后得pc:
pc = np.reshape(pc_2d, [dof, -1]).T

• 定義b

b的公式如下：

計算時由于矩陣的格式問題要先進行一些變換：
這里的shapeMU也是只包含了68個特征點的
將格式為(68*3,1)的shapeMU轉換為格式（3，68）：
mu_3d = np.resize(shapeMU, [n, 3]).T
這里的shape = shapePC.dot(sp)，即$\sum _{i=1}^m{a}_i{S}_i$：
shape_3d = shape

到這里b就可以按照公式計算出來了，得到的b格式為(2,68)
b = A.dot(mu_3d + shape_3d) + np.tile(t2d[:, np.newaxis], [1, n])
然后將b轉換為格式(68*2,1)
b = np.reshape(b.T, [-1, 1])

求取$\beta$

完成A、pc和b的定義和計算之后$X_{projection}$的公式就可以寫成：

帶入pc的式子可以寫成：

將$X_{projection}$的公式帶入能量方程：

得到

對$\beta$進行求導，得到導數為零時$\beta$的取值。
L2范數求導可以使用公式：

得到

化簡得到

接著就是求取$\beta$的代碼：

equation_left = np.dot(pc.T, pc) + lamb * np.diagflat(1/sigma**2) x = np.reshape(x.T, [-1, 1]) equation_right = np.dot(pc.T, x - b) exp_para = np.dot(np.linalg.inv(equation_left), equation_right)

(4).將(2)和(3)中求出的$\alpha$代入能量方程，解得$\alpha$

同理，求取$\alpha$的代碼如下：

expression = expPC.dot(ep) expression = np.reshape(expression, [int(len(expression)/3), 3]).T sp = estimate_shape(x, shapeMU, shapePC, model['shapeEV'][:n_sp,:], expression, s, R, t[:2], lamb = 0.004)

算法過程與求取$\beta$相同，但是帶入的$\beta$是經過上面計算后的新值。

(5)更新$\alpha ,\beta$的值，重復(2)-(4)進行迭代更新。

到這里，循環迭代部分的代碼告一段落，經過多次迭代計算（程序中給的迭代次數為三次），獲得了所需要的sp, ep, s, R, t。

回到例程

回到bfm.fit從
fitted_sp, fitted_ep, s, R, t = fit.fit_points(x, X_ind, self.model, n_sp = self.n_shape_para, n_ep = self.n_exp_para, max_iter = max_iter)
繼續向下執行:

def fit(self, x, X_ind, max_iter = 4, isShow = False):''' fit 3dmm & pose parametersArgs:x: (n, 2) image pointsX_ind: (n,) corresponding Model vertex indicesmax_iter: iterationisShow: whether to reserve middle results for showReturns:fitted_sp: (n_sp, 1). shape parametersfitted_ep: (n_ep, 1). exp parameterss, angles, t'''if isShow:fitted_sp, fitted_ep, s, R, t = fit.fit_points_for_show(x, X_ind, self.model, n_sp = self.n_shape_para, n_ep = self.n_exp_para, max_iter = max_iter)angles = np.zeros((R.shape[0], 3))for i in range(R.shape[0]):angles[i] = mesh.transform.matrix2angle(R[i])else:fitted_sp, fitted_ep, s, R, t = fit.fit_points(x, X_ind, self.model, n_sp = self.n_shape_para, n_ep = self.n_exp_para, max_iter = max_iter)angles = mesh.transform.matrix2angle(R)return fitted_sp, fitted_ep, s, angles, t

將旋轉矩陣轉換為XYZ角度angles = mesh.transform.matrix2angle(R)
之后返回fitted_sp, fitted_ep, s, angles, t。

回到3DMM例程
fitted_sp, fitted_ep, fitted_s, fitted_angles, fitted_t = bfm.fit(x, X_ind, max_iter = 3)
部分執行完畢，繼續向下執行：

x = projected_vertices[bfm.kpt_ind, :2] # 2d keypoint, which can be detected from image X_ind = bfm.kpt_ind # index of keypoints in 3DMM. fixed.# fit fitted_sp, fitted_ep, fitted_s, fitted_angles, fitted_t = bfm.fit(x, X_ind, max_iter = 3)# verify fitted parameters fitted_vertices = bfm.generate_vertices(fitted_sp, fitted_ep) transformed_vertices = bfm.transform(fitted_vertices, fitted_s, fitted_angles, fitted_t)image_vertices = mesh.transform.to_image(transformed_vertices, h, w) fitted_image = mesh.render.render_colors(image_vertices, bfm.triangles, colors, h, w)

接下來是
fitted_vertices = bfm.generate_vertices(fitted_sp, fitted_ep)
根據計算出的$\alpha ,\beta$代入

算出$S_{newModel}$,對應的源碼如下：

def generate_vertices(self, shape_para, exp_para):'''Args:shape_para: (n_shape_para, 1)exp_para: (n_exp_para, 1) Returns:vertices: (nver, 3)'''vertices = self.model['shapeMU'] + \self.model['shapePC'].dot(shape_para) + \self.model['expPC'].dot(exp_para)vertices = np.reshape(vertices, [int(3), int(len(vertices)/3)], 'F').Treturn vertices

算出$S_{newModel}$后再對三維模型進行相似變換：

transformed_vertices = bfm.transform(fitted_vertices, fitted_s, fitted_angles, fitted_t)

即$S_{tranformed}=s?R?S_{newModel}+t_{3d}$

def transform(self, vertices, s, angles, t3d):R = mesh.transform.angle2matrix(angles)return mesh.transform.similarity_transform(vertices, s, R, t3d)

之后的代碼：

image_vertices = mesh.transform.to_image(transformed_vertices, h, w) fitted_image = mesh.render.render_colors(image_vertices, bfm.triangles, colors, h, w)

分別是將三維模型轉為二維圖像的格式并用自帶的顏色信息對模型上色，這里不做過多解讀了。

結果展示

生成的新的人臉圖片被保存到results/3dmm目錄下：

# ------------- print & show print('pose, groudtruth: \n', s, angles[0], angles[1], angles[2], t[0], t[1]) print('pose, fitted: \n', fitted_s, fitted_angles[0], fitted_angles[1], fitted_angles[2], fitted_t[0], fitted_t[1])save_folder = 'results/3dmm' if not os.path.exists(save_folder):os.mkdir(save_folder)io.imsave('{}/generated.jpg'.format(save_folder), image) io.imsave('{}/fitted.jpg'.format(save_folder), fitted_image)

還可以通過生成一個gif展示特征點擬合的過程：

# fitfitted_sp, fitted_ep, fitted_s, fitted_angles, fitted_t = bfm.fit(x, X_ind, max_iter = 3, isShow = True)# verify fitted parameters for i in range(fitted_sp.shape[0]):fitted_vertices = bfm.generate_vertices(fitted_sp[i], fitted_ep[i])transformed_vertices = bfm.transform(fitted_vertices, fitted_s[i], fitted_angles[i], fitted_t[i])image_vertices = mesh.transform.to_image(transformed_vertices, h, w)fitted_image = mesh.render.render_colors(image_vertices, bfm.triangles, colors, h, w)io.imsave('{}/show_{:0>2d}.jpg'.format(save_folder, i), fitted_image)options = '-delay 20 -loop 0 -layers optimize' # gif. need ImageMagick. subprocess.call('convert {} {}/show_*.jpg {}'.format(options, save_folder, save_folder + '/3dmm.gif'), shell=True) subprocess.call('rm {}/show_*.jpg'.format(save_folder), shell=True)

下面是結果展示：

• generated.jpg
  
• fitted.jpg

• 3dmm.gif

當然，每次執行程序生成的新的隨機模型的樣子也會不同。

結語

歷時將近一周，終于將3DMM例程部分從頭到尾分析了一遍，多了許多收獲同時也有不少疑惑，我以后也會不斷更新去完善這幾篇文章，期待大家的關注。

總結

以上是生活随笔為你收集整理的Face3D学习笔记（6）3DMM示例源码解析【下】从二维图片的特征点重建三维模型的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。