CN103440510A - 一种面部图像中特征点的定位方法 - Google Patents

Abstract

本发明公开了一种面部图像中特征点的定位方法，属于信号处理技术领域。本发明包括模型训练阶段和图像搜索阶段；在模型训练阶段，选取足够数量的正常人脸图像和不同光照下的人脸图像作为训练集合，手工标定好每幅训练图像准确的特征点位置，并建立活动形状模型和局部纹理模型，该模型训练阶段仅需执行一次；在图像搜索阶段，首先确定目标图像准确的瞳孔位置，以初始化活动形状模型的初始位置，然后通过迭代使模型拟合到准确的目标特征点位置。

Description

一种面部图像中特征点的定位方法

技术领域

本发明涉及一种面部图像中特征点的定位方法，属于信号处理技术领域。

背景技术

人的面部特征点定位是人脸识别、表情识别、性别识别、人脸美化等图像处理应用中的关键技术，它在视频监控、公共安全、数字娱乐等领域具有广泛应用。局部面部特征点定位方法是一种定位局部器官的方法，该方法计算量较小，但只能给出很少的特征点，定位精度较低，且不太稳健。基于活动形状模型(Active Shape Models,ASM)的特征点定位方法是一种全局定位方法，它可同时定位面部多个特征点，且定位速度较快，因而得到广泛应用，但该方法受光照影响较大，且对初始位置敏感。

发明内容

针对以上技术问题，本发明提出了一种面部图像中特征点的定位方法。

本发明采取的技术方案如下：

本发明包括模型训练阶段和图像搜索阶段；在模型训练阶段，选取足够数量的正常人脸图像和不同光照下的人脸图像作为训练集合，手工标定好每幅训练图像准确的特征点位置，并建立活动形状模型和局部纹理模型，该模型训练阶段仅需执行一次；在图像搜索阶段，首先确定目标图像准确的瞳孔位置，以初始化活动形状模型的初始位置，然后通过迭代使模型拟合到准确的目标特征点位置。

本发明技术方案带来的有益效果：

模型训练阶段，为每个特征点建立局部纹理模型时采用边长为25的正方形；图像搜索阶段，搜索邻域的长度为25个像素点，计算搜索邻域中每个点的LBP直方图时同样采用边长为25的正方形，权重wt、wθ、ws和wb的值都为1。用YefeiChen和Jianbo Su在论文“Fast eye localization based on a new Haar-like feature”（The10th World Congress on Intelligent Control and Automation,Beijing,China.2012,4825-4830）中所述的方法定位瞳孔位置，结果如图12所示。图13为用瞳孔位置得到的活动形状模型的初始位置的例子。在得到活动形状模型的初始位置后，用局部纹理模型与直方图相似性判别公式，对活动形状模型的特征点进行调整，使其收敛至准确的特征点位置，本发明方法的定位结果如图14所示。

附图说明

图1本发明的技术方案功能框图。

图2活动形状模型方法流程图。

图3原始LBP编码值计算的示意图。

图4扩展LBP算子的例子。

图5解梅和魏云龙的专利《一种基于改进的ASM算法的人脸特征定位方法》（中国专利公开号：201110167408.4，2011）方法流程图。

图6Jin Q,Zhao J Y和Zhang Y的论文“Facial feature extraction with a depthAAM algorithm”（International Conference on Fuzzy Systems andKnowledge Discovery,Sichuan,China.2012:1792-1796）中使用深度AAM算法的面部特征提取方法流程图。

图7训练样本对齐模块流程图。

图8活动形状模型建立模块流程图。

图9特征点j的统计LBP直方图的选取、计算示意图。

图10特征点搜索方案流程图。

图11使用统计LBP直方图进行搜索示意图。

图12瞳孔定位结果例子。

图13活动形状模型的初始位置例子。

图14面部特征点定位结果例子。

具体实施方式

使用基于统计学的ASM模型，该模型分为两个阶段：模型训练阶段和图像搜索阶段。在模型训练阶段，选取足够数量的正常人脸图像和不同光照下的人脸图像作为训练集，手工标定好每幅训练图像准确的特征点位置，并建立活动形状模型和局部纹理模型。在图像搜索阶段，首先确定目标图像准确的瞳孔位置，以初始化活动形状模型的初始位置，然后通过迭代使模型拟合到准确的目标特征点位置。本发明的技术方案功能框图如图1所示。

(1)活动形状模型

活动形状模型(Active Shape Model,ASM)是由英国曼彻斯特大学的Cootes等在1994年提出，它与活动轮廓模型(Active Contour Model，ACM)类似，是一种基于统计学的模型。ASM方法应用广泛，是一种有效的面部特征点定位方法。其主要思想是：首先在训练集合的每幅图像上手工标定准确的面部特征点，以此作为形状向量；将所有形状向量的集合作为形状样本集，并据此样本集建立活动形状模型和局部纹理模型；然后利用建立好的活动形状模型和局部纹理模型对目标图像进行定位，模型通过不断地变形来对目标形状进行拟合，直到收敛到准确的面部特征点位置。ASM方法的流程如图2所示。

ASM方法对模型的初始位置很敏感，如果初始模型中特征点位置靠近实际特征点位置，ASM方法能快速且准确地找到所有特征点；如果初始位置远离实际特征点位置，ASM方法常给出错误的定位结果。此外，ASM方法受光照影响较大，如果光照改变了图像灰度，则ASM方法也会定位错误。

(2)LBP算子

Ojala等提出的局部二值模式(Local Binary Patterns，LBP)算子，是一种图像纹理分析的有效方法。该方法可有效地描述纹理，且对光照变化具有鲁棒性，近年逐渐应用于人脸识别。原始LBP算子用一个3×3的滑动窗口扫描图像，图像中的每个像素点都对应一个3×3邻域，该像素点位于窗口中心。将邻域内的所有像素点与中心点比较，若该点像素点小于中心点，则取0，否则取1；然后将得到的0、1数字串顺时针连接，形成一个二进制的编码，即为该中心点的原始LBP编码值。原始LBP编码值计算的示意图如图3所示。

对于给定中心点(x_c,y_c)，设v_c为中心点的灰度值，v_i为窗口中其余各点的灰度值，则该点的原始LBP十进制编码值为

$\mathrm{LBP}(x_c,y_c)=\underset{i=0}{\overset{7}{\mathrm{\Σ}}}\mathrm{fun}(v_i-v_c)2^i---\left(1\right)$

其中阶跃函数fun(x)定义为

$\mathrm{fun}\left(x\right)=\left\{\begin{array}{cc}1,& x\&GreaterEqual;0\\ 0,& x<0\end{array}\right.$

当光照变化时，整个邻域内的像素灰度都随之变化，但邻域内的相对灰度信息保持不变，即原始LBP编码值不变，由此可见，LBP对光照有一定鲁棒性。但受原始LBP算子尺寸过小的限制，其无法描述大尺寸图像的特征。

Ojala T,Pietikainen M和Maenpaa T在论文“Multiresolution gray-scale androtation invariant texture classification with local binary patterns”(IEEE Transationon Pattern Analysis and Machine Intelligence,2002,24:971-987)中将原始LBP算子扩展为不同半径的圆形区域。图4给出了扩展LBP算子(记为LBP_P,R)的例子，其表示在半径为R的圆上均匀地选择P个像素点。设(x_i,y_i)为邻域点的坐标，(x_c,y_c)为中心点的坐标，则其邻域点坐标为

$(x_i,y_i)=(x_c+R\sin \left(\frac{2\mathrm{\&pi;i}}{P}\right),y_c+R\cos \left(\frac{2\mathrm{\&pi;i}}{P}\right))---\left(3\right)$

人脸特征点定位多用LBP统计直方图描述面部特征的局部纹理信息。但LBP_P,R直方图的数量有2^P种，如此高的维数无法在实际中应用。LBP算子的进一步扩展是用统一模式，它可降低LBP纹理描述的维数，即减少LBP直方图数量。经LBP_P,R算子编码后，得到的二进制串属于统一模式类的充分条件是：该二进制串最多只包含两次从1到0或者从0到1的变化。而不属于统一模式类的编码都是混合模式类。例如，00000000，11111001等最多只包含两次变化，属于统一模式类；而10010001和01110110等包含了四次变化，属于混合模式类。LBP_P,R总共有(P(P-1)+2)种统一模式类，而其余为混合模式类，于是统一LBP算子(记为

)可重新定义为

其中

$g\left(\mathrm{LB}{P}_{P,R}^{u2}\right)=|\mathrm{fun}(v_{P-1}-v_c)-\mathrm{fun}(v_0-v_c)|+\underset{i=1}{\overset{P-1}{\mathrm{\&Sigma;}}}|\mathrm{fun}(v_i-v_c)-\mathrm{fun}(v_{i-1}-v_c)|---\left(5\right)$

这里，v_c为中心点(x_c,y_c)的灰度值，v_i为其余各点的灰度值，g(x)表示该二进制串的跳变次数。

显然，统一

算子的模式种类由2^P种变成了P(P-1)+3种，数量得到极大减少。经统一

算子计算出的模式是人脸图像中最重要的模式，它能够描述图像中的大部分有效信息。

(3)LBP直方图

直方图是LBP算子描述局部纹理特征的方法，其定义为

$H_i=\underset{x,y}{\mathrm{\&Sigma;}}I\{L(x,y)=i\},i=\mathrm{0,1},\&CenterDot;\&CenterDot;\&CenterDot;Q-1---\left(6\right)$

其中，Q为LBP模式的数量，即直方图的数量；L(x,y)表示图像I在点(x,y)处的LBP编码值；I(X)定义为

LBP直方图可较好地描述图像的纹理特征，但如果直接对较大的区域进行统计，会损失空间位置信息，且极易丢失局部微观信息。为了解决该问题，一般将较大区域进行划分，使其变为若干个较小子区域的组合，然后统计每个子区域的LBP直方图，最后将每个子区域的直方图组合成原区域的直方图，即为原区域的LBP特征。

用LBP直方图描述纹理特征时，需要用LBP直方图相似性度量方法来确定两个特征点的相似度。假设H_i是某个特征点的LBP直方图，H_j为另一特征点的LBP直方图，则H_i与H_j的相似性度量为

${\mathrm{\&chi;}}^2\left(H_{i,}{H}_j\right)=\underset{k}{\mathrm{\&Sigma;}}\frac{{[H_i\left(k\right)-H_j\left(k\right)]}^2}{(H_i\left(k\right)+H_j\left(k\right))}---\left(8\right)$

χ²(H_i,H_j)的值越小，说明两个直方图越相似。

现有技术一的技术方案

解梅和魏云龙在公开号为201110167408.4的中国发明专利“一种基于改进的ASM算法的人脸特征定位方法”中，提出了一种对面部表情鲁棒性较高的特征点定位方法。该发明针对传统ASM方法在人脸包含表情变化情况下其特征点定位遇到的困难，通过将人脸面部特征按变化相关度进行区域划分，形成上、下两个形状区域，对这两个区域独立进行统计形状模型和局部灰度模型的建模。特征点匹配搜索过程是通过引入一个能量函数来对上、下模型生成实例的综合形状进行误差约束，最终得到准确的特征点定位结果。该方法流程图如图5所示。

现有技术一的缺点

技术一的主要缺点是：(1)仅用AdaBoost算法进行人脸检测，其人脸位置的初始模型不够准确；(2)对于光照等情况，其鲁棒性较低。

现有技术二的技术方案

Jin Q,Zhao J Y和Zhang Y在论文“Facial feature extraction with a depth AAMalgorithm”（International Conference on Fuzzy Systems and Knowledge Discovery,Sichuan,China.2012:1792-1796）中提出了深度AAM算法。该方法首先利用人脸检测算法定位出头部并估计头部姿态，然后利用头部位置和姿态初始化AAM的全局形状姿态，这保证模型能够拟合到正确位置；在AAM的全局能量函数中自适应地调节RGB信息和深度信息的权重，从而提高算法的鲁棒性和定位精度；此外，该方法还使用了图像金字塔算法和反向组合算法来加快迭代速度。该方法流程图如图6所示。

现有技术二的缺点

技术二的方法只解决了头部姿态对面部特征点定位的影响，其对光照的鲁棒性仍然较低。

本发明技术方案的详细阐述

本发明所要解决的技术问题

用统一模式的LBP算子对传统ASM算法进行改进，进一步提高传统ASM算法的脸部特征定位精度及对光照的鲁棒性。

本发明提供的完整技术方案

ASM模型是一种基于统计学的模型，它在使用中分为两个阶段：模型训练阶段和图像搜索阶段。在模型训练阶段，选取足够数量的正常人脸图像和不同光照下的人脸图像作为训练集，手工标定好每幅训练图像准确的特征点位置，并建立活动形状模型和局部纹理模型。在图像搜索阶段，首先确定目标图像准确的瞳孔位置，以初始化活动形状模型的初始位置，然后通过迭代使模型拟合到准确的目标特征点位置。本发明的技术方案如图1所示。

模型训练阶段

在模型训练阶段，从人脸数据库中选择N幅训练图像作为输入。CAS-PEAL-R1人脸图像数据库由中国科学院计算技术研究所创建，其包含了1040名中国人共99450幅头肩部图像，所有图像在专门的采集环境中采集，涵盖了姿态、表情、饰物和光照4种主要变化条件。IMM人脸库由丹麦技术大学创建，它包含了240张不同姿态、表情、光照的人脸图像。本发明从IMM与CAS-PEAL-R1人脸库中选取了200幅正常人脸图像和200幅不同光照下的人脸图像，则训练图像数量N=400。手工标定好每幅训练图像准确的68(n=68)个特征点的位置，记录每个特征点的坐标，合成N个形状向量X_u，即

X_ui＝(x_i0,y_i0,x_i1,y_i1,…,x_ik,y_ik,…,x_i(n-1),y_i(n-1))^T     (9)其中，x_ik为第i幅训练图像中特征点k的横坐标，y_ik为第i幅图像中特征点k的纵坐标。

(1)训练样本对齐模块：训练样本对齐模块的处理步骤如下：

该模块的输入为手工标定好的未对齐的N个形状向量X_u。这些向量不能体现出人脸形状变化的规律，因此，需要对这N个形状向量进行对齐，输出对齐后的形状向量X。图7为训练样本对齐模块流程图。

训练样本对齐的具体步骤如下：

(a)计算权重矩阵W。

首先，计算第i个形状中特征点m与特征点l的距离R_iml

$R_{\mathrm{iml}}=\sqrt{{(x_{\mathrm{im}}-x_{\mathrm{il}})}^2+{(y_{\mathrm{im}}-y_{\mathrm{il}})}^2}-----\left(10\right)$

其中，(x_im,y_im)为第i个形状中特征点m的坐标，(x_il,y_il)为第i个形状中特征点l的坐标。

然后，根据计算特征点m与特征点l距离的方差V_Rml，并计算特征点m的加权值w_m

$\left\{\begin{array}{c}\stackrel{\&OverBar;}{R_{\mathrm{iml}}}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{R}_{\mathrm{iml}}\\ V_{\mathrm{Rml}}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}(R_{\mathrm{iml}}-\stackrel{\&OverBar;}{R_{\mathrm{iml}}})\end{array}\right.---\left(11\right)$

$w_m={\left(\underset{l=0}{\overset{n-1}{\mathrm{\&Sigma;}}}{V}_{\mathrm{Rml}}\right)}^{-1}---\left(12\right)$

其中，R_iml表示第i个形状中特征点m与特征点l的距离，w_m表示特征点m的稳定程度。

最后，以w_m为对角线计算对角矩阵W，并将其作为权重矩阵。

$W=\left(\begin{array}{ccccc}{w}_0& 0& 0& 0& 0\\ 0& w_1& 0& 0& 0\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ 0& 0& 0& w_{m-2}& 0\\ 0& 0& 0& 0& w_{m-1}\end{array}\right)---\left(13\right)$

(b)选取X_u中的第一个形状向量作为初始平均形状

(c)将每个形状向量X_ui向平均形状对齐。首先根据式(14)计算仿射变换所需的四个参数：缩放尺度s、旋转角度θ、平移向量t，平移向量t包括x轴平移量t_x，y轴平移量t_y

$M^{\&prime;}=\left(\begin{array}{c}{a}_x\\ a_y\\ t_x\\ t_y\end{array}\right)={\left(B^T{W}^T\mathrm{WB}\right)}^{-1}{B}^T{W}^{T}W\stackrel{\&OverBar;}{X_u}---\left(14\right)$

其中， $B=\left(\begin{array}{cccc}{x}_{\mathrm{ui}0}& -y_{\mathrm{ui}0}& 1& 0\\ y_{\mathrm{ui}0}& x_{\mathrm{ui}0}& 0& 1\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ x_{\mathrm{ui}(n-1)}& -y_{\mathrm{ui}(n-1)}& 1& 0\\ y_{\mathrm{ui}(n-1)}& x_{\mathrm{ui}(n-1)}& 0& 1\end{array}\right),\left\{\begin{array}{c}{a}_x=s\cos \mathrm{\&theta;}\\ a_y=s\sin \mathrm{\&theta;}\end{array}\right.$

根据求得的四个仿射变换参数，通过式(15)对每一个形状向量X_ui进行平移、旋转和缩放操作，使它们向平均形状对齐。

$X_{\mathrm{ui}}^{\&prime;}=M(\mathrm{\&theta;},s)X_{\mathrm{ui}}+t=\left(\begin{array}{cccc}{x}_{\mathrm{ui}0}& -y_{\mathrm{ui}0}& 1& 0\\ y_{\mathrm{ui}0}& x_{\mathrm{ui}0}& 0& 1\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ x_{\mathrm{ui}(n-1)}& -y_{\mathrm{ui}(n-1)}& 1& 0\\ y_{\mathrm{ui}(n-1)}& x_{\mathrm{ui}(n-1)}& 0& 1\end{array}\right)\left(\begin{array}{c}s\cos \mathrm{\&theta;}\\ s\sin \mathrm{\&theta;}\\ t_x\\ t_y\end{array}\right)----\left(15\right)$

(d)计算对齐后所有形状向量的平均形状。

$\stackrel{\&OverBar;}{X_u^{\&prime;}}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{X}_{\mathrm{ui}}^{\&prime;}---\left(16\right)$

(e)重复上述(c)、(d)两步，直到两个相邻的平均形状收敛。在对齐结束后，输出对齐后的形状向量X。

判断两相邻平均形状收敛的方法，即将两个平均形状代入式(14)，求出s、θ、t，如果|s-1|<0.001，|θ|<0.001π/180，|t|<0.01这三个条件同时满足，则两个平均形状收敛；否则，不收敛。

(2)活动形状模型建立模块：活动形状模型建立模块的处理步骤如下：

该模块的输入为对齐后的形状向量X，利用对齐后的数据进行活动形状建模，建模结束后，输出活动形状模型。活动形状模型建立模块流程图如图8所示。具体步骤如下：

(a)计算出对齐后的N个形状向量的平均值，即平均形状

$\stackrel{\&OverBar;}{X}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{X}_i----\left(17\right)$

(b)求出N个形状向量的协方差矩阵

$F=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}(X_i-\stackrel{\&OverBar;}{X}){(X_i-\stackrel{\&OverBar;}{X})}^T---\left(18\right)$

(c)求出该协方差阵的特征值λi(从大到小排序)及其对应的特征向量p_i，其中，i=0,1,…,2n-1。

(d)利用式(19)选取前k个最大的特征值，即前k个特征值满足

$\frac{\underset{t=0}{\overset{k}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_t}{\underset{t=0}{\overset{2n-1}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_t}>0.95---\left(19\right)$

并将前k个特征值对应的特征向量构成主成分矩阵P=(p₀,p₁,…,p_k-1)。

(e)利用式(20)构建活动形状模型。

$X=\stackrel{\&OverBar;}{X}+\mathrm{Pb}---\left(20\right)$

其中，

为平均形状，P为主成分矩阵，b为t×1的形状参数。如果b的变化太大，那么式(20)计算出的新形状有可能不是人脸形状，因此需要利用式(4)对b进行约束。

$-3\sqrt{{\mathrm{\&lambda;}}_i}\&le;b_i\&le;3\sqrt{{\mathrm{\&lambda;}}_i}---\left(4\right)$ 其中，λ_i为协方差矩阵的特征值。如果，则将b_i约束为；如果b_i>

则将b_i约束为

(3)局部纹理模型建立模块：局部纹理模型建立模块的处理步骤如下：

该模块的输入为未对齐的N个形状向量X_u以及它们对应的N幅训练图像；建模之后，输出每个特征点的平均LBP直方图，即局部纹理模型。其具体步骤如下：

(a)在第i幅图像的特征点j处，以该特征点为中心选取一个25×25的正方形(该正方形的两边分别平行于x方向和y方向)，并将该正方形划分为4个子区域，如图9所示。

(b)利用

算子计算子区域中每个点的LBP编码值，然后根据式(6)统计每个子区域的LBP直方图，并根据图9所示的顺序将每个子区域的LBP直方图进行级联，作为该特征点的LBP直方图。

(c)使用相同方法计算出N幅训练图像特征点j的LBP直方图，并求取它们的平均直方图作为特征点j的LBP直方图，即特征点j的局部纹理模型。

计算出局部纹理模型，则目标图像中某个特征点与活动形状模型中特征点j的相似度可以用下式表示：

${\mathrm{\&chi;}}^2(H_{\mathrm{ij}},\stackrel{\&OverBar;}{H_i})=\underset{k}{\mathrm{\&Sigma;}}\frac{{(H_{\mathrm{ij}}\left(k\right)-\stackrel{\&OverBar;}{H_i}\left(k\right))}^2}{[H_{\mathrm{ij}}\left(k\right)+\stackrel{\&OverBar;}{H_i}\left(k\right)]}---\left(22\right)$

其中，H_ij为目标图像中特征点i邻域中第j个位置的LBP直方图，

为局部纹理模型中特征点i的平均直方图。

值越小，说明目标特征点与特征点j的相似度越大，所以在图像搜索过程中可以据此来确定特征点j的最佳匹配位置。

图像搜索阶段

在图像搜索阶段，输入待检测的目标图像，并给定活动形状模型一个初始位置，使其覆盖在目标图像上。由于活动形状模型的初始位置不可能处于准确的特征点位置，因此需要通过不断拟合，直至其收敛至准确的目标特征点处。

(1)瞳孔定位模块

该模块的输入为目标图像，本发明使用Yefei Chen和Jianbo Su在论文“Fast eyelocalization based on a new Haar-like feature”（The10th World Congress onIntelligent Control and Automation,Beijing,China.2012,4825-4830）中提出的基于新的类Haar特征的快速人眼定位方法来确定瞳孔位置。该模块的输出为目标图像左瞳孔位置(x₁,y₁)和右瞳孔位置(x₂,y₂)。

(2)特征点拟合模块

该模块的输入为目标图像、目标图像左右瞳孔的位置坐标(x₁,y₁)和(x₂,y₂)、活动形状模型以及局部纹理模型，其功能是根据局部纹理模型和式(4)所示的直方图相似性的判断方法，通过反复迭代，使得活动形状模型收敛至准确的目标特征点位置。方案流程如图10所示。

特征点拟合模块的主要步骤如下：

(a)计算活动形状模型的初始位置。首先取出主动形状模型中左右瞳孔的坐标，分别为(x₁′,y₁′)和(x₂′,y₂′)，并利用式(23)计算出缩放尺度s和旋转角度θ

$\left\{\begin{array}{c}s=\frac{d}{d^{\&prime;}}=\frac{\sqrt{{(x_1-x_2)}^2+{(y_1-y_2)}^2}}{\sqrt{{(x_1^{\&prime;}-x_2^{\&prime;})}^2+{(y_1^{\&prime;}-y_2^{\&prime;})}^2}}\\ \mathrm{\&theta;}=\arctan \frac{y_2-y_1}{x_2-x_1}\end{array}\right.---\left(23\right)$

其中，(x₁,y₁)和(x₂,y₂)分别为目标图像中左右瞳孔坐标。

令平移向量t=0，根据式(15)对活动形状模型进行缩放和旋转，得到一个临时形状，且临时形状左瞳孔的坐标为(x₁″,y₁″)。然后计算出平移向量t＝(t_x,t_y)^T＝(x_1-x₁″,y₁-y₁″)^T，且令s=0、θ=0，继续用式(15)对临时形状进行平移，得到模型的最终初始位置x_t。

(b)初始化姿态参数s=0、θ=0、t=0以及形状参数b=0。

(c)每次循环开始，模型对应的初始位置为X＝M(s,θ)x+t，其中，M(s,θ)的具体形式见式(15)，x＝x_t+Pb。

(d)计算每个特征点的新位置。首先将当前循环活动形状模型的初始位置覆盖在图像上，对于模型中第j个特征点，沿着特征点法线方向(即沿前一个特征点与后一个特征点连线的垂直方向)在该特征点左右两侧各取12点，加上特征点j构成一个长度为25的搜索邻域；然后，对于搜索邻域中的每一个点，以该点为中心选取一个25×25的正方形(该正方形的两边分别平行于x方向和y方向)，并将该正方形划分为4个子区域，根据式(6)计算每个子区域的统计LBP直方图，将每个子区域的LBP直方图进行级联，作为该点的统计LBP直方图；最后，利用式(4)计算搜索邻域中每个点的统计LBP直方图与特征点j的平均直方图的相似度，选取与平均直方图最接近的直方图所对应的点作为特征点j的最佳更新位置。用统计LBP直方图进行特征点最佳位置搜索的示意图如图11所示，从图11可以看出，特征点的更新必然会产生一个位移，所有特征点的位移可组成向量dX=(dX₀,dX₁,dX₂,…,dX_n-1)。

(e)更新姿态参数(包括缩放尺度s，旋转角度θ，平移向量t)及形状参数b。此次循环的初始位置为X=M(s,θ)x+t，将X向X+dX对齐，得到新的姿态参数t+dt，s(1+ds)，θ+dθ，从而可以得到如下公式：

$M(s(1+\mathrm{ds}),\mathrm{\&theta;}+\mathrm{d\&theta;})(\stackrel{\&OverBar;}{X}+\mathrm{Pb}+\mathrm{dx})+t+\mathrm{dt}=X+\mathrm{dX}---\left(24\right)$

根据式(25)求出dx

$\left\{\begin{array}{c}y=M(s,\mathrm{\&theta;})(\stackrel{\&OverBar;}{X}+\mathrm{Pb})+\mathrm{dX}-(t+\mathrm{dt})\\ \mathrm{dx}=M\{{\left[s(1+\mathrm{ds})\right]}^{-1},-(\mathrm{\&theta;}+\mathrm{d\&theta;})\}y-(\stackrel{\&OverBar;}{X}+\mathrm{Pb})\end{array}\right.---\left(25\right)$

最后根据式(26)～(27)求出形状参数b的变化量db，即

$x=\stackrel{\&OverBar;}{X}+\mathrm{Pb}+\mathrm{dx}\&ap;\stackrel{\&OverBar;}{X}+P(b+\mathrm{db})---\left(26\right)$

db≈P^Tdx        (27)得到参数的变化量之后，利用式(28)对参数进行更新

$\left\{\begin{array}{c}t\&RightArrow;t+w_t\mathrm{dt}\\ \mathrm{\&theta;}\&RightArrow;\mathrm{\&theta;}+w_{\mathrm{\&theta;}}\mathrm{d\&theta;}\\ s\&RightArrow;s(1+w_s\mathrm{ds})\\ b\&RightArrow;b+w_b\mathrm{db}\end{array}\right.----\left(28\right)$

其中，w_t、w_θ、w_s和w_b为控制参数变化的权重。另外，更新b后，需用式(4)对其进行调整。

(f)根据更新后的形状参数b，得到新形状x=x_t+Pb，若该形状与原形状的差别足够小(即db小于阈值0.0001)，或者循环次数已经达到最大次数(40次)，则搜索完成；否则，回到第(c)步，继续进行处理。

本发明技术方案带来的有益效果：

模型训练阶段，为每个特征点建立局部纹理模型时采用边长为25的正方形；图像搜索阶段，搜索邻域的长度为25个像素点，计算搜索邻域中每个点的LBP直方图时同样采用边长为25的正方形，权重w_t、w_θ、w_s和w_b的值都为1。用YefeiChen和Jianbo Su在论文“Fast eye localization based on a new Haar-like feature”（The10th World Congress on Intelligent Control and Automation,Beijing,China.2012,4825-4830）中所述的方法定位瞳孔位置，结果如图12所示。图13为用瞳孔位置得到的活动形状模型的初始位置的例子。在得到活动形状模型的初始位置后，用局部纹理模型与直方图相似性判别公式，对活动形状模型的特征点进行调整，使其收敛至准确的特征点位置，本发明方法的定位结果如图14所示。

为进一步验证本发明的准确性和鲁棒性，从IMM人脸库与CAS-PEAL-R1人脸库中，选取完全不同于训练集的100张正常条件下的图像及100张不同光照下的图像作为测试集，利用该测试集对本发明方法进行测试。测试过程中，只要准确定位的特征点的数量超过90%，就认为该图像定位准确。测试结果如表1所示。由表1可见，对于测试集中的100张正常人脸图像，本发明方法能够准确定位的数量为96张；对于测试集中100张不同光照变化的人脸图像，能够准确定位的数量为75张。本发明综合定位准确率为85.5%。

表1 定位准确率测试结果

选取17个有代表意义的特征点，如表2所示，并据式(29)求出光照变化图像中这些特征点的平均欧式距离误差，即

$\mathrm{Err}=\frac{1}{N_c}\underset{i=1}{\overset{N_c}{\mathrm{\&Sigma;}}}\sqrt{{(x_i-x_i^{\&prime;})}^2+{(y_i-y_i^{\&prime;})}^2}$

其中，N_c为测试所需的光照图像的数量，(x_i,y_i)为标定的特征点坐标，(x_i′,y_i′)为本发明算法得到的相应特征点。表3给出了这17个特征点的平均欧式距离误差，从表3可以看出，本发明提出的方法能够有效地光照的干扰，定位鲁棒性较高。

表2 有代表意义的特征点

表3 光照图像的平均欧式距离误差

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。

本发明涉及的缩略语和关键术语定义：

AAM：Active Appearance Model，活动表观模型。

ASM：Active Shape Model，活动形状模型。

LBP：Local Binary Patterns，局部二值模式。

Claims (3)

1.一种面部图像中特征点的定位方法，其特征在于：包括模型训练阶段和图像搜索阶段；在模型训练阶段，选取足够数量的正常人脸图像和不同光照下的人脸图像作为训练集合，手工标定好每幅训练图像准确的特征点位置，并建立活动形状模型和局部纹理模型，该模型训练阶段仅需执行一次；在图像搜索阶段，首先确定目标图像准确的瞳孔位置，以初始化活动形状模型的初始位置，然后通过迭代使模型拟合到准确的目标特征点位置。

2.根据权利要求1所述的一种面部图像中特征点的定位方法，其特征在于：模型训练阶段的处理方法为：从CAS-PEAL-R1人脸数据库中选择N幅训练图像作为输入；从IMM与CAS-PEAL-R1人脸库中选取了200幅正常人脸图像和200幅不同光照下的人脸图像，构成人脸图像训练集合，则训练图像数量N=400；手工标定好每幅训练图像准确的n个特征点的位置，n=68，记录每个特征点的坐标，合成N个形状向量X_u，即

X_ui＝(x_i0,y_i0,x_i1,y_i1,…,x_ik,y_ik,…,x_i(n-1),y_i(n-1))^T，i=0,1,…,N-1其中，x_ik、y_ik分别为第i幅训练图像中特征点k的横、纵坐标，0≤k<n；所述的训练样本对齐模块的输入为形状向量X_u，输出端与活动形状模型建立模块的输入端相连接，所述的活动形状模型建立模块输出端与特征点拟合模块相连接，所述的局部纹理模型建立模块的输入为形状向量X_u，输出端与特征点拟合模块相连接；形状向量X_u经过训练样本对齐模块处理后，输出对齐后的形状向量X，活动形状模型建立模块根据对齐后的形状向量X生成活动形状模型，供特征点拟合模块使用；局部纹理模型建立模块根据形状向量X_u生成局部纹理模型，供特征点拟合模块使用，所述的训练样本对齐模块、活动形状模型建立模块、局部纹理模型建立模块仅在活动形状和局部纹理模型建立时执行一次；

训练样本对齐模块的处理步骤如下：

(a)计算权重矩阵W首先，计算第i个形状中特征点m与特征点l的距离R_iml

$R_{\mathrm{iml}}=\sqrt{{(x_{\mathrm{im}}-x_{\mathrm{il}})}^2+{(y_{\mathrm{im}}-y_{\mathrm{il}})}^2}$

其中，(x_im,y_im)为第i个形状中特征点m的坐标，(x_il,y_il)为第i个形状中特征点l的坐标；

然后，根据计算特征点m与特征点l距离的方差V_Rml，并计算特征点m的加权值w_m

$\left\{\begin{array}{c}\stackrel{\&OverBar;}{R_{\mathrm{iml}}}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{R}_{\mathrm{iml}}\\ V_{\mathrm{Rml}}=\frac{1}{n}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}(R_{\mathrm{iml}}-\stackrel{\&OverBar;}{R_{\mathrm{iml}}})\end{array}\right.$

$w_m={\left(\underset{l=0}{\overset{n-1}{\mathrm{\&Sigma;}}}{V}_{\mathrm{Rml}}\right)}^{-1}$

其中，R_iml表示第i个形状中特征点m与特征点l的距离，w_m表示特征点m的稳定程度；

最后，以w_m为对角线计算对角矩阵W，并将其作为权重矩阵；

$W=\left(\begin{array}{ccccc}{w}_0& 0& 0& 0& 0\\ 0& w_1& 0& 0& 0\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ 0& 0& 0& w_{m-2}& 0\\ 0& 0& 0& 0& w_{m-1}\end{array}\right)$

(b)选取X_u中的第一个形状向量作为初始平均形状

(c)将每个形状向量X_ui向平均形状对齐；首先根据式(1)计算仿射变换所需的四个参数：缩放尺度s、旋转角度θ、平移向量t，平移向量t包括x轴平移量t_x，y轴平移量t_y，

$M^{\&prime;}=\left(\begin{array}{c}{a}_x\\ a_y\\ t_x\\ t_y\end{array}\right)={\left(B^T{W}^T\mathrm{WB}\right)}^{-1}{B}^T{W}^{T}W\stackrel{\&OverBar;}{X_u}---\left(1\right)$ 其中， $B=\left(\begin{array}{cccc}{x}_{\mathrm{ui}0}& -y_{\mathrm{ui}0}& 1& 0\\ y_{\mathrm{ui}0}& x_{\mathrm{ui}0}& 0& 1\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ x_{\mathrm{ui}(n-1)}& -y_{\mathrm{ui}(n-1)}& 1& 0\\ y_{\mathrm{ui}(n-1)}& x_{\mathrm{ui}(n-1)}& 0& 1\end{array}\right),\left\{\begin{array}{c}{a}_x=s\cos \mathrm{\&theta;}\\ a_y=s\sin \mathrm{\&theta;}\end{array}\right.$

根据仿射变换的四个参数M′，通过式(2)对每一个形状向量X_ui进行平移、旋转和缩放操作，使它们向平均形状对齐；

$X_{\mathrm{ui}}^{\&prime;}=M(\mathrm{\&theta;},s)X_{\mathrm{ui}}+t=\left(\begin{array}{cccc}{x}_{\mathrm{ui}0}& -y_{\mathrm{ui}0}& 1& 0\\ y_{\mathrm{ui}0}& x_{\mathrm{ui}0}& 0& 1\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ x_{\mathrm{ui}(n-1)}& -y_{\mathrm{ui}(n-1)}& 1& 0\\ y_{\mathrm{ui}(n-1)}& x_{\mathrm{ui}(n-1)}& 0& 1\end{array}\right)\left(\begin{array}{c}s\cos \mathrm{\&theta;}\\ s\sin \mathrm{\&theta;}\\ t_x\\ t_y\end{array}\right)---\left(2\right)$

(d)计算对齐后所有形状向量的平均形状；

$\stackrel{\&OverBar;}{X_u^{\&prime;}}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{X}_{\mathrm{ui}}^{\&prime;}$

(e)重复(c)、(d)两步，直到两个相邻的平均形状收敛；在对齐结束后，输出对齐后的形状向量X；

判断两相邻平均形状收敛的方法，即将两个平均形状代入式(1)，求出s、θ、t，如果|s-1|<0.001，|θ|<0.001π/180，|t|<0.01这三个条件同时满足，则认为两个平均形状收敛；否则，认为尚未收敛；

活动形状模型建立模块的处理步骤如下：

该模块的输入为对齐后的形状向量X，用对齐后的数据进行活动形状建模，建模结束后，输出活动形状模型，具体步骤如下：

(a)计算出对齐后的N个形状向量的平均值，即平均形状

$\stackrel{\&OverBar;}{X}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{X}_i$

(b)求出N个形状向量的协方差矩阵

$F=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}(X_i-\stackrel{\&OverBar;}{X}){(X_i-\stackrel{\&OverBar;}{X})}^T$

(c)求出该协方差阵的特征值λ_i，从大到小排序及其对应的特征向量p_i，i=0,1,…,2n-1；

(d)根据条件

选取前k个最大的特征值，并将前k个特征值对应的特征向量构成主成分矩阵P=(p₀,p₁,…,p_k-1)；

(e)利用式(3)构建活动形状模型；

$X=\stackrel{\&OverBar;}{X}+\mathrm{Pb}---\left(3\right)$

其中，

为平均形状，P为主成分矩阵，b为t×1的形状参数；如果b的变化太大，式(3)计算出的新形状有可能不是人脸形状，因此需用式(4)对b进行约束；

$-3\sqrt{{\mathrm{\&lambda;}}_i}\&le;b_i\&le;3\sqrt{{\mathrm{\&lambda;}}_i}---\left(4\right)$ 其中，λ_i为协方差矩阵F的特征值；如果

，则将b_i约束为

如果b_i>

则将b_i约束为

局部纹理模型建立模块的处理步骤如下：

该模块的输入为未对齐的N个形状向量X_u以及它们对应的N幅训练图像；建模之后，输出每个特征点的平均LBP直方图，即局部纹理模型；其具体步骤如下：

(a)在第i幅图像的特征点j处，以该特征点为中心选取一个25×25的正方形，该正方形的两边分别平行于x方向和y方向,并将该正方形划分为4个子区域；

(b)利用统一算子计算子区域中每个点的LBP编码值，然后根据式(6)统计每个子区域的LBP直方图；

统一LBP算子，记为

定义为

其中，v_c为中心点(x_c,y_c)的灰度值，v_i为其余各点的灰度值，g(x)表示该二进制串的跳变次数， $g\left({\mathrm{LBP}}_{P,R}^{u2}\right)=|\mathrm{fun}(v_{P-1}-v_c)-\mathrm{fun}(v_0-v_c)|+\underset{i=1}{\overset{P-1}{\mathrm{\&Sigma;}}}|\mathrm{fun}(v_i-v_c)-\mathrm{fun}(v_{i-1}-v_c)|,$ $\mathrm{fun}\left(x\right)=\left\{\begin{array}{cc}1,& x\&GreaterEqual;0\\ 0,& x<0\end{array}\right.,$

直方图是LBP算子描述局部纹理特征的方法，其定义为

$H_i=\underset{x,y}{\mathrm{\&Sigma;}}I\{L(x,y)=i\},i=\mathrm{0,1},\&CenterDot;\&CenterDot;\&CenterDot;,Q-1---\left(6\right)$ 其中，Q为LBP模式的数量，即直方图的数量；L(x,y)表示图像I在点(x,y)处的LBP编码值；I(X)定义为

并根据在图像上由左至右，由上至下的顺序将每个子区域的LBP直方图进行级联，作为该特征点的LBP直方图；

(c)使用相同方法计算出N幅训练图像特征点j的LBP直方图，并求取它们的平均直方图作为特征点j的LBP直方图，即特征点j的局部纹理模型；

计算出局部纹理模型，则目标图像中某个特征点与活动形状模型中特征点j相似度可以用下式表示：

${\mathrm{\&chi;}}^2\left(H_{\mathrm{ij},}\stackrel{\&OverBar;}{H_i}\right)=\underset{k}{\mathrm{\&Sigma;}}\frac{{(H_{\mathrm{ij}}\left(k\right)-\stackrel{\&OverBar;}{H_i}\left(k\right))}^2}{(H_{\mathrm{ij}}\left(k\right)+\stackrel{\&OverBar;}{H_i}\left(k\right))}$

其中，H_ij为目标图像中特征点i邻域中第j个位置的LBP直方图，

为局部纹理模型中特征点i的平均直方图；

值越小，说明目标特征点与特征点j的相似度越大，所以在图像搜索过程中可以据此来确定特征点j的最佳匹配位置。

3.根据权利要求1所述的一种面部图像中特征点的定位方法，其特征在于：图像搜索阶段的步骤包括如下：

(1)瞳孔定位模块

(2)特征点拟合模块

该模块的输入为目标图像、目标图像左右瞳孔的位置坐标(x₁,y₁)和(x₂,y₂)、活动形状模型以及局部纹理模型，其功能是根据局部纹理模型和式(4)所示的直方图相似性的判断方法，通过反复迭代，使得活动形状模型收敛至准确的目标特征点位置；特征点拟合模块的主要步骤如下：

(a)计算活动形状模型的初始位置；首先取出主动形状模型中左右瞳孔的坐标，分别为(x₁′,y₁′)和(x₂′,y₂′)，并计算出缩放尺度s和旋转角度θ，

$\left\{\begin{array}{c}s=\frac{d}{d^{\&prime;}}=\frac{\sqrt{{(x_1-x_2)}^2+{(y_1-y_2)}^2}}{\sqrt{{(x_1^{\&prime;}-x_2^{\&prime;})}^2+{(y_1^{\&prime;}-y_2^{\&prime;})}^2}}\\ \mathrm{\&theta;}=\arctan \frac{y_2-y_1}{x_2-x_1}\end{array}\right.$

其中，(x₁,y₁)和(x₂,y₂)分别为目标图像中左右瞳孔坐标；

令平移向量t=0，根据式(15)对活动形状模型进行缩放和旋转，得到一个临时形状，且临时形状左瞳孔的坐标为(x₁″,y₁″)，然后计算出平移向量

t＝(t_x,t_y)^T＝(x₁-x₁″,y₁-y₁″)^T

且令s=0，θ=0，继续用式(2)对临时形状进行平移，得到模型的最终初始位置x_t；

(b)初始化姿态参数s=0、θ=0、t=0以及形状参数b=0；

(c)每次循环开始，模型对应的初始位置为X=M(s,θ)+t，其中，M(s,θ)的计算见式(15)，x=x_t+Pb；

(d)计算每个特征点的新位置，首先将当前循环活动形状模型的初始位置覆盖在图像上，对于模型中第j个特征点，沿着特征点法线方向，即沿前一个特征点与后一个特征点连线的垂直方向，在该特征点左右两侧各取12点，加上特征点j构成一个长度为25的搜索邻域；然后，对于搜索邻域中的每一个点，以该点为中心选取一个25×25的正方形，该正方形的两边分别平行于x轴方向和y轴方向；并将该正方形划分为4个子区域，根据式(6)计算每个子区域的统计LBP直方图，将每个子区域的LBP直方图进行级联，作为该点的统计LBP直方图；最后，利用式(4)计算搜索邻域中每个点的统计LBP直方图与特征点j的平均直方图的相似度，选取与平均直方图最接近的直方图所对应的点作为特征点j的最佳更新位置；特征点的更新必然会产生一个位移，所有特征点的位移可组成向量dX=(dX₀,dX₁,dX₂,…,dX_n-1)；

(e)更新姿态参数，包括缩放尺度s，旋转角度θ，平移向量t及形状参数b；此次循环的初始位置为X=M(s,θ)x+t，将X向X+dX对齐，得到新的姿态参数t+dt，s(1+ds)，θ+dθ，从而可以得到如下公式：

$M(s(1+\mathrm{ds}),\mathrm{\&theta;}+\mathrm{d\&theta;})(\stackrel{\&OverBar;}{X}+\mathrm{Pb}+\mathrm{dx})+t+\mathrm{dt}=X+\mathrm{dX}$

求出dx

$\left\{\begin{array}{c}y=M(s,\mathrm{\&theta;})(\stackrel{\&OverBar;}{X}+\mathrm{Pb})+\mathrm{dX}-(t-\mathrm{dt})\\ \mathrm{dx}=M\{{\left[s(1+\mathrm{ds})\right]}^{-1},-(\mathrm{\&theta;}+\mathrm{d\&theta;})\}y-(\stackrel{\&OverBar;}{X}+\mathrm{Pb})\end{array}\right.$

最后求出形状参数b的变化量db

$x=\stackrel{\&OverBar;}{X}+\mathrm{Pb}+\mathrm{dx}\&ap;\stackrel{\&OverBar;}{X}+P(b+\mathrm{db})$

db≈P^Tdx，

得到参数的变化量之后，利用

$\left\{\begin{array}{c}t\&RightArrow;t+w_t\mathrm{dt}\\ \mathrm{\&theta;}\&RightArrow;\mathrm{\&theta;}+w_{\mathrm{\&theta;}}\mathrm{d\&theta;}\\ s\&RightArrow;s(1+w_s\mathrm{ds})\\ b\&RightArrow;b+w_b\mathrm{db}\end{array}\right.$

对参数进行更新；其中，w_t、w_θ、w_s和w_b为控制参数变化的权重；另外，更新b后，需用式(4)对其进行调整；

(f)根据更新后的形状参数b，得到新形状x=x_t+Pb，若该形状与原形状的差别足够小即db小于阈值0.0001，或者循环次数已经达到最大次数40次，则搜索完成；否则，回到第(c)步，继续进行处理。

Images