前言
关于使用3DMM进行人脸重建的方法已经有不少,相关开源的代码也有不少,但少有用python来写的,突然在github上发现这个face3d,感觉不错分享一下
github地址:https://github.com/YadiraF/face3d
在使用3DMM方法人脸建模的时候,最大的问题便是shape系数以及exp系数的确定,但是在论文中往往是一两句话便略过,使得像我这样的初学者总是一头雾水。正好遇到这个face3d,对里面的fitting部分进行分析,以便以后使用。
估计目标
这里先对3DMM方法进行一个简介,具体可以参考99年的这篇论文https://blog.csdn.net/likewind1993/article/details/79177566,论文里提出了人脸的一种线性表示方法,从而一个新的人脸形状可以由以下的方法得到(原线性表示中还有一个纹理部分,但是拟合效果往往不够好,一般直接从照片中提取纹理进行贴合,因此这里只考虑重建人脸形状的部分):
S
n
e
w
M
o
d
e
l
=
S
ˉ
+
∑
i
=
1
m
−
1
α
i
s
i
S_{newModel} = \bar S + \sum_{i=1}^{m-1} \alpha_{i}s_{i}
SnewModel=Sˉ+i=1∑m−1αisi
其中
S
ˉ
\bar S
Sˉ表示平均脸部模型,
s
i
s_i
si表示shape对应的PCA部分,
α
i
\alpha_i
αi表示相应的系数,这样一张新人脸形状就可依照上式得到。在14年的时候,FacewareHouse这篇论文提出并公开了一个人脸表情数据库,使得3DMM有了更强的表现力。从而人脸模型的线性表示可以扩充为:
S
n
e
w
M
o
d
e
l
=
S
ˉ
+
∑
i
=
1
m
−
1
α
i
s
i
+
∑
i
=
1
n
−
1
β
i
e
i
S_{newModel} = \bar S + \sum_{i=1}^{m-1} \alpha_{i}s_{i}+ \sum_{i=1}^{n-1} \beta_{i}e_{i}
SnewModel=Sˉ+i=1∑m−1αisi+i=1∑n−1βiei
加入了e(expression,表情)。
于是人脸重建问题转为了求 α , β \alpha , \beta α,β系数的问题。
得到一张单张正脸照片,可以从里面得到人脸的68个特征点坐标( X X X),在BFM模型中有对应的68个特征点( X 3 d X_{3d} X3d),根据这些信息便可以求出 α , β \alpha , \beta α,β系数,将平均脸模型与照片中的脸部进行拟合。
具体求解过程如下:
X
p
r
o
j
e
c
t
i
o
n
=
s
∗
P
o
r
t
h
∗
R
∗
(
S
ˉ
+
∑
i
=
1
m
−
1
α
i
s
i
+
∑
i
=
1
n
−
1
β
i
e
i
)
+
t
2
d
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}
Xprojection=s∗Porth∗R∗(Sˉ+i=1∑m−1αisi+i=1∑n−1βiei)+t2d
这里
X
p
r
o
j
e
c
t
i
o
n
X_{projection}
Xprojection是三维模型投影到二维平面的点,
P
o
r
t
h
=
[
[
1
,
0
,
0
]
,
[
0
,
1
,
0
]
]
P_{orth} = [[1, 0, 0], [0, 1, 0]]
Porth=[[1,0,0],[0,1,0]]为正交投影矩阵,R(3,3)为旋转矩阵,
t
2
d
t_{2d}
t2d为位移矩阵。
于是乎,再一次三维求解问题又可以转化为求解满足以下的能量方程的系数( s , R , t 2 d {s, R, t_{2d}} s,R,t2d, 以及 α \alpha α和 β \beta β)
a
r
g
m
i
n
∣
∣
X
p
r
o
j
e
c
t
i
o
n
−
X
∣
∣
2
+
λ
∑
i
=
1
(
γ
i
σ
i
)
2
arg min ||X_{projection} - X||^2 + \lambda \sum_{i=1}(\frac{\gamma_i}{\sigma_i})^2
argmin∣∣Xprojection−X∣∣2+λi=1∑(σiγi)2
这里加了正则化部分,其中
γ
\gamma
γ是PCA系数(包括形状系数
α
\alpha
α以及表情系数
β
\beta
β),
σ
\sigma
σ表示对应的主成分偏差。
即,由上式求解使得三维模型中的68特征点投影到二维平面上的值与二维平面原68个特征点距离相差最小的系数。
一般论文到这里便戛然而止,对于如何求解并没有给出详细的过程,顶多只给出一个“使用最小二乘法”进行求解。
这里我们便往下进行深挖一下如何对这个问题进行求解:
盘点下,我们需要求得的参数主要有 s , R , t 2 d , α , β {s, R, t_{2d}, \alpha, \beta} s,R,t2d,α,β,
这里可以把参数分为三个部分, s , R , t 2 d {s, R, t_{2d}} s,R,t2d, 以及 α \alpha α和 β \beta β。
求解方法如下:
- 将 α \alpha α以及 β \beta β初始化为0
- 求出 s , R , t 2 d {s, R, t_{2d}} s,R,t2d
- 将上一步求出的 s , R , t 2 d {s, R, t_{2d}} s,R,t2d代入,求出 α \alpha α
- 将之前求出的 s , R , t 2 d {s, R, t_{2d}} s,R,t2d以及 α \alpha α代入,求出 β \beta β
- 利用求得的 α \alpha α以及 β \beta β,重复2-4步骤进行迭代。
这里的求解方法又可以分为两类,一类是确定 s , R , t 2 d {s, R, t_{2d}} s,R,t2d(根据对应特征点可以得到一个仿射矩阵,对矩阵进行分解可以得到 s , R , t 2 d {s, R, t_{2d}} s,R,t2d),另一类是确定 α \alpha α和 β \beta β(这两种系数虽然不同,但在式子中的位置、属性差不多,都可以用最小二乘法求解)。
估计 s , R , t 2 d {s, R, t_{2d}} s,R,t2d
人脸模型的三维点以及对应照片中的二维点存在映射关系,这个可以由一个3x4的仿射矩阵进行表示。
即
x
2
d
=
P
∗
X
3
d
x_{2d} = P * X_{3d}
x2d=P∗X3d
P
P
P即是我们需要求得仿射矩阵,作用在三维坐标点上可以得到对应二维点的坐标。
这里使用黄金标准算法(Gold Standard algorithm来求解该仿射矩阵。
算法如下:
目标:
给定n>=4组3维(
X
i
X_i
Xi)到2维(
x
i
x_i
xi)的图像点对应,确定仿射摄像机投影矩阵的最大似然估计。
- 归一化,对于二维点( x i x_i xi),计算一个相似变换T,使得 x ˉ = T x i \bar x = T x_{i} xˉ=Txi,同样的对于三维点,计算 X ˉ = U X i \bar X = U X_{i} Xˉ=UXi
- 对于每组对应点 x i x_{i} xi~ X i X_{i} Xi,都有上图方程的对应关系存在,形如 A x = b A x = b Ax=b
- 求出A的伪逆
- 去掉归一化,得到仿射矩阵。
(这里讲的过于抽象,具体归一化做了哪些操作以及相应的推导可以参考“Multiple View Geometry in Computer Vision”这本书,以及相关的代码实现,这里就不细谈了)
在得到仿射矩阵 P P P之后,需要将 P P P进行分解,得到缩放系数s,旋转矩阵R,以及位移矩阵t,下面是face3d里具体的实现代码:
def P2sRt(P):
''' decompositing camera matrix P
Args:
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.0
r1 = 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
估计 α \alpha α和 β \beta β
下面的部分来源于face3d代码,为了将这个问题说的更明白点,这里搬运一下(顺便对原式进行了简化,去掉了累加和等)。
令
A
=
s
∗
P
∗
R
A = s * P * R
A=s∗P∗R
b
=
s
∗
P
∗
R
∗
(
S
ˉ
+
e
β
)
+
t
2
d
b = s * P * R * ( \bar S + e \beta) + t_2d
b=s∗P∗R∗(Sˉ+eβ)+t2d
可以得到,
X
p
r
o
j
e
c
t
i
o
n
=
A
∗
s
α
+
b
X_{projection} = A * s \alpha + b
Xprojection=A∗sα+b
Note:以下部分较多的涉及到各个变量的维数变化,显得有些啰嗦甚至多余,但正是这些部分往往对我们解这些方程的时候起了较大的阻碍。
以BFM模型为例,设模型顶点总数为n, 得到各个参数的维数为:
s
s
s(3n x 199)
α
\alpha
α (199 x 1)
A
A
A (2 x 3)
在进行
A
∗
s
α
A * s \alpha
A∗sα的时候,需要后者进行一下变换,否则不满足矩阵相乘要求。
s
h
a
p
e
=
r
e
s
h
a
p
e
(
s
α
)
shape = reshape(s \alpha)
shape=reshape(sα)
使得
(
3
n
(3n
(3n x
1
1
1)===》(
3
3
3 x
n
)
n)
n)
A
A
A *
s
h
a
p
e
shape
shape的维数为(2 x n)
为了求解方便,对
A
A
A * shape,
b
b
b进行reshape
使其维数变为(2n x 1)
得到
x
f
l
a
t
t
e
n
=
A
∗
s
h
a
p
e
+
b
f
l
a
t
t
e
n
x_{flatten} = A * shape + b_{flatten}
xflatten=A∗shape+bflatten
这里重新定义
p
c
2
d
=
A
∗
r
e
s
h
a
p
e
(
s
)
pc_{2d} = A * reshape(s)
pc2d=A∗reshape(s)
这步在具体实现的时候,由于
s
s
s的维数为(3n, 199),进行reshape将维数调整为(199n, 3),这样可以进行如下的运算
r
e
s
h
a
p
e
(
s
)
∗
A
.
T
reshape(s) * A.T
reshape(s)∗A.T,得到
p
c
2
d
pc_{2d}
pc2d,其维数为(199n, 2)。
p
c
2
d
f
l
a
t
t
e
n
=
r
e
s
h
a
p
e
(
p
c
2
d
)
pc_{2d_flatten} = reshape(pc_2d)
pc2dflatten=reshape(pc2d)
这里再次对
p
c
2
d
pc_{2d}
pc2d进行维数的调整,调整为(2n, 199),可以与维数为(199, 1)的
α
\alpha
α直接相乘。
得到
x
f
l
a
t
t
e
n
=
p
c
2
d
f
l
a
t
t
e
n
α
+
b
f
l
a
t
t
e
n
x_{flatten} = pc_{2d_flatten} \alpha + b_{flatten}
xflatten=pc2dflattenα+bflatten
转回到原有的能量方程:
注:这里后加入项是正则项,
γ
\gamma
γ中有
α
,
β
\alpha,\beta
α,β系数,底下的
σ
\sigma
σ为对应系数的方差。正则项加入能使求解能量方程得到的
α
,
β
\alpha, \beta
α,β系数趋于稳定,避免过拟合。
a
r
g
m
i
n
∣
∣
X
p
r
o
j
e
c
t
i
o
n
−
X
∣
∣
2
+
λ
∑
i
=
1
(
γ
i
σ
i
)
2
arg min ||X_{projection} - X||^2 + \lambda \sum_{i=1} (\frac{\gamma_i}{\sigma_i})^2
argmin∣∣Xprojection−X∣∣2+λi=1∑(σiγi)2
将这里已经推导好的式子代入,
得到
E
=
∣
∣
p
c
2
d
f
l
a
t
t
e
n
α
+
b
f
l
a
t
t
e
n
−
X
∣
∣
2
+
λ
∑
i
=
1
(
γ
i
σ
i
)
2
E = ||pc_{2d_flatten} \alpha + b_{flatten} - X||^2 + \lambda \sum_{i=1} (\frac{\gamma_i}{\sigma_i})^2
E=∣∣pc2dflattenα+bflatten−X∣∣2+λi=1∑(σiγi)2
对
α
\alpha
α进行求导,得到导数为零时
α
\alpha
α的取值即为我们要求的。
d
(
E
)
/
d
(
α
)
=
0
d(E)/d(\alpha)= 0
d(E)/d(α)=0
(这里涉及到了L2范数的求导,类比这里的结论,很容易的得到以下的结果,后面的 α \alpha α是对 γ \gamma γ进行求导的结果)
即
2
∗
p
c
2
d
f
l
a
t
t
e
n
.
T
∗
(
p
c
2
d
f
l
a
t
t
e
n
α
+
b
f
l
a
t
t
e
n
−
X
)
+
2
∗
λ
∗
α
σ
.
T
∗
σ
=
0
2 * pc_{2d_flatten}.T*(pc_{2d_flatten} \alpha + b_{flatten} - X)+2 * \lambda *\frac{\alpha}{\sigma.T* \sigma} = 0
2∗pc2dflatten.T∗(pc2dflattenα+bflatten−X)+2∗λ∗σ.T∗σα=0
整理得(为了简化式子,这里去掉flatten等后缀):
(
p
c
.
T
∗
p
c
+
λ
σ
.
T
∗
σ
)
α
=
p
c
.
T
∗
(
x
−
b
)
(pc.T* pc + \frac{\lambda }{\sigma.T* \sigma}) \alpha = pc.T * (x-b)
(pc.T∗pc+σ.T∗σλ)α=pc.T∗(x−b)
解出这个方程即可得到
α
\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: scale
R: (3, 3). rotation matrix
t2d: (2,). 2d translation
lambda: regulation coefficient
Returns:
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 = shapeEV
t2d = np.array(t2d)
P = np.array([[1, 0, 0], [0, 1, 0]], dtype = np.float32)
A = s*P.dot(R)
# --- calc pc
pc_3d = np.resize(shapePC.T, [dof, n, 3]) # 199 x n x 3
pc_3d = np.reshape(pc_3d, [dof*n, 3]) # 199n x 3
pc_2d = pc_3d.dot(A.T.copy()) # A.T 3 x 2 199 x n x 2
pc = np.reshape(pc_2d, [dof, -1]).T # 2n x 199
# --- calc b
# shapeMU
mu_3d = np.resize(shapeMU, [n, 3]).T # 3 x n
# expression
exp_3d = expression
#
b = A.dot(mu_3d + exp_3d) + np.tile(t2d[:, np.newaxis], [1, n]) # 2 x n
b = np.reshape(b.T, [-1, 1]) # 2n x 1
# --- solve
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)
shape_para = np.dot(np.linalg.inv(equation_left), equation_right)
return shape_para
得到 α \alpha α的值后,将值代入继续进行 β \beta β的求解。重复迭代,即可求出相应的解。
后记
虽然3DMM方法自94年提出以来到现在已经有二十多年了,但是对于刚入门三维人脸重建的新手来说,作为练手仍然是不错的一个小项目。