前言

關(guān)于使用3DMM進(jìn)行人臉重建的方法已經(jīng)有不少，相關(guān)開(kāi)源的代碼也有不少，但少有用python來(lái)寫(xiě)的，突然在github上發(fā)現(xiàn)這個(gè)face3d，感覺(jué)不錯(cuò)分享一下
github地址：https://github.com/YadiraF/face3d
在使用3DMM方法人臉建模的時(shí)候，最大的問(wèn)題便是shape系數(shù)以及exp系數(shù)的確定，但是在論文中往往是一兩句話便略過(guò)，使得像我這樣的初學(xué)者總是一頭霧水。正好遇到這個(gè)face3d，對(duì)里面的fitting部分進(jìn)行分析，以便以后使用。

估計(jì)目標(biāo)

這里先對(duì)3DMM方法進(jìn)行一個(gè)簡(jiǎn)介，具體可以參考99年的這篇論文https://blog.csdn.net/likewind1993/article/details/79177566，論文里提出了人臉的一種線性表示方法，從而一個(gè)新的人臉形狀可以由以下的方法得到（原線性表示中還有一個(gè)紋理部分，但是擬合效果往往不夠好，一般直接從照片中提取紋理進(jìn)行貼合，因此這里只考慮重建人臉形狀的部分）：

$$S_{newModel}=\bar{S}+\sum _{i=1}^{m?1}{\alpha }_i{s}_i$$
其中$\bar{S}$表示平均臉部模型，$s_i$表示shape對(duì)應(yīng)的PCA部分，${\alpha }_i$表示相應(yīng)的系數(shù)，這樣一張新人臉形狀就可依照上式得到。在14年的時(shí)候，FacewareHouse這篇論文提出并公開(kāi)了一個(gè)人臉表情數(shù)據(jù)庫(kù)，使得3DMM有了更強(qiáng)的表現(xiàn)力。從而人臉模型的線性表示可以擴(kuò)充為：

$$S_{newModel}=\bar{S}+\sum _{i=1}^{m?1}{\alpha }_i{s}_i+\sum _{i=1}^{n?1}{\beta }_i{e}_i$$
加入了e（expression，表情）。

于是人臉重建問(wèn)題轉(zhuǎn)為了求$\alpha ,\beta$系數(shù)的問(wèn)題。

得到一張單張正臉照片，可以從里面得到人臉的68個(gè)特征點(diǎn)坐標(biāo)（$X$)，在BFM模型中有對(duì)應(yīng)的68個(gè)特征點(diǎn)（$X_{3d}$），根據(jù)這些信息便可以求出$\alpha ,\beta$系數(shù)，將平均臉模型與照片中的臉部進(jìn)行擬合。

具體求解過(guò)程如下：

$$X_{projection}=s?P_{orth}?R?(\bar{S}+\sum _{i=1}^{m?1}{\alpha }_i{s}_i+\sum _{i=1}^{n?1}{\beta }_i{e}_i)+t_{2d}$$
這里$X_{projection}$是三維模型投影到二維平面的點(diǎn)，$P_{orth}=[[1,0,0],[0,1,0]]$為正交投影矩陣，R（3，3）為旋轉(zhuǎn)矩陣，$t_{2d}$為位移矩陣。

于是乎，再一次三維求解問(wèn)題又可以轉(zhuǎn)化為求解滿足以下的能量方程的系數(shù)（$s,R,t_{2d}$, 以及$\alpha$和$\beta$）

$$argmin||X_{projection}?X|{|}^2+\lambda \sum _{i=1}（\frac{{\gamma }_i}{{\sigma }_i}{）}^2$$
這里加了正則化部分，其中$\gamma$是PCA系數(shù)（包括形狀系數(shù)$\alpha$以及表情系數(shù)$\beta$），$\sigma$表示對(duì)應(yīng)的主成分偏差。

即，由上式求解使得三維模型中的68特征點(diǎn)投影到二維平面上的值與二維平面原68個(gè)特征點(diǎn)距離相差最小的系數(shù)。

一般論文到這里便戛然而止，對(duì)于如何求解并沒(méi)有給出詳細(xì)的過(guò)程，頂多只給出一個(gè)“使用最小二乘法”進(jìn)行求解。

這里我們便往下進(jìn)行深挖一下如何對(duì)這個(gè)問(wèn)題進(jìn)行求解：

盤(pán)點(diǎn)下，我們需要求得的參數(shù)主要有$s,R,t_{2d},\alpha ,\beta$，

這里可以把參數(shù)分為三個(gè)部分，$s,R,t_{2d}$, 以及$\alpha$和$\beta$。

求解方法如下：

將$\alpha$以及$\beta$初始化為0

求出$s,R,t_{2d}$

將上一步求出的$s,R,t_{2d}$代入，求出$\alpha$

將之前求出的$s,R,t_{2d}$以及$\alpha$代入，求出$\beta$

利用求得的$\alpha$以及$\beta$，重復(fù)2-4步驟進(jìn)行迭代。

這里的求解方法又可以分為兩類(lèi)，一類(lèi)是確定$s,R,t_{2d}$（根據(jù)對(duì)應(yīng)特征點(diǎn)可以得到一個(gè)仿射矩陣，對(duì)矩陣進(jìn)行分解可以得到$s,R,t_{2d}$）,另一類(lèi)是確定$\alpha$和$\beta$（這兩種系數(shù)雖然不同，但在式子中的位置、屬性差不多，都可以用最小二乘法求解）。

估計(jì)$s,R,t_{2d}$

人臉模型的三維點(diǎn)以及對(duì)應(yīng)照片中的二維點(diǎn)存在映射關(guān)系，這個(gè)可以由一個(gè)3x4的仿射矩陣進(jìn)行表示。
即
$$x_{2d}=P?X_{3d}$$
$P$即是我們需要求得仿射矩陣，作用在三維坐標(biāo)點(diǎn)上可以得到對(duì)應(yīng)二維點(diǎn)的坐標(biāo)。
這里使用黃金標(biāo)準(zhǔn)算法（Gold Standard algorithm來(lái)求解該仿射矩陣。

算法如下：
目標(biāo)：
給定n>=4組3維（$X_i$）到2維($x_i$)的圖像點(diǎn)對(duì)應(yīng)，確定仿射攝像機(jī)投影矩陣的最大似然估計(jì)。

• 歸一化，對(duì)于二維點(diǎn)($x_i$)，計(jì)算一個(gè)相似變換T，使得$\bar{x}=T{x}_i$，同樣的對(duì)于三維點(diǎn)，計(jì)算$\bar{X}=U{X}_i$
• 對(duì)于每組對(duì)應(yīng)點(diǎn)$x_i$~$X_i$，都有上圖方程的對(duì)應(yīng)關(guān)系存在，形如$Ax=b$
• 求出A的偽逆
• 去掉歸一化，得到仿射矩陣。

（這里講的過(guò)于抽象，具體歸一化做了哪些操作以及相應(yīng)的推導(dǎo)可以參考“Multiple View Geometry in Computer Vision”這本書(shū)，以及相關(guān)的代碼實(shí)現(xiàn)，這里就不細(xì)談了）

在得到仿射矩陣$P$之后，需要將$P$進(jìn)行分解，得到縮放系數(shù)s，旋轉(zhuǎn)矩陣R，以及位移矩陣t，下面是face3d里具體的實(shí)現(xiàn)代碼：

def P2sRt(P):''' decompositing camera matrix PArgs: P: (3, 4). Affine Camera Matrix.Returns:s: scale factor.R: (3, 3). rotation matrix.t: (3,). translation. '''t = P[:, 3]R1 = P[0:1, :3]R2 = P[1:2, :3]s = (np.linalg.norm(R1) + np.linalg.norm(R2))/2.0r1 = R1/np.linalg.norm(R1)r2 = R2/np.linalg.norm(R2)r3 = np.cross(r1, r2)R = np.concatenate((r1, r2, r3), 0)return s, R, t

估計(jì)$\alpha$和$\beta$

下面的部分來(lái)源于face3d代碼，為了將這個(gè)問(wèn)題說(shuō)的更明白點(diǎn)，這里搬運(yùn)一下（順便對(duì)原式進(jìn)行了簡(jiǎn)化，去掉了累加和等）。
令$$A=s?P?R$$$$b=s?P?R?(\bar{S}+e\beta )+t_{2}d$$
可以得到，
$$X_{projection}=A?s\alpha +b$$
Note:以下部分較多的涉及到各個(gè)變量的維數(shù)變化，顯得有些啰嗦甚至多余，但正是這些部分往往對(duì)我們解這些方程的時(shí)候起了較大的阻礙。

以BFM模型為例，設(shè)模型頂點(diǎn)總數(shù)為n, 得到各個(gè)參數(shù)的維數(shù)為：
$s$（3n x 199）
$\alpha$ (199 x 1)
$A$ (2 x 3)
在進(jìn)行$A?s\alpha$的時(shí)候，需要后者進(jìn)行一下變換，否則不滿足矩陣相乘要求。
$$shape=reshape(s\alpha )$$
使得$(3n$ x $1$)===》($3$ x $n)$
$A$ * $shape$的維數(shù)為（2 x n）
為了求解方便，對(duì)$A$ * shape， $b$進(jìn)行reshape
使其維數(shù)變?yōu)?#xff08;2n x 1）
得到
$$x_{flatten}=A?shape+b_{flatten}$$
這里重新定義

$$p{c}_{2d}=A?reshape(s)$$
這步在具體實(shí)現(xiàn)的時(shí)候，由于$s$的維數(shù)為（3n, 199），進(jìn)行reshape將維數(shù)調(diào)整為（199n, 3），這樣可以進(jìn)行如下的運(yùn)算$reshape(s)?A.T$，得到$p{c}_{2d}$，其維數(shù)為（199n, 2）。

$$p{c}_{2d_{f}latten}=reshape(p{c}_{2}d)$$
這里再次對(duì)$p{c}_{2d}$進(jìn)行維數(shù)的調(diào)整，調(diào)整為（2n, 199）,可以與維數(shù)為(199, 1)的$\alpha$直接相乘。
得到
$$x_{flatten}=p{c}_{2d_{f}latten}\alpha +b_{flatten}$$

轉(zhuǎn)回到原有的能量方程：
注：這里后加入項(xiàng)是正則項(xiàng)，$\gamma$中有$\alpha ，\beta$系數(shù)，底下的$\sigma$為對(duì)應(yīng)系數(shù)的方差。正則項(xiàng)加入能使求解能量方程得到的$\alpha ,\beta$系數(shù)趨于穩(wěn)定，避免過(guò)擬合。

$$argmin||X_{projection}?X|{|}^2+\lambda \sum _{i=1}（\frac{{\gamma }_i}{{\sigma }_i}{）}^2$$
將這里已經(jīng)推導(dǎo)好的式子代入，
得到$$E=||p{c}_{2d_{f}latten}\alpha +b_{flatten}?X|{|}^2+\lambda \sum _{i=1}（\frac{{\gamma }_i}{{\sigma }_i}{）}^2$$
對(duì)$\alpha$進(jìn)行求導(dǎo)，得到導(dǎo)數(shù)為零時(shí)$\alpha$的取值即為我們要求的。
$d(E)/d(\alpha )=0$

（這里涉及到了L2范數(shù)的求導(dǎo)，類(lèi)比這里的結(jié)論，很容易的得到以下的結(jié)果，后面的$\alpha$是對(duì)$\gamma$進(jìn)行求導(dǎo)的結(jié)果）

即
$$2?p{c}_{2d_{f}latten}.T?(p{c}_{2d_{f}latten}\alpha +b_{flatten}?X)+2?\lambda ?\frac{\alpha }{\sigma .T?\sigma }=0$$
整理得(為了簡(jiǎn)化式子，這里去掉flatten等后綴)：
$$(pc.T?pc+\frac{\lambda }{\sigma .T?\sigma })\alpha =pc.T?(x?b)$$
解出這個(gè)方程即可得到$\alpha$的值。

face3d中的代碼如下：

def estimate_shape(x, shapeMU, shapePC, shapeEV, expression, s, R, t2d, lamb = 3000):'''Args:x: (2, n). image points (to be fitted)shapeMU: (3n, 1)shapePC: (3n, n_sp)shapeEV: (n_sp, 1)expression: (3, n)s: scaleR: (3, 3). rotation matrixt2d: (2,). 2d translationlambda: regulation coefficientReturns:shape_para: (n_sp, 1) shape parameters(coefficients)'''x = x.copy()assert(shapeMU.shape[0] == shapePC.shape[0])assert(shapeMU.shape[0] == x.shape[1]*3)dof = shapePC.shape[1]n = x.shape[1]sigma = shapeEVt2d = np.array(t2d)P = np.array([[1, 0, 0], [0, 1, 0]], dtype = np.float32)A = s*P.dot(R)# --- calc pcpc_3d = np.resize(shapePC.T, [dof, n, 3]) # 199 x n x 3pc_3d = np.reshape(pc_3d, [dof*n, 3]) # 199n x 3pc_2d = pc_3d.dot(A.T.copy()) # A.T 3 x 2 199 x n x 2pc = np.reshape(pc_2d, [dof, -1]).T # 2n x 199# --- calc b# shapeMUmu_3d = np.resize(shapeMU, [n, 3]).T # 3 x n# expressionexp_3d = expression# b = A.dot(mu_3d + exp_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)shape_para = np.dot(np.linalg.inv(equation_left), equation_right)return shape_para

得到$\alpha$的值后，將值代入繼續(xù)進(jìn)行$\beta$的求解。重復(fù)迭代，即可求出相應(yīng)的解。

后記

雖然3DMM方法自94年提出以來(lái)到現(xiàn)在已經(jīng)有二十多年了，但是對(duì)于剛?cè)腴T(mén)三維人臉重建的新手來(lái)說(shuō)，作為練手仍然是不錯(cuò)的一個(gè)小項(xiàng)目。

總結(jié)

以上是生活随笔為你收集整理的使用3DMM进行人脸重建中的配准方法的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

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