寫在前面

• 為了保證整個(gè)示例項(xiàng)目更加直觀，方便理解，在展示一些函數(shù)的源碼時(shí)會(huì)使用numpy版本進(jìn)行展示，而在示例程序中并未使用numpy版本的庫(kù)，在Cython版本與numpy版本出現(xiàn)差異的原碼前會(huì)有標(biāo)注，希望讀者留意。
• 3DMM實(shí)例程序的jupyter版本后續(xù)會(huì)更新，完全免費(fèi)，歡迎大家下載

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

代碼解析

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

在上一篇文章我們通過(guò)黃金標(biāo)準(zhǔn)算法求解出了仿射矩陣$P_A$并將它分解的到了$s,R,t_{2d}$，這部分將繼續(xù)按照求解步驟進(jìn)行源碼分析:

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

對(duì)于公式

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

shape的格式為（159645，1），再通過(guò)
shape = np.reshape(shape, [int(len(shape)/3), 3]).T
將shape的XYZ坐標(biāo)分開(kāi)，轉(zhuǎn)為（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

數(shù)據(jù)處理

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)

首先是確認(rèn)輸入的格式正確：
assert(shapeMU.shape[0] == expPC.shape[0])
assert(shapeMU.shape[0] == x.shape[1]*3)
然后此時(shí)輸入的表情主成分expPC格式為(159645,29)
令dof=29
dof = expPC.shape[1]
令n=68
n = x.shape[1]
另sigma=expEV 即表情主成分方差$\sigma$
sigma = expEV
$t_{2d}$轉(zhuǎn)化為array數(shù)組
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

已知公式

定義和計(jì)算A、pc和b：

• 定義A:

A = s*P.dot(R)

• 計(jì)算pc：

這里的pc計(jì)算相當(dāng)于下式

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

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

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

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

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

• 定義b

b的公式如下：

計(jì)算時(shí)由于矩陣的格式問(wèn)題要先進(jìn)行一些變換：
這里的shapeMU也是只包含了68個(gè)特征點(diǎn)的
將格式為(68*3,1)的shapeMU轉(zhuǎn)換為格式（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就可以按照公式計(jì)算出來(lái)了，得到的b格式為(2,68)
b = A.dot(mu_3d + shape_3d) + np.tile(t2d[:, np.newaxis], [1, n])
然后將b轉(zhuǎn)換為格式(68*2,1)
b = np.reshape(b.T, [-1, 1])

求取$\beta$

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

帶入pc的式子可以寫成：

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

得到

對(duì)$\beta$進(jìn)行求導(dǎo)，得到導(dǎo)數(shù)為零時(shí)$\beta$的取值。
L2范數(shù)求導(dǎo)可以使用公式：

得到

化簡(jiǎn)得到

接著就是求取$\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)

算法過(guò)程與求取$\beta$相同，但是帶入的$\beta$是經(jīng)過(guò)上面計(jì)算后的新值。

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

到這里，循環(huán)迭代部分的代碼告一段落，經(jīng)過(guò)多次迭代計(jì)算（程序中給的迭代次數(shù)為三次），獲得了所需要的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)
繼續(xù)向下執(zhí)行:

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

將旋轉(zhuǎn)矩陣轉(zhuǎn)換為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)
部分執(zhí)行完畢，繼續(xù)向下執(zhí)行：

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)

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

算出$S_{newModel}$,對(duì)應(yīng)的源碼如下：

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}$后再對(duì)三維模型進(jìn)行相似變換：

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)

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

結(jié)果展示

生成的新的人臉圖片被保存到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)

還可以通過(guò)生成一個(gè)gif展示特征點(diǎn)擬合的過(guò)程：

# 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)

下面是結(jié)果展示：

• generated.jpg
  
• fitted.jpg

• 3dmm.gif

當(dāng)然，每次執(zhí)行程序生成的新的隨機(jī)模型的樣子也會(huì)不同。

結(jié)語(yǔ)

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

總結(jié)

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

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。