人脸对齐 matlab,人脸对齐中的相似性变换

人脸对齐是大多数人脸分析算法中的一个关键模块，在人脸识别、表情识别、人脸建模等领域有着广泛的应用。人脸对齐获取图像中人脸的几何结构，基于平移、缩放和旋转得到对齐后的标准人脸。

在欧式几何中，如果两个物体具有相同的形状，或者其中一个物体的形状与另一个物体的镜像相同，那么这两个物体是相似的。更准确地说，可以通过均匀缩放(放大或缩小)并叠加必要的平移、旋转和反射来获得另一个。这意味着任意物体都可以重新缩放、重新定位和反射，以便与另一物体精确重合。如果两个物体相似，则每个物体都与另一个物体的特定均匀缩放结果一致。

目前大多数人脸对齐方案在定位关键点后采用 Umeyama 对齐算法进行处理。例如 insightface 中的 similarTransform，dlib 中的 get_face_chip_details。相似性变换由平移变换、旋转变换以及尺度变换组合而成。下面以 skimage.transform.SimilarityTransform 为例进行介绍。

SimilarityTransform

2D 相似性变换。具有以下形式：

X = a0 * x - b0 * y + a1 =

= s * x * cos(rotation) - s * y * sin(rotation) + a1

Y = b0 * x + a0 * y + b1 =

= s * x * sin(rotation) + s * y * cos(rotation) + b1

其中s是比例因子，齐次变换矩阵为:

[ a 0 b 0 a 1 b 0 a 0 b 1 0 0 1 ] \begin{bmatrix} a_0 & b_0 & a_1 \\ b_0 & a_0 & b_1 \\ 0 & 0 & 1 \end{bmatrix}⎣⎡​a0​b0​0​b0​a0​0​a1​b1​1​⎦⎤​

除了旋转和平移参数外，相似性变换还使用单个缩放因子扩展了欧几里德变换。

参数：

matrix：(3，3)数组，可选。齐次变换矩阵。

scale：float，可选。比例因子。

rotation：float，可选。逆时针旋转角度(以弧度表示)。

translation：(tx, ty)作为数组、列表或元组，可选。x，y 方向平移参数。

属性：

params：(3，3)数组，齐次变换矩阵。

[ s cos ⁡ ( θ ) sin ⁡ ( θ ) t x sin ⁡ ( θ ) s cos ⁡ ( θ ) t y 0 0 1 ] \begin{bmatrix} s \cos(\theta) & \sin(\theta) & t_x\\ \sin(\theta) & s\cos(\theta) & t_y\\ 0 & 0 & 1 \end{bmatrix}⎣⎡​scos(θ)sin(θ)0​sin(θ)scos(θ)0​tx​ty​1​⎦⎤​

def __init__(self, matrix=None, scale=None, rotation=None,

translation=None):

params = any(param is not None

for param in (scale, rotation, translation))

if params and matrix is not None:

raise ValueError("You cannot specify the transformation matrix and"

" the implicit parameters at the same time.")

elif matrix is not None:

if matrix.shape != (3, 3):

raise ValueError("Invalid shape of transformation matrix.")

self.params = matrix

elif params:

if scale is None:

scale = 1

if rotation is None:

rotation = 0

if translation is None:

translation = (0, 0)

self.params = np.array([

[math.cos(rotation), - math.sin(rotation), 0],

[math.sin(rotation), math.cos(rotation), 0],

[ 0, 0, 1]

])

self.params[0:2, 0:2] *= scale

self.params[0:2, 2] = translation

else:

# default to an identity transform

self.params = np.eye(3)

estimate

从一组对应点估计变换。可以使用总体最小二乘法确定过定、确定和欠定的参数。源坐标和目标坐标的数量必须匹配。

参数：

src：(N，2)数组,源坐标。

dst：(N，2)数组，目标坐标。

返回值：

success：布尔型，如果模型估计成功，则返回True。

self.params = _umeyama(src, dst, True)

return True

scale

@property

def scale(self):

if abs(math.cos(self.rotation)) < np.spacing(1):

# sin(self.rotation) == 1

scale = self.params[1, 0]

else:

scale = self.params[0, 0] / math.cos(self.rotation)

return scale

_umeyama

估算是否具有缩放比例的N-D相似度变换。

估计有或无标度的N-D相似变换。

参数：

src：(M，N)数组，源坐标。

dst：(M，N)数组，目标坐标。

estimate_scale：布尔，是否估计比例因子。

返回值：

T：(N + 1，N + 1)，齐次相似性变化矩阵。仅当问题条件不完善时，矩阵才会包含NaN值。

参考文献

[1] “Least-squares estimation of transformation parameters between two point patterns”, Shinji Umeyama, PAMI 1991, :DOI:10.1109/34.88573

μ x = 1 n ∑ i = 1 n x i (34) μ y = 1 n ∑ i = 1 n y i (35) σ x 2 = 1 n ∑ i = 1 n ∥ x i − μ x ∥ 2 (36) σ y 2 = 1 n ∑ i = 1 n ∥ y i − μ y ∥ 2 (37) Σ x y = 1 n ∑ i = 1 n ( y i − μ y ) ( x i − μ x ) ⊤ (38) \begin{aligned} \mu_x &= \frac{1}{n}\sum_{i=1}^n x_i \qquad\qquad\qquad \text{(34)}\\ \mu_y &= \frac{1}{n}\sum_{i=1}^n y_i \qquad\qquad\qquad \text{(35)}\\ \sigma_x^2 &= \frac{1}{n}\sum_{i=1}^n \| x_i-\mu_x\|^2 \qquad\qquad \text{(36)}\\ \sigma_y^2 &= \frac{1}{n}\sum_{i=1}^n \| y_i-\mu_y\|^2 \qquad\qquad \text{(37)}\\ \Sigma_{xy} &= \frac{1}{n}\sum_{i=1}^n (y_i-\mu_y)(x_i-\mu_x)^\top \quad \text{(38)} \end{aligned}μx​μy​σx2​σy2​Σxy​​=n1​i=1∑n​xi​(34)=n1​i=1∑n​yi​(35)=n1​i=1∑n​∥xi​−μx​∥2(36)=n1​i=1∑n​∥yi​−μy​∥2(37)=n1​i=1∑n​(yi​−μy​)(xi​−μx​)⊤(38)​

Σ x y \Sigma_{xy}Σxy​ 是 X XX 和 Y YY 的协方差矩阵，μ x \mu_xμx​ 和 μ y \mu_yμy​ 分别是 X XX 和 Y YY 的均值向量，σ x 2 \sigma_x^2σx2​ 和 σ y 2 \sigma_y^2σy2​ 是 X XX 和 Y YY 的方差。

num是点的数量，dim为点的坐标维度。

num = src.shape[0]

dim = src.shape[1]

# Compute mean of src and dst.

src_mean = src.mean(axis=0)

dst_mean = dst.mean(axis=0)

# Subtract mean from src and dst.

src_demean = src - src_mean

dst_demean = dst - dst_mean

# Eq. (38).

A = dst_demean.T @ src_demean / num

对协方差矩阵 Σ x y \Sigma_{xy}Σxy​ 进行奇异值分解 U D V ⊤ UDV^\topUDV⊤，其中 D = d i a g ( d i ) , d 1 ≥ d 2 ≥ ⋯ ≥ d m ≥ 0 D=\mathrm{diag}(d_i), d_1 \ge d_2 \ge\cdots \ge d_m \ge 0D=diag(di​),d1​≥d2​≥⋯≥dm​≥0。

S = { I if  d e t ( Σ x y ) ≥ 0 d i a g ( 1 , 1 , … , 1 , − 1 ) if  d e t ( Σ x y ) < 0 (39) S = \begin{cases} I &\text{if } \mathrm{det}(\Sigma_{xy}) \ge 0 \\ \mathrm{diag}(1,1,\dots,1,-1) &\text{if } \mathrm{det}(\Sigma_{xy}) < 0 \end{cases}\quad \text{(39)}S={Idiag(1,1,…,1,−1)​ifdet(Σxy​)≥0ifdet(Σxy​)<0​(39)

numpy.linalg.det 计算数组的行列式(determinant)。

numpy.linalg.svd 对数组进行奇异值分解。

S对应分解的奇异值 d i a g ( d i ) \mathrm{diag}(d_i)diag(di​)。

d为行向量，对应公式中的 S SS。

# Eq. (39).

d = np.ones((dim,), dtype=np.double)

if np.linalg.det(A) < 0:

d[dim - 1] = -1

T = np.eye(dim + 1, dtype=np.double)

U, S, V = np.linalg.svd(A)

numpy.linalg.matrix_rank 使用 SVD 方法返回数组的矩阵秩。

当 rank ( Σ x y ) = m \text{rank}{(\Sigma_{xy})}= mrank(Σxy​)=m 时，

S = { I if  d e t ( U ) d e t ( V ) = 1 d i a g ( 1 , 1 , … , 1 , − 1 ) if  d e t ( U ) d e t ( V ) = − 1 (43) S = \begin{cases} I &\text{if } \mathrm{det}(U)\mathrm{det}(V)=1 \\ \mathrm{diag}(1,1,\dots,1,-1) &\text{if } \mathrm{det}(U)\mathrm{det}(V)=-1 \end{cases} \quad \text{(43)}S={Idiag(1,1,…,1,−1)​ifdet(U)det(V)=1ifdet(U)det(V)=−1​(43)

当 rank ( Σ x y ) > m \text{rank}{(\Sigma_{xy})}> mrank(Σxy​)>m 时，最变换参数为：

R = U S V ⊤ (40) t = μ y − c R μ x (41) c = 1 σ x 2 t r ( D S ) (42) \begin{aligned} R &= USV^\top \qquad \text{(40)}\\ t &= \mu_y - cR\mu_x \quad \text{(41)}\\ c &= \frac{1}{\sigma_x^2}\mathrm{tr}(DS) \quad \text{(42)} \end{aligned}Rtc​=USV⊤(40)=μy​−cRμx​(41)=σx2​1​tr(DS)(42)​

s临时保存。

# Eq. (40) and (43).

rank = np.linalg.matrix_rank(A)

if rank == 0:

return np.nan * T

elif rank == dim - 1:

if np.linalg.det(U) * np.linalg.det(V) > 0:

T[:dim, :dim] = U @ V

else:

s = d[dim - 1]

d[dim - 1] = -1

T[:dim, :dim] = U @ np.diag(d) @ V

d[dim - 1] = s

else:

T[:dim, :dim] = U @ np.diag(d) @ V

首先计算尺度参数 c cc (公式42)；

然后计算平移参数 t tt(公式41)；

最后旋转参数 R RR 和尺度参数 c cc 融合到一起。

T的形式为：

[ s cos ⁡ ( θ ) sin ⁡ ( θ ) t x sin ⁡ ( θ ) s cos ⁡ ( θ ) t y 0 0 1 ] \begin{bmatrix} s \cos(\theta) & \sin(\theta) & t_x\\ \sin(\theta) & s\cos(\theta) & t_y\\ 0 & 0 & 1 \end{bmatrix}⎣⎡​scos(θ)sin(θ)0​sin(θ)scos(θ)0​tx​ty​1​⎦⎤​

if estimate_scale:

# Eq. (41) and (42).

scale = 1.0 / src_demean.var(axis=0).sum() * (S @ d)

else:

scale = 1.0

T[:dim, dim] = dst_mean - scale * (T[:dim, :dim] @ src_mean.T)

T[:dim, :dim] *= scale

return T

参考资料：

What does the “at” (@) symbol do in Python?

face alignment[Ordinary Procrustes Analysis]

Similarity (geometry)

How does dlib face aligment works? #1382

传统算法和深度学习的结合和实践，解读与优化 deepfake

Umeyama算法

similarity transform matrix in c++ is different from python #481

C++ using

4.3 Planar Graphs

Find All Cycles (Faces) In a Graph

Chinese Whispers

Chinese Whispers - an Efficient Graph Clustering Algorithm and its Application to Natural Language Processing Problems

Applying Affine transform on an image using dlib

dlib人脸关键点代码解析

Get Face Landmarks

dlib人脸对齐源码详解

dlib 人脸对齐 基本原理

对mtcnn的人脸对齐的理解

Review of similarity transformation and Singular Value Decomposition

图像的等距变换，相似变换，仿射变换，射影变换及其matlab实现

skimage库的transform.SimilarityTransform()用法

How to compute the similarity transformation matrix

Average Face : OpenCV ( C++ / Python ) Tutorial

SimilarityTransform

Aligning Face Images

Face Alignment

face alignment algorithm on images

dlib数据结构matrix

Dlib源码解析之一 matrix和array2d和image_view

Dlib Element Specific Operations

maketform

dlib.net/matrix_ex.cpp

deepfakes/faceswap/lib/umeyama.py

ethz-asl/maplab/test/end-to-end-common/python/end_to_end_common/umeyama.py

CarloNicolini/ralign

Singular Value Decomposition (SVD) tutorial