融合多视角、多线索二维信息的人脸三维模型的建立方法
专利汇可以提供融合多视角、多线索二维信息的人脸三维模型的建立方法专利检索，专利查询，专利分析的服务。并且融合多视 角 、多线索二维信息的人脸三维模型的建立方法属于计算机人脸图像重建技术领域，其特征在于：首先用标定好的上下两台立体摄像机拍摄建模对象的人脸，得到人脸在无表情变化时从 正面 逐步转到侧面的视频序列；再用全自动的 姿态 估计 算法 得到序列中各时刻的精确的姿态参数；接着在模型和姿态二者的初始化之后，根据两类二维线索分别 抽取 三维信息：一类是基于图象间对应点进行立体重建，另一类是模型投影轮廓与图象轮廓的匹配、修正；最后，根据融合在多视角中提取出的、对应于不同二维线索的三维信息建立人脸模型。它得到的人脸模型有更细致、更精确的形状，并能做MPEG4(运动图像压缩编码组织标准4)兼容的人脸动画。，下面是融合多视角、多线索二维信息的人脸三维模型的建立方法专利的具体信息内容。
1.融合多视角、多线索二维信息的人脸三维模型的建立方法，含有二维对应创建人脸模 型的方法，其特征在于：它首先用与人脸上、下位置相对应的两台相同且外极线方向沿 垂直方向的摄像机在人脸无表情变化的情况下拍摄建模对象的人脸从正面逐步转约90度到 侧面的视频序列，再用两块相同的采集卡来采集上下两路同步视频并输入计算机中；然后， 从人脸的通用模型Faceg开始，根据获得的建模对象的二维信息逐步进行修改，最终得到特定 人的人脸模型Faces，其中修改的只是人脸上的三维结点列表P中的三维点坐标，而人脸上的 三角形面片列表F保持不变，即只改变模型形状P而模型的拓扑结构F不变；具体而言，它 依次含有如下步骤：
(1).输入立体视频序列IL，t和IR，t：
L，R对应于两台摄像机，t＝1，...，N，表示N帧不同时刻的图像：I则是由像素组成的二维 数组，每个像素包括R、G、B三个分量；t＝1的时刻对应于人脸正面姿态，通过人的判断手 工选择；
(2).姿态估计，即求解人头旋转过程中的姿态参数，即刚体变换系数：
用下述刚体变换公式表示上述图象序列中人脸的帧间运动：
Mt+1＝RtMt+Tt，t＝1，..，N-1，
Mt、Mt+1分别表示第t帧、第t+1帧时人脸上任一点在L摄像机坐标系中的坐标；Mt、 Mt+1是人脸上的同一位置在不同时刻的三维坐标，即是帧间三维对应点；用可靠三维点的运 动就可估计姿态参数；
{Rt，Tt，t＝1，...，N-1}是第t帧到第t+1帧姿态参数的旋转和平移分量；
姿态估计依次含有如下步骤：
(2.1)，起始时刻：先进行二维特征的检测和匹配，再通过立体重建得到起始时刻的三维 点：
首先，设置特征检测帧t0＝t，并在IL，t图上作二维特征检测，其步骤如下：
用已有的KLT算法计算IL，t图上特征点p(x＝i，y＝j)处的图像变化程度G(i，j)：
为待测点p处图像的梯度向量：
I(i+1，j)为该点处图像的灰度，余类推；
以待测点p为中心的窗口W的尺寸为(2*blockx+1)*(2*blocky+1)；然后再用已有的matlab 工具包求出G的最小奇异值sing_min，再选择sing_min很大的点p为检测出的特征点；把在 t时刻检测出的特征点记为{pL，t(k)}，L，t表示图象序号，K(t)是特征点总数，k为其标号， 1≤k≤K(t)；
其次，对同一时刻两个摄像机L，R拍摄同一场景得到的两幅图像即立体图对作二维特征 匹配；同一个三维点在不同图像中的投影都为二维点，它们之间互称为二维对应，立体图对 间的二维对应称为匹配：
对于IL，t图中某一点pL，t(k)＝(xL，t(k)yL，t(k))T，它的匹配点pR，t(k)＝(xR，t(k)yR，t(k))T应位 于一条由pL，t(k)决定的直线上，在这个一维空间中，选择使一个局部窗口W内总灰度差异最 小的点作为pL，t(k)的匹配点pR，t(k)，总灰度差异的表达式：
(i，j)是直线上的任一点；记时刻t的匹配结果为{pR，t(k)}，总数目为K(t)个， 每点都与pL，t(k)相对应，在下述公式表示的直线上搜索，使Diff(i，j)最小的位置：
其中：
fL，fR，uOL，vOL，uOR，vOR为 摄像机内参数，已知量；Rc，Tc＝[Tc，x Tc，y Tc，z]T分别为旋转和平移分量，属摄像机外参数， 已知量；
再次，在获得同一时刻立体图对间的匹配结果pL，t和pR，t后，通过立体重建得到用摄像机 坐标系表示的三维点Mm，t，m为摄像机L或R，Mm，t为t时刻三维点在m摄像机坐标系中的 坐标；
根据两摄像机的配置关系式：
ML＝RcMR+Tc，ML、MR分别是同一点在这两个摄像机坐标系中的坐标；以及摄像机 透视投影公式：
sm，即sL或sR，是待定的比例因子由稍后的公式求出，
得到三维点在摄像机L、R中的坐标ML，t和MR，t：
而
在姿态估计的余下部分，三维点都用L摄像机坐标系中的坐标来表示，因此省略下标L， 即将ML，t简记为Mt；
(2.2)，后继时刻：对t+1时刻的立体图对用KLT算法进行二维特征跟踪，再同样通过立 体重建得到各后继时刻的三维点；跟踪是指同一摄像机拍摄的不同时刻图像间的二维对应：
在二维特征跟踪时，设检测的特征点在图像Im，t(p)，m为L或R上，p(x＝i，y＝j)是任一特 征点，则其在图像Im，t+1上的跟踪结果为其中 G已如上述，而
记跟踪得到的特征点为{pm，t+1(k)}，m表示L或R，t+1为时刻：
对于摄像机L为：
对于摄像机R为：
根据任一对匹配的跟踪点计算三维点{Mt+1(k)|k＝1，..，K(t)}，如前所述，这里省略了下标 L；
其中
(2.3)，姿态参数初始化：根据两个时刻t和t+1之间的一组三维对应点{Mt(k)|k＝1，...，K(t)} 和{Mt+1(k)|k＝1，...，K(t)}，和表示可靠性的序列{bReliablek}，k＝1，...，K(t)从式 Mt+1(k)＝RtMt(k)+Tt，t＝1，...，N-1中用基于最小中值的鲁棒估计算法求解Rt，Tt：
设：定义并不断更新布尔数组{bReliableZk}为可靠度量序列，Yes表示可靠，No表示不可 靠，其中取值为Yes的项的标号为{ks}，s＝1，...，Kr(t)，Kr(t)是取值为Yes的项的总数，从其 中挑出Num_subset个子集，每个子集含有三对三维对应点，即 {Setn＝{N1，n，N2，n，N3，n}}，n＝1，...，Num_subset，1≤N1，n，N2，n，N3，n≤Kr(t)，每个三元子集中三维对 应点为： 和 可简记为 {(Mim，Mim’)|im＝0，1，2}，它表示三对不共线的三维对应点，Rt，Tt满足 Mim′＝RtMim+Tt|im＝0，1，2；然后从下列公式中求出R和T，略去下标t：
其中
是旋转的中心轴，称旋转轴，3*1的向量；
θ是绕的旋转角；
[r]x是根据3*1向量r＝[rx ry rz]T定义的反对称矩阵： 而 是单位矩阵；
平移分量T由下式估计出：
对每个子集Setn按上述公式估计出姿态参数，记为Rt，n，Tt，n；再对所有的可靠三维对应 {Mt(ks)，Mt+1(ks)}，s＝1，...，Kr(t)按下式计算对应点Mt(ks)，Mt+1(ks)与姿态参数Rt，n，Tt，n的符合 程度εn，s，该值越小，则两者越符合：
εn，s＝‖Mt+1(ks)-Rt，nMt(ks)-Tt，n‖2；
再从集合{Rt，n，Tt，n}，n＝1，...，Num_subset中选取Rt，Tt：
表示针对下标s取集合{εn，s|s＝1，...，Kr(t)}的中值，它对不同的n构成一个新的数组 而则是在这个中值数组中选择取最小值的那个n；
(2.4)，在t＝t+1时刻，设置：
是否需要重新调用特征检测和匹配模块的布尔标志量bReDetect；
是否需要进行姿态求精的布尔标志量bRefine；
(2.5)，当t＝N时，姿态估计结束；
(3).模型初始化和姿态初始化：
模型初始化是通用模型特定化方法的初始步骤，它从通用模型Faceg得到一个初始化模型 Faceinit，即对通用模型进行尺度变换，可用下式进行：
Mg和Minit分别是Faceg和Faceinit上的任一组对应点；sx是建 模对象和Faceg的两内眼角距离Seg_E长度的比值，sy是两眼中心到嘴中心距离Seg_EM长 度的比值；sz＝sy；再根据此对Faceg中所有点按上式变换；
姿态初始化是为了获得模型和图象间的对应，即计算从人脸generic坐标系到L摄像机坐 标系的刚体变换；即从下式求解Rg，1，Tg，1：
M1＝Rg，1Minit+Tg，1，M1和Minit分别是时刻t＝1时人脸和Faceinit间的任一组对应点；本 发明取上述两个模型在LE即左眼、RE即右眼和Mid_M即嘴中心这三个位置处的对应点， 按上述、θ、R、[r]x、T各式依次求得Rg，1和Tg，1；
再按下式从Rt，Tt，t＝1，...，N-1求出Rg，t，Tg，t，t＝2，...，N：
Rg，t＝Rt-1Rg，t-1，Tg，t＝Rt-1*Tg，t-1+Tt-1，t＝2，...，N；
(4).利用三维信息修正人脸模型：
(4.1)，手工拾取侧面特征点：
设：手工标出的侧面特征点为 kk＝1，...，NumSF，NumSF＝24为侧 面特征点的数目，side表示对应于侧面视图的时刻；
(4.2)，形状修正A：把Faceinit变形为Face1；我们先确定人脸三维结点列表P中一个子集 Subset的新位置，再通过一个已知的径向基函数算法来确定P中所有结点的坐标；
径向基函数算法如下所述：设要将Facestart变形为Facenew，已知一个Subset中的结点在 Facenew中的坐标为{New1，…，News}，记其在Facestart中的对应点为{Sta1，…，Stas}，S为Subset 中点的数目；现在对于模型中的任一位置pt，其在Facenew中的坐标可由下面的已知公 式求出：
是pt点在Facestart中的坐标，而{C1，…，Cs}可按下面的公式求出
φ是已知的径向基函数；
“形状修正A”在利用上面的径向基函数算法时，设置Facestart＝Faceinit， Facenew＝Face1，而Subset中的点包括两类：一类是人脸上的左眼内角LE、右眼内角RE、 左外嘴唇LM、右外嘴唇RM四个三维显著特征点，由二维显著特征点通过上述立体重建方 法得出；另一类是上述24个侧面特征点，它们的三维位置按下面的步骤从手工标出的二维点 中获得：
设：侧面特征点在人脸generic坐标系中的三维坐标为
则对于任一kk，有
其中 是姿态参数的分量表示，为已知量；可直接用线 性方程组的求解方法得到中的两个未知数和
(4.3)，形状修正B：按前述的径向基函数算法把Face1变形为Face2；
设Facestart＝Face1，Facenew＝Face2，Subset中除了包括“形状修正A”里已经确定了的 点之外，还有一些新添加的点，即从正面姿态下手工标出的人脸轮廓恢复出来的三维点，得 到这些三维点的基本方法是把当前人脸模型Face1向正面姿态下的图象平面作透视投影并计 算三维模型投影结果的二维轮廓点，最后依赖于这些轮廓点和上述图象上的人脸轮廓即闭合 多边形的匹配来计算其三维坐标：
首先，计算出Face1在t＝1时刻向手工选择的正面视图的图象平面上的投影轮廓 是模型三维结点列表P中点的标号， 而nNum1是Cont1中顶点的数目，它们即形状修正B中向Subset新添加的点；它们是通过判 断投影线与模型的交点来选择的，即在某一轮廓点上，投影线与模型不能有其他交点，其算 法如下：Face1为三维模型，透视投影中心为 已知量，则对模型中所有顶点 计算过Center点的投影线和模型的交点，据此判断该顶点是否在投影轮廓上，若Mnow为Face1 中任一点，Planejp为Face1中任一面片构成的平面，计算过Center点和Mnow的直线，即投影 线，与Planejp的交点为Cnow；如果对某一面片1≤jp≤nMesh而言，点Cnow在线段CMnow的内 部，则Mnow不在投影轮廓上；否则即为投影轮廓点；
再计算Cont1中的任一点在Face2的三维坐标：设上述图象上的人脸轮廓是 它是由二维点列表构成的多边形， nNumCI1是点的数目；则 在Face2的三维坐标
其中， 是pti点在Face1中的坐标；经点的投影方向为vpn为 已知量；vn为pti点在Faceg中的法向；
参数tline可根据二维直线 与闭合多边形Cont_img1的交点求 出：
设Cont_img1中的任一线段为seg＝(x0+sdx y0+sdy)T，0≤s≤1，s为线段参数，其他都 是已知量；先按下面的公式计算
对所有的seg求出这两个数值si和tline；对于很多seg上面的公式无解，即公式里的矩阵 不可逆；有解且得到的si满足0≤si≤1的seg有一个或两个，取与点 距离近的 那一个seg，它得到的tline即为所求；
(4.4)，形状修正C：按前述的径向基函数算法把Face2变形为Face3；
设Facestar，＝Face2，Facenew＝Face3，Subset中除了包括“形状修正A、B”里已经确定 了的点之外，还有一些新添加的点，即从手工标出的正面人脸特征轮廓点：眼睛、鼻孔和嘴 的轮廓点恢复出来的三维点：
设某特征轮廓点的二维坐标是 该点在Face2中的坐标为 设该点的Z坐标已经计算正确，即它在Face3中的Z坐标等于 则它在Face3中的坐标为 即这些特征轮 廓点在正面姿态下的投影等于实际图象上的特征点；
(4.5)，形状修正D：按前述的径向基函数算法把Face3变形为Face4；
设Facestart＝Face3，Facenew＝Face4，Subset中除了包括“形状修正A、B、C”里已经 确定了的点之外，还有一些新添加的点，即从手工标出的中间视图1即使Rg，t的旋转角最接 近25度的时刻t的人脸轮廓中恢复出来的三维点，具体步骤如下：设该视图对应于时刻int1， 与形状修正B所述的相同，先算出Face3在int1时的投影轮廓 是P中点的标号，nNum_int1是投影轮廓中的顶 点数目；对于任一点 它在Face3中的三维坐标 而 如前所述是L摄像机光心在generic坐标系中的坐标，则该点在Face4中的坐 标应满足 参数tline2同样根据直 线 与闭合多边形Cont_imgint1的交点求出；
(4.6)，形状修正E：按前述的径向基函数算法把Face4变形为Face5，即最终的形状模型 Faces；
设Facestart＝Face4，Facenew＝Face5，Subset中除了包括“形状修正A、B、C、D”里 已经确定了的点之外，还有一些新添加的点，即从手工标出的中间视图2即使Rg，t的旋转角 最接近40度的时刻t的人脸轮廓中恢复出来的三维点，具体的恢复步骤与形状修正D相同；
(5).纹理映射：
目前的纹理是通过L摄像机拍摄的正面和侧面视图，即IL，1和IL，side，生成，纹理映射是 要为创建的模型Faces生成看起来有相片一样真实感的柱面纹理图Textures，即要把采集的图 象转换到一个统一的柱面坐标系中并进行融合，它依次含有如下步骤：
(5.1)，生成柱面展开图：
首先对形状修正E的最终结果Faces中的任一三维点 按下式作柱 面映射：
映射的结果是一个二维平面，其上的点为{Cynpt}＝{(lgpt ltpt)}；这 里lgpt是柱面上的经度位置，用弧度值表示；ltpt则是柱面旋转轴向的坐标位置；
(5.2)，生成正面纹理图：
把所有点在已知的正面姿态参数Rg，1，Tg，1下按下式在正面视图IL，1上投影， 是Faces中三维点在IL，1上的二维投影：
其中 是旋转矩阵和平移向量的分量 表示；
接着，把IL，1中的各像素映射到Cyn平面上，即对Faces中的任一面片 facet＝(pt1 pt2 pt3)T，mproj＝Δ(p1，p2，p3)表示这三点在IL，1平面上构成的三角形， mcyn＝Δ(p1′，p2′，p3′)是Cyn平面上的对应三角形，即 iT1＝1，2，3；再对三角形mcyn内的点[lg lt]按下式计 算[x y]：
a1～a6由下式求出：
再令正面纹理图
(5.3)，生成侧面纹理图：
把所有点在已知的侧面姿态参数Rg，side，Tg，side下按下式在侧面视图IL，side上投影， 是Faces中三维点在IL，side上的二维投影：
再对Faces中的任一面片facet＝(pt1 pt2 pt3)T，设msproj＝Δ(ps1，ps2，ps3)， iT2＝1，2，3；再对三角形mcyn内的点[lg lt]按下式计算[xs ys]：
as1～as6由下式求出：
再令侧面纹理图
(5.4)，反射侧面纹理图以得到另一侧的纹理图，这是因为Faces的拓扑结构是完全对称的；
设：直接获得纹理的那一侧脸上的任一三角面片mcyn＝Δ(p1，p2，p3)，它在另一侧脸上的 对称面片是mcyn′＝Δ(p1′，p2′，p3′)，其中 iT3＝1，2，3；
对mcyn中的任一点p＝(lg lt)T按下式计算[lg′ lt′]，
rs1～rs6按下式求出：
再从直接有纹理的一侧反射生成另一侧：
(5.5)，正面的和侧面的纹理图的融合，得到最终纹理图Textures；
设定阈值lgmin和lgmax，Textures在任一位置Cyn＝(lg lt)T处的取值按下式给出：
2.根据权利要求1所述的融合多视角、多线索二维信息的人脸三维模型的建立方法，其 特征在于：所述的姿态估计步骤中的更新可靠度量序列{bReliablek}的方法，它依次含有如下 步骤：
(1).根据对序列Errt的统计分析计算正确对应的数目K_True；
(1.1)，设：K_True的初始值为Kr(t)*ratio，ratio为预先设定的常数； Errt＝{εs＝‖Mt+1(ks)-RtMt(ks)-Tt‖2|s＝1，...，Kr(t)}中按从小到大排序的前K_True个元素组 成的子序列为E；
(1.2)，计算E的均值 和标准差
(1.3)，在Errt中寻找满足‖εs-ε‖>3*σ的点，设数目为InValid_T；
(1.4)，若0≤Kr(t)-K_True-InValid_T≤1，即K_True个对应正确，终止；否则，转步 骤(1.5)；
(1.5)，若K-K_True-InValid_T>1，令K_True加1，转步骤(1.2)；否则，转步骤(1.6)；
(1.6)，若K-K_True-InValid_T<0，令K_True减1，转步骤(1.2)；
(2).设Errt中前K_True个较小的εs对应的标号为{su}，u＝1，...，K_True，把{b Reliableks} 中对应于这些标号的位置赋值为Yes，其他赋值为No。
3.根据权利要求1所述的融合多视角、多线索二维信息的人脸三维模型的建立方法，其 特征在于：所述的姿态估计步骤中的设立是否要重新调用特征检测和匹配模块的布尔标志量 bReDetect、是否要进行姿态求精的布尔标志量bRefine这两者的方法如下：
(1).若t＝N，即当前图象是序列中的最后一帧，则bReDetect＝No，bRefine＝Yes，结束； 否则，转步骤(2)；
(2).若当前Rt的旋转角θt超过预设的阈值θ_bound，则bReDetect＝Yes，bRefine＝Yes， 结束；否则，转步骤(3)；
(3).若可靠三维点的数目Kr(t)(即{bReliablek}，k＝1，...，K(t)中取值为Yes的点数)小于 预设的阈值K_bound，则bReDetect＝Yes，bRefine＝Yes；否则，转步骤(4)；
(4).bReDetect＝No；如果t—t0是预设的常数Inc_refine的倍数，则bRefine＝Yes，即求精 是每隔Inc_refine就做一次；否则bRefine＝No；
(5).令t＝t+1，进入下一时刻的处理。
4.根据权利要求3所述的融合多视角、多线索二维信息的人脸三维模型的建立方法，其特 征在于：所述的姿态估计步骤中，bRefine为Yes表示在当前时刻作姿态求精，即得出更精确 的姿态参数：
设选择在连续Num＝t-tStart+1个时刻始终正确的三维对应来求精所有相关的姿态参 数{RtStart+τ，TtStart+τ}，τ＝0，...，t-tStart-1，tStart为求精运算的起始时刻，现取t时刻之前的第二 个设置bRefine为Yes的时刻，若tStart
其中RatStart+τ，TatStart+τ是从tStart时刻到tStart+τ时刻的刚体变换系数； {MtStart(kpp)}，pp＝1，..，Kr(t)，是在这Num个时刻中都可靠的三维点数组，共Kr(t)个，pp为 三维点的标号；m为L或R，表示摄像机；是在m号摄像机、tStart+τ时刻图象 中的二维检测或跟踪的结果，即
Am为AL或AR；是pp号三维点在tStart+τ时刻， 投影到m摄像机上的结果，即
其中 和 是从 L摄像机坐标系到m摄像机坐标系的刚体变换，即 另外，
优化f_obj就能得到姿态参数的精确值；先用Levenberg-Marquardt算法求解 {RatStart+τ，TatStart+τ，τ＝1，...，t-tStart；MtStart(kpp)，pp＝1，...，Kr(t)}；再用下式把 {RatStart+τ，TatStart+τ}，τ＝1，...，t-tStart转换到{RtStart+τ，TtStart+τ}，τ＝0，...，t-tStart-1：
RtStart＝RatStart+1，TtStart＝TatStart+1，
技术领域
融合多视角、多线索二维信息的人脸三维模型的建立方法属于计算机人脸图像重建技术 领域。
背景技术
在“利用桌面摄像机克隆人脸”(《第八届计算机视觉国际会议论文集》，2001，第二 卷，pp.745-745)一文中，作者张正友、Z.Liu等人提出了一种利用单目视频中的二维对应创 建人脸模型的方法。该方法的原理是：手工选择两幅接近正面视图的图象，通过摄像机自标 定和简化的人脸形状模型，基于二维对应用非线性求解的办法得到特定人的模型。该方法虽 然利用了多视角信息，但只是使用了多幅图象上对应的二维点信息而没有充分利用其他的形 状约束，得到的模型不够准确，尤其是在纹理不丰富的区域更是如此。
这些方法有一些共同的缺点：(1)模型不够精细；(2)不能做MPEG4(运动图像压缩 编码组织标准4)兼容的动画。
发明内容
本发明的目的在于提供一种融合多视角、多线索二维信息的人脸三维模型的建立方法。
本发明涉及一种可动画三维人脸的建模方法，首先用标定好的立体摄像机拍摄建模对象 的人脸从正面逐步转到侧面的视频序列(拍摄过程人脸应无表情变化)；然后是一个全自动 的姿态估计算法以得到序列中各时刻的精确的姿态参数；接着在模型和姿态初始化之后，根 据两类二维线索通过不同的处理方案抽取三维信息：一类是基于图象间对应点进行立体重建， 另一类是模型投影轮廓与图象轮廓的匹配、修正；最后，融合在多视角中提取出的、对应于 不同二维线索的三维信息建立人脸模型。本方法与现有方法比，得到的人脸模型有更细致、 更精确的形状，并能做MPEG4兼容的人脸动画。
本发明的特征在于：
1.它首先用与人脸上、下位置相对应的两台相同且外极线方向基本沿垂直方向的摄像机 在人脸无表情变化的情况下拍摄建模对象的人脸从正面逐步转约90度到侧面的视频序列，再 用两块相同的采集卡来采集上下两路同步视频并输入计算机中；然后，从人脸的通用模型 Faceg开始，根据获得的建模对象的二维信息逐步进行修改，最终得到特定人的人脸模型Faces， 其中修改的只是人脸上的三维结点列表P中的三维点坐标，而人脸上的三角形面片列表F保 持不变，即只改变模型形状P而模型的拓扑结构F不变；具体而言，它依次含有如下步骤：
(1).输入立体视频序列IL，t和IR，t：
L，R对应于两台摄像机，t＝1，...，N，表示N帧不同时刻的图像：I则是由像素组成的二维 数组，每个像素包括R、G、B三个分量；t＝1的时刻对应于人脸正面姿态，通过人的判断手 工选择；
(2).姿态估计，即求解人头旋转过程中的姿态参数(即刚体变换系数)：
用下述刚体变换公式表示上述图象序列中人脸的帧间运动：
Mt+1＝RtMt+Tt，t＝1，...，N-1，
Mt、Mt+1分别表示第t帧、第t+1帧时人脸上任一点在L摄像机坐标系中的坐标；Mt、 Mt+1是人脸上的同一位置在不同时刻的三维坐标，即是帧间三维对应点；用可靠三维点的运 动就可估计姿态参数；
{Rt，Tt，t＝1，...，N-1}是第t帧到第t+1帧姿态参数的旋转和平移分量；
姿态估计依次含有如下步骤：
(2.1)，起始时刻：先进行二维特征的检测和匹配，再通过立体重建得到起始时刻的三维 点：
首先，设置特征检测帧t0＝t，并在IL，t图上作二维特征检测，其步骤如下：
用已有的KLT算法计算IL，t图上特征点p(x＝i，y＝j)处的图像变化程度G(i，j)：
为待测点p处图像的梯度向量：
I(i+1，j)为该点处图像的灰度，余类推；
以待测点p为中心的窗口W的尺寸为(2*blockx+1)*(2*blocky+1)；然后再用已有的matlab 工具包求出G的最小奇异值sing_min，再选择sing_min很大的点p为检测出的特征点；把在 t时刻检测出的特征点记为{pL，t(k)}，L，t表示图象序号，K(t)是特征点总数，k为其标号， 1≤k≤K(t)；
其次，对同一时刻两个摄像机L，R拍摄同一场景得到的两幅图像即立体图对作二维特征 匹配；同一个三维点在不同图像中的投影都为二维点，它们之间互称为二维对应，立体图对 间的二维对应称为匹配：
对于IL，t图中某一点pL，t(k)＝(xL，t(k)yL，t(k))T，它的匹配点pR，t(k)＝(xR，t(k)yR，t(k))T应位 于一条由pL，t(k)决定的直线上，在这个一维空间中，选择使一个局部窗口W内总灰度差异最 小的点作为pL，t(k)的匹配点pR，t(k)，总灰度差异的表达式：
(i，j)是直线上的任一点；记时刻t的匹配结果为{pR，t(k)}，总数目为K(t)个， 每点都与pL，t(k)相对应，在下述公式表示的直线上搜索，使Diff(i，j)最小的位置：
其中：
fL，fR，u0L，v0L，u0R，v0R为 摄像机内参数，已知量；Rc，Tc＝[Tc，xTc，yTc，z]T分别为旋转和平移分量，属摄像机外参数， 已知量；
再次，在获得同一时刻立体图对间的匹配结果pL，t和pR，t后，通过立体重建得到用摄像机 坐标系表示的三维点Mm，t，m为摄像机L或R，Mm，t为t时刻三维点在m摄像机坐标系中的 坐标；
根据两摄像机的配置关系式：
ML＝RcMR+Tc，ML、MR分别是同一点在这两个摄像机坐标系中的坐标；以及摄像机 透视投影公式：
sm(即sL或sR)是待定的比例因子(由稍后的公式求出)，
得到三维点在摄像机L、R中的坐标ML，t和MR，t：
而
在姿态估计的余下部分，三维点都用L摄像机坐标系中的坐标来表示，因此省略下标L， 即将ML，t简记为Mt；
(2.2)，后继时刻：对t+1时刻的立体图对用KLT算法进行二维特征跟踪，再同样通过立 体重建得到各后继时刻的三维点；跟踪是指同一摄像机拍摄的不同时刻图像间的二维对应：
在二维特征跟踪时，设检测的特征点在图像Im，t(p)(m为L或R)上，p(x＝i，y＝j)是任一 特征点，则其在图像Im，t+1上的跟踪结果为其中 G已如上述，而
记跟踪得到的特征点为{pm，t+1(k)}，m表示L或R，t+1为时刻：
对于摄像机L为：
对于摄像机R为：
根据任一对匹配的跟踪点计算三维点{Mt+1(k)|k＝1，...，K(t)}(如前所述，这里省略了下标 L)；
其中
(2.3)，姿态参数初始化：根据两个时刻t和t+1之间的一组三维对应点{Mt(k)|k＝1，...，K(t)} 和{Mt+1(k)|k＝1，...，K(t)}，和表示可靠性的序列{bReliablek}，k＝1，...，K(t)从式 Mt+1(k)＝RtMt(k)+Tt，t＝1，...，N-1中用基于最小中值的鲁棒估计算法求解Rt，Tt：
设：定义并不断更新布尔数组{bReliablek}为可靠度量序列，Yes表示可靠，No表示不可 靠，其中取值为Yes的项的标号为{ks}，s＝1，...，Kr(t)，Kr(t)是取值为Yes的项的总数，从其 中挑出Num_subset个子集，每个子集含有三对三维对应点，即 {Setn＝{N1，n，N2，n，N3，n}}，n＝1，...，Num_subset，1≤N1，n，N2，n，N3，n≤Kr(t)，每个三元子集中三维对 应点为： 和 可简记为 {(Mim，Mim’)|im＝0，1，2}，它表示三对不共线的三维对应点，Rt，Tt满足 Mim′＝RtMim+Tt|im＝0，1，2；然后从下列公式中求出R和T(略去下标t)：
其中
是旋转的中心轴，称旋转轴，3*1的向量；
θ是绕的旋转角；
[r]x是根据3*1向量r＝[rxryrz]T定义的反对称矩阵： 而 是单位矩阵；
平移分量T由下式估计出：
对每个子集Setn按上述公式估计出姿态参数，记为Rt，n，Tt，n；再对所有的可靠三维对应 {Mt(ks)，Mt+1(ks)}，s＝1，...，Kr(t)按下式计算对应点Mt(ks)，Mt+1(ks)与姿态参数Rt，n，Tt，n的符合 程度εn，s，该值越小，则两者越符合：
εn，s＝‖Mt+1(ks)-Rt，nMt(ks)-Tt，n‖2；
再从集合{Rt，n，Tt，n}，n＝1，...，Num_subset中选取Rt，Tt：
{Rt，Tt}＝{Rn_m，Tn_m}，
表示针对下标s取集合{εn，s|s＝1，...，Kr(t)}的中值，它对不同的n构成一个新的数组 而则是在这个中值数组中选择取最小值的那个n；
(2.4)，在t＝t+1时刻，设置：
是否需要重新调用特征检测和匹配模块的布尔标志量bReDetect；
是否需要进行姿态求精的布尔标志量bRefine
(2.5)，当t＝N时，姿态估计结束；
(3).模型初始化和姿态初始化：
模型初始化是通用模型特定化方法的初始步骤，它从通用模型Faceg得到一个初始化模型 Faceinit，即对通用模型进行尺度变换，可用下式进行：
Mg和Minit分别是Faceg和Faceinit上的任一组对应点；sx是建 模对象和Faceg的两内眼角距离Seg_E长度的比值，sy是两眼中心到嘴中心距离Seg_EM长 度的比值；sz＝sy；再根据此对Faceg中所有点按上式变换；
姿态初始化是为了获得模型和图象间的对应，即计算从人脸generic坐标系到L摄像机坐 标系的刚体变换；即从下式求解Rg，1，Tg，1：
M1＝Rg，1Minit+Tg，1，M1和Minit分别是时刻t＝1时人脸和Faceinit间的任一组对应点；本 发明取上述两个模型在LE(左眼)、RE(右眼)和Mid_M(嘴中心)这三个位置处的对应 点，按上述、θ、R、[r]x、T各式依次求得Rg，1和Tg，1；
再按下式从Rt，Tt，t＝1，...，N-1求出Rg，t，Tg，t，t＝2，...，N：
Rg，t＝Rt-1Rg，t-1，Tg，t＝Rt-1*Tg，t-1+Tt-1，t＝2，...，N；
(4).利用三维信息修正人脸模型：
(4.1)，手工拾取侧面特征点：
设：手工标出的侧面特征点为 kk=1,...,NumSF,NumSF＝24为侧 面特征点的数目，side表示对应于侧面视图的时刻；
(4.2)，形状修正A：把Faceinit变形为Face1；我们先确定人脸三维结点列表P中一个子集 Subset的新位置，再通过一个已知的径向基函数算法来确定P中所有结点的坐标；
径向基函数算法如下所述：设要将Facestart变形为Facenew，已知一个Subset中的结点在 Facenew中的坐标为{New1，…，News}，记其在Facestart中的对应点为{Sta1，…，Stas}，S为Subset 中点的数目；现在对于模型中的任一位置pt，其在Facenew中的坐标可由下面的已知公 式求出：
是pt点在Facestart中的坐标，而{C1，…，CS}可按下面的公式求出
φ是已知的径向基函数；
“形状修正A”在利用上面的径向基函数算法时，设置Facestart＝Faceinit， Facenew＝Face1。而Subset中的点包括两类：一类是人脸上的左眼内角LE、右眼内角RE、 左外嘴唇LM、右外嘴唇RM四个三维显著特征点，由二维显著特征点通过上述立体重建方 法得出；另一类是上述24个侧面特征点，它们的三维位置按下面的步骤从手工标出的二维点 中获得：
设：侧面特征点在人脸generic坐标系中的三维坐标为
则对于任一kk，有
其中 是姿态参数的分量表示，为已知量；可直接用线 性方程组的求解方法得到中的两个未知数和
(4.3)，形状修正B：按前述的径向基函数算法把Face1变形为Face2；
设Facestart＝Face1，Facenew＝Face2，Subset中除了包括“形状修正A”里已经确定了的 点之外，还有一些新添加的点，即从正面姿态下手工标出的人脸轮廓恢复出来的三维点，得 到这些三维点的基本方法是把当前人脸模型Face1向正面姿态下的图象平面作透视投影并计 算三维模型投影结果的二维轮廓点，最后依赖于这些轮廓点和上述图象上的人脸轮廓即闭合 多边形的匹配来计算其三维坐标：
首先，计算出Face1在t＝1时刻向手工选择的正面视图的图象平面上的投影轮廓 是模型三维结点列表P中点的标号， 而nNum1是Cont1中顶点的数目，它们即形状修正B中向Subset新添加的点；它们是通过判 断投影线与模型的交点来选择的，即在某一轮廓点上，投影线与模型不能有其他交点，其算 法如下：Face1为三维模型，透视投影中心为 已知量，则对模型中所有顶点 计算过Center点的投影线和模型的交点，据此判断该顶点是否在投影轮廓上，若Mnow为Face1 中任一点，Planejp为Face1中任一面片构成的平面，计算过Center点和Mnow的直线，即投影 线，与Planejp的交点为Cnow；如果对某一面片1≤jp≤nMesh而言，点Cnow在线段CMnow的内 部，则Mnow不在投影轮廓上；否则即为投影轮廓点；
再计算Cont1中的任一点在Face2的三维坐标：设上述图象上的人脸轮廓是 它是由二维点列表构成的多边形， nNumCI1是点的数目；则 在Face2的三维坐标
其中， 是pti点在Face1中的坐标；经点的投影方向为vpn为 已知量；vn为pti点在Faceg中的法向；
参数tline可根据二维直线 与闭合多边形Cont_img1的交点求 出：
设Cont_img1中的任一线段为seg＝(x0+sdx y0+sdy)T，0≤s≤1，s为线段参数，其他都 是已知量；先按下面的公式计算
对所有的seg求出这两个数值si和tline；对于很多seg上面的公式无解，即公式里的矩阵 不可逆；有解且得到的si满足0≤si≤1的seg有一个或两个，取与点 距离近的 那一个seg，它得到的tline即为所求；
(4.4)，形状修正C：按前述的径向基函数算法把Face2变形为Face3；
设Facestart＝Face2，Facenew＝Face3，Subset中除了包括“形状修正A、B”里已经确定 了的点之外，还有一些新添加的点，即从手工标出的正面人脸特征轮廓点：眼睛、鼻孔和嘴 的轮廓点恢复出来的三维点：
设某特征轮廓点的二维坐标是 该点在Face2中的坐标为 设该点的Z坐标已经计算正确，即它在Face3中的Z坐标等于 则它在Face3中的坐标为 即这些特征轮 廓点在正面姿态下的投影等于实际图象上的特征点；
(4.5)，形状修正D：按前述的径向基函数算法把Face3变形为Face4；
设Facestart＝Face3，Facenew＝Face4，Subset中除了包括“形状修正A、B、C”里已经 确定了的点之外，还有一些新添加的点，即从手工标出的中间视图1即使Rg，t的旋转角最接 近25度的时刻t的人脸轮廓中恢复出来的三维点，具体步骤如下：设该视图对应于时刻int1， 与形状修正B所述的相同，先算出Face3在int1时的投影轮廓 是P中点的标号，nNum_int1是投影轮廓中的顶 点数目；对于任一点 它在Face3中的三维坐标 而 如前所述是L摄像机光心在generic坐标系中的坐标，则该点在Face4中的坐 标应满足 参数tline2同样根据直 线 与闭合多边形Cont_imgint1的交点求出；
(4.6)，形状修正E：按前述的径向基函数算法把Face4变形为Face5(即最终的形状模型 Faces)；
设Facestart＝Face4，Facenew＝Face5，Subset中除了包括“形状修正A、B、C、D”里 已经确定了的点之外，还有一些新添加的点，即从手工标出的中间视图2即使Rg，t的旋转角 最接近40度的时刻t的人脸轮廓中恢复出来的三维点，具体的恢复步骤与形状修正D相同；
(5).纹理映射：
目前的纹理是通过L摄像机拍摄的正面和侧面视图(IL，1和IL，side)生成，纹理映射是要 为创建的模型Faces生成看起来有相片一样真实感的柱面纹理图Textures，即要把采集的图象 转换到一个统一的柱面坐标系中并进行融合，它依次含有如下步骤：
(5.1)，生成柱面展开图：
首先对形状修正E的最终结果Faces中的任一三维点 按下式作柱 面映射：
映射的结果是一个二维平面，其上的点为{Cynpt}＝{(lgpt ltpt)}；这 里lgpt是柱面上的经度位置，用弧度值表示；ltpt则是柱面旋转轴向的坐标位置；
(5.2)，生成正面纹理图：
把所有点在已知的正面姿态参数Rg，1，Tg，1下按下式在正面视图IL，1上投影， 是Faces中三维点在IL，1上的二维投影：
其中 是旋转矩阵和平移向量的分量 表示；
接着，把IL，1中的各像素映射到Cyn平面上，即对Faces中的任一面片 facet＝(pt1pt2pt3)T，mproj＝Δ(p1，p2，p3)表示这三点在IL，1平面上构成的三角形， mcyn＝Δ(p1′，p2′，p3′)是Cyn平面上的对应三角形，即 iT1＝1，2，3；再对三角形mcyn内的点[lg lt]按下式计 算[x y]：
a1～a6由下式求出：
再令正面纹理图
(5.3)，生成侧面纹理图：
把所有点在已知的侧面姿态参数Rg，side，Tg，side下按下式在侧面视图IL，side上投影， 是Faces中三维点在IL，side上的二维投影：
再对Faces中的任一面片facet＝(pt1 pt2 pt3)T，设msproj＝Δ(ps1，ps2，ps3)， iT2＝1，2，3；再对三角形mcyn内的点[lg lt]按下式计算[xs ys]：
as1～as6由下式求出：
再令侧面纹理图
(5.4)，反射侧面纹理图以得到另一侧的纹理图，这是因为Faces的拓扑结构是完全对称的；
设：直接获得纹理的那一侧脸上的任一三角面片mcyn＝Δ(p1，p2，p3)，它在另一侧脸上的 对称面片是mcyn′＝Δ(p1′，p2′，p3′)，其中 iT3＝1，2，3；
对mcyn中的任一点p＝(lg lt)T按下式计算[lg′lt′]，
rs1～rs6按下式求出：
再从直接有纹理的一侧反射生成另一侧：
(5.5)，正面的和侧面的纹理图的融合，得到最终纹理图Textures；
设定阈值lgmin和lgmax，Textures在任一位置Cyn＝(lglt)T处的取值按下式给出：
2.所述的姿态估计步骤中的更新可靠度量序列{bReliablek}的方法，它依次含有如下步骤：
(1).根据对序列Errt的统计分析计算正确对应的数目K_True；
(1.1)，设：K_True的初始值为Kr(t)*ratio，ratio为预先设定的常数； Errt＝{εs＝‖Mt+1(ks)-RtMt(ks)-Tt‖2|s＝1，...，Kr(t)}中按从小到大排序的前K_True个元素组 成的子序列为E；
(1.2)，计算E的均值 和标准差
(1.3)，在Errt中寻找满足‖εs-ε‖>3*σ的点，设数目为InValid_T；
(1.4)，若0≤Kr(t)-K_True-InValid_T≤1，即K_True个对应正确，终止；否则，转步 骤(1.5)；
(1.5)，若K-K_True-InValid_T>1，令K_True加1，转步骤(1.2)；否则，转步骤(1.6)；
(1.6)，若K-K_True-InValid_T<0，令K_True减1，转步骤(1.2)；
(2).设Errt中前K_True个较小的εs对应的标号为{su}，u＝1，...，K_True，把 中对应于这些标号的位置赋值为Yes，其他赋值为No；
3.所述的姿态估计步骤中的设立是否要重新调用特征检测和匹配模块的布尔标志量 bReDetect、是否要进行姿态求精的布尔标志量bRefine这两者的方法如下：
(1).若t＝N，即当前图象是序列中的最后一帧，则bReDetect＝No，bRefine＝Yes，结束； 否则，转步骤(2)；
(2).若当前Rt的旋转角θt超过预设的阈值θ_bound，则bReDetect＝Yes，bRefine＝Yes， 结束；否则，转步骤(3)；
(3).若可靠三维点的数目Kr(t)(即{b Reliablek}，k＝1，...，K(t)中取值为Yes的点数)小于 预设的阈值K_bound，则bReDetect＝Yes，bRefine＝Yes；否则，转步骤(4)；
(4).bReDetect＝No；如果t—t0是预设的常数Inc_refine的倍数，则bRefine＝Yes，即求精 是每隔Inc_refine就做一次；否则bRefine＝No；
(5).令t＝t+1，进入下一时刻的处理；
4.所述的姿态估计步骤中，bRefine为Yes表示在当前时刻作姿态求精，即得出更精确的 姿态参数：
设选择在连续Num＝t-tStart+1个时刻始终正确的三维对应来求精所有相关的姿态参 数{RtStart+τ，TtStart+τ}，τ＝0，...，t-tStart-1，tStart为求精运算的起始时刻，现取t时刻之前的第二 个设置bRefine为Yes的时刻，若tStart
其中RatStart+τ，TatStart+τ是从tStart时刻到tStart+τ时刻的刚体变换系数； {MtStart(kpp)}，pp＝1，...，Kr(t)，是在这Num个时刻中都可靠的三维点数组，共Kr(t)个，pp为 三维点的标号；m为L或R，表示摄像机；是在m号摄像机、tStart+τ时刻图象 中的二维检测或跟踪的结果，即
Am为AL或AR；是pp号三维点在tStart+τ时刻， 投影到m摄像机上的结果，即
其中 和 是从 L摄像机坐标系到m摄像机坐标系的刚体变换，即 另外，
优化f_obj就能得到姿态参数的精确值；先用Levenberg-Marquardt算法求解 {RatStart+τ，TatStart+τ，τ＝1，...，t-tStart；MtStart(kpp)，pp＝1，...，Kr(t)}；再用下式把 {RatStart+τ，TatStart+τ}，τ＝1，...，t-tStart转换到{RtStart+τ，TtStart+τ}，τ＝0，...，t-tStart-1：
RtStart＝RatStart+1，TtStart＝TatStart+1，
试验证明：它达到了预期目的。
附图说明
图1.人脸建模流程图
图2.姿态估计流程图
图3.纹理映射流程图
图4.摄像机坐标系与像平面上二维坐标系的关系
图5.正面人脸的显著特征点及线段
图6.典型的人脸侧面投影以及一些侧面特征点
图7.叠加显示的正面人脸轮廓
图8.叠加显示的正面人脸特征轮廓点
图9.叠加显示的中间视图1人脸轮廓(25度角)
图10.叠加显示的中间视图2人脸轮廓(40度角)
图11.上下放置的立体摄像机配置
图12.系统采集的典型立体图对(上下两行分别对应于top，bottom摄像机，每一列是同 一时刻的图象，标出的数字是在整个序列中的帧序号)
图13.通用模型的中性和动画状态(从左到右依次是中性线段图、中性面片图、动画线 段图和动画面片图)
图14.不同姿态下特征检测和匹配的实验结果(第一行是检测结果，第二行是对应的立 体图对间的匹配结果；标出的数字是在整个序列中的帧序号)
图15.正面姿态的立体图对上标出的二维显著特征点(特征点叠加显示在原始图象上)
图16.侧面特征点的叠加显示
图17.实际得到的Texture1
图18.第三十五步直接得到的Textureside
图19.最终得到的Textureside
图20.柱面纹理图Textures
图21.映射了Textures的建模结果Faces(从不同视角查看)
具体实施方式
下面我们将详细说明建模步骤，并给出一个具体实施实例。
人脸建模的总体流程如图1所示，各个步骤的前后顺序以及输入和输出都明确标出。在 “处理模块”这部分，框图中的数字编号是本节中详细对该模块进行阐述部分的小节号。此 外，“输入”指那些建模之前通过手工干预准备好的数据(视频数据、Faceg，以及需要手工 标出的二维点)，“中间结果”则是系统在建模流程中自动获得的数据。
输入的立体视频序列记为IL，t和IR，t，t＝1，...，N表示不同的时刻，N则是视频序列的总帧 数，L，R对应于两台摄像机，I则是由像素组成的二维数组，每个像素包括R、G、B三个分 量。t＝1的时刻对应于人脸正面姿态，通过人的判断手工选择。
下面我们详细描述各步骤的具体处理方法。
1.1姿态估计
这一步要求解人头旋转过程中的姿态参数(也称为刚体变换系数)。记为 {Rt，Tt，t＝1，...，N-1}，其中每个Rt，Tt是第t帧到第t+1帧的人脸刚体变换系数，即对于t时刻 和t+1时刻的任一对三维对应点Mt和Mt+1(所谓对应，是指它们是人脸上的同一位置在不同 时刻的三维坐标)，有下面的公式
Mt+1＝RtMt+Tt，t＝1，...，N-1
其中三维点所在的坐标系是L摄像机坐标系。摄像机坐标系的定义如图4所示，以透视 摄像机的光心C为坐标原点，Z轴沿光轴方向，X，Y轴分别和像平面上二维坐标系的x，y轴 平行。
姿态估计的整个处理流程如图2所示，下面我们先对算法的原理作简单说明，在详细介 绍各步骤。
我们的思路是利用可靠三维点的运动来估计参数，而三维点需要根据投影到不同图象上 的二维对应点(在下面被称为“二维特征”。关于“二维对应”，是指同一个三维点在不同 图象中都投影为二维点，它们之间互称为二维对应点；其中立体图对间的二维对应称为“匹 配”，而同一摄像机拍摄的不同时刻图象间的二维对应则称为“跟踪”。这里“立体图对” 指同一时刻L和R摄像机拍摄的一对图像)计算出来。具体地，我们首先要得到三维点，这 在起始时刻是通过二维特征检测和匹配，然后通过立体重建得到；而对于后续时刻则是在L 和R两个图象序列中分别进行跟踪，接着同样通过立体重建得到。有了不同时刻的三维点后 我们就可以估计姿态。从未知量的个数来考虑，只需要三对两个时刻的对应三维点就可完成， 但要求重建出的三维点是可靠的，而实际实验中很难保证所有的三维点都可靠。为了评价三 维点的可靠性，我们定义并不断更新数组{b Reliablek}，这是一个布尔数组，取值为Yes(表 示该点可靠)或No(表示不可靠)。该数组会随着跟踪的进行而被更新，刚开始时初始化所 有点都是可靠的，后来有些点则会被判断为不可靠的，我们作出这种判断的方法将在下面详 细说明。
直接根据三维点在相邻两时刻的运动计算出的姿态参数误差较大，我们会在某些时刻根 据一段较长时间内的所有可靠三维点进行姿态求精以减小误差。何时进行这种求精运算是由 布尔标志bRefine来控制的，它取值Yes表示当前时刻进行姿态求精，No则相反。
随着跟踪的进行，可靠的三维点会越来越少。为了使姿态估计可以继续进行并得到精确 的结果，我们有必要重新进行二维特征检测和匹配以及立体重建，找出新的可靠三维点。何 时进行这种操作由布尔标志bReDetect控制，它取值Yes表示当前时刻需要进行找出新三维 点的工作，No则相反。
bRefine和bReDetect这两个布尔标志的设置方法也将在下面详细介绍。
下面我们详细阐述完整的过程：
1.1.1二维特征的检测
我们利用已有的KLT算法在灰度图上完成二维特征的检测。其原理是：设t时刻的m摄 像机(m为L或R，表示两个摄像机中的某一个)得到的图象为Im，t(p)，其中p＝(x，y)表示 图象上的任一点；该点处图象的梯度向量为 计算2*2的 矩阵 其中W是以p为中心的一个矩形窗口。我们求出G的最小奇异值sing_min (利用已有的matlab工具包)，并选择那些sing_min很大的p点作为检测出的特征点。从物 理意义上看，G实际是对图象上局部灰度变化丰富程度的表示：如果局部的灰度恒定，这时 梯度向量g为0，G为0，所有的奇异值都是0；如果局部的灰度变化有固定的方向(例如图 象上比较强的单一边缘处)，梯度向量g基本都垂直于该方向，这时G最小的奇异值接近于 0；只有在沿各个方向的灰度变化都比较显著的局部(例如两条边的交点)才满足G的最小奇 异值很大。因此，我们实际是选择那些局部灰度变化丰富的点作为特征点，这样的点在跟踪 和匹配时可靠程度较高。
计算机实现时我们需要对上面的公式作离散化处理，具体的计算公式如下：
这里blockx，blocky决定了窗口W包含的网格点数目。实际上，该窗口以待测点(i，j)为 中心，网格点总数为(2*blockx+1)*(2*blocky+1)。
在图2中，t时刻检测出的特征点记为{pL，t(k)}，L，t表示图象序号。特征点的总数目为 K(t)，k满足1≤k≤K(t)，表示点的标号。
这时我们令t0＝t，表示将要被跟踪的这些特征在哪一时刻被检测出来。
1.1.2二维特征匹配
匹配问题根据已知方法完成。简单介绍其原理如下：
设两摄像机间的配置关系用下式表示
ML＝RcMR+Tc (1.1.3)
L和R表示两个摄像机坐标系，ML和MR是同一三维点在这两个坐标系中的坐标表示； Rc和Tc＝[Tc，xTc，yTc，z]T分别是旋转和平移分量(合称立体摄像机的外参数，是已知量)， 前者是一3*3的矩阵，后者是三维向量。
有了(1.1.3)后，令 (该矩阵 被称为基础矩阵)，其中fL，fR，u0L，v0L，u0R，v0R是摄像机的内参数，都是已知量；其矩阵表示 是
则对于L，t图中的某一点pL，t(k)＝(xL，t(k)yL，t(k))T，它在R，t图中的匹配点 pR，t(k)＝(xR，t(k)yR，t(k))T应满足下面的已知公式(1.1.4)：
(1.1.4)式的物理意义是pR，t(k)位于由pL，t(k)决定的一条直线上，这样为寻找该点的搜索空 间从整个图象降为该直线，即从二维降到了一维。在这个一维空间里，选择使一个局部窗口 Wm内总灰度差异最小的点作为匹配结果，即要使下面的(1.1.5)式最小：
Diff＝∫[IR，t(pR，t(k))-IL，t(pL，t(k))]2dp (1.1.5)
离散化后，实际计算时的公式是：
其中(i，j)是直线上的任一点。
图2中t时刻的匹配结果为{pR，t(k)}，总数目也为K(t)，每个点都和pL，t(k)相对应。对每 一pL，t(k)的匹配过程是独立的，即在公式(1.1.4)表示的直线上搜索使(1.1.6)式最小的位置。
获得匹配后用立体重建的方法获得三维点序列{Mt(k)|k＝1，...，K(t)}。立体重建的过程对每 一对匹配点是独立的，原理如下：
首先，我们有摄像机的透视投影原理，即下面的已知公式， 其中m为L或R，表示两个摄像机中的任一个；Mm，t为t时刻三维 点在m摄像机坐标系中的坐标；sm是待求的一个比例因子。我们省略了(k)来表示任一点。
再结合两摄像机间的配置关系(1.1.3)可求解出三维点，ML，t的表达式为
其中 而
而MR，t则可按下面的公式求出
如果不特别说明，在本文中提到的立体重建出的三维点都是指在L摄像机坐标系中的 ML，t。为方便起见，在陈述中我们去掉L，简写为Mt。对所有的特征点重复上面的立体重建 步骤，我们得到三维点序列{Mt(k)|k＝1，...，K(t)}。
从物理意义上看，每一二维投影点实际是确定了三维空间的一条直线，这样根据同一三 维点在两个摄像机中的一对二维点，通过求解两条三维直线的交点，就可以确定该三维点。
同时，我们还定义布尔变量序列{bReliablek}，k＝1，...，K(t)，它表示各三维点是否可靠，取 值为Yes或No，目前序列中的所有值都初始化为Yes，表示所有点都是可靠的。
1.1.3二维特征的跟踪
灰度图上的二维特征跟踪同样利用已有的KLT算法完成。其原理是：设检测的特征点在 图象Im，t(p)(m同样为L或R)上，其中p是任一特征点；现在需要在图象Im，t+1上跟踪该点。 基于局部纹理信息的匹配，令 (G的计算如前所述)， 则即为图象Im，t+1上对p的跟踪结果。
实际计算时同样需要对的求解作离散化处理：
跟踪得到的特征点在图2中表示为{pm，t+1(k)}。
1.1.4姿态参数初始化
对于跟踪的结果用和匹配时相同的立体重建方法可获得t+1时刻的三维点序列，这样我 们就获得了这两个时刻(t和t+1)的一组三维对应点{Mp(k)|k＝1，...，K(t)}和 {Mt+1(k)|k＝1，...，K(t)}。下面根据这组三维对应和表示其可靠性的序列{bReliablek}，k＝1，...，K(t) (该序列被称为“可靠度量序列”)来求解第t帧到第t+1帧的姿态变换参数Rt，Tt，同时更 新可靠度量序列。
首先，从未知量的个数来考虑，三对三维对应点就足以计算出姿态参数。简单介绍该已 知算法的原理如下：
设{(Mim，Mim’)|im＝0，1，2}是三对不共线的三维对应点(Mim和Mim’是前后两帧中同一三维 点在L摄像机坐标系中的坐标)，姿态参数Rt，Tt满足Mim′RtMim+Tt|im＝0，1，2。为确定其 中的旋转分量Rt，我们需要确定旋转的中心轴(称为旋转轴，3*1的向量)和绕该轴的旋转 角度θ。定义单位向量 和 则
在得到旋转轴后，θ可以按下式进行求解：
得到了和θ之后，旋转矩阵Rt可由下式给出：
其中 [r]x是根据3*1向量r＝[rx ry rz]T定义的反对称矩阵，即 而 是单位矩阵。
估计出旋转分量Rt之后，平移分量Tt可由下式得到：
上面的算法虽然在理论上能得到准确的结果，但是数值稳定性很差，即在其输入 {(Mim，Mim’)|im＝0，1，2}有较小误差时，算法的求解结果会很不准确。而输入数据一般是有误差 的，甚至还有一些是错误的。为了在这种情况下得到准确的结果，我们采用了常用的基于最 小中值(Least-Median-Square)的鲁棒估计算法。其基本原理是：随机地从当前的可靠三维 对应(即{bReliablek}，k＝1，...，K(t)中取值为Yes的项。设这些项的标号为{ks}，s＝1，...，Kr(t)， 其中Kr(t)是可靠点的总数目)中挑出Num_subset个子集，每个子集Setn，n＝1，...，Num_subset 含有三对三维对应点，据此按上面的算法估计出姿态参数Rt，n，Tt，n，然后对所有的可靠三维对 应{Mt(ks)，Mt+1(ks)}，s＝1，...，Kr(t)按下面的式(1.1.8)计算εn，s：
εn，s＝‖Mt+1(ks)-Rt，nMt(ks)-Tt，n||2 (1.1.8)
该εn，s的物理意义是对应点Mt(ks)，Mt+1(ks)与姿态参数Rt，n，Tt，n的符合程度，该值越小， 说明两者越符合。我们按式(1.1.9)来从集合{Rt，n，Tt，n}，n＝1，...，Num_subset中选择Rt，Tt，即 {Rt，Tt}＝{Rn_m，Tn_m}，其中
其中表示针对下标s取集合{εn，s|s＝1，...，Kr(t)}的中值，它对不同的n构成一个新的 数组 而则是在这个中值数组中选择取最小值的那个n。
上面求解Rt，Tt的算法能保证在存在很多错误的三维对应时仍得到正确的结果，这是因 为：第一，根据{Mt(k)|k＝1，...，K(t)}和{Mt+1(k)|k＝1，...，K(t)}中的三个准确对应就可以计算出姿 态参数；第二，按式(1.1.9)选择姿态参数的依据是整个序列的中值，因此即使有很多三维对 应是错误的，造成该序列里的很多数值较大；只要正确的三维对应数目超过K(t)的一半，算 法仍能得到正确的姿态估计结果。
接着我们需要更新可靠度量序列{bReliablek}，即判断哪些三维对应是正确的。具体方法 如下：令Errt＝{εs＝‖Mt+1(ks)-RtMt(ks)-Tt‖2|s＝1，...，Kr(t)}表示当前的可靠三维点对Rt，Tt的 支持程度，设这Kr(t)个点中的K_True(计算方法如下所述)个是正确的，我们设Errt中前 K_True个较小的εs对应的标号为{su}，u＝1，...，K_True，则中对应于这些标号的 位置取值Yes，其他的位置是No。
K_True的计算步骤是：
(1)设置K_True的初始值为Kr(t)*ratio(ratio为一预先设定的常数)；记Errt中按从小 到大排序的前K_True个元素组成的子序列为E；
(2)计算E的均值 和标准差
(3)在Errt中寻找满足‖εs-ε‖>3*σ的点，设数目为In Valid_T；
(4)如果0≤Kr(t)-K_True-InValid_T≤1，说明当前对误差和错误的判断是相容的，终 止，当前的K_True即为最终结果；否则，转步骤5；
(5)如果K-K_True-InValid_T>1，K_True加1(即E中添加一点)，转步骤2；否 则，转步骤6；
(6)K_True减1(即E中减少一点)，转步骤2。
该方法的原理是误差和错误在统计意义上的分界位置。
1.1.5设置标志bReDetect和bRefine
这是两个布尔型的标志量，取值为Yes或No。bReDetect表示是否要重新调用特征检测 和匹配模块，bRefine则表示是否要进行姿态求精。它们的设置方式如下：
(1)如果t＝N，即当前图象是序列中的最后一幅，则bReDetect＝No，bRefine＝Yes；完毕； 否则，转步骤2；
(2)如果当前Rt的旋转角θt超过了预设的阈值θ_bound，则bReDetect＝Yes， bRefine＝Yes；完毕；否则，转步骤3；
(3)如果可靠三维点的数目Kr(t)(即{bReliablek}，k＝1，...，K(t)中取值为Yes的点的数目) 小于预设的阈值K_bound，则bReDetect＝Yes，bRefine＝Yes；否则，转步骤4；
(4)bReDetect＝No；如果t—t0是预设的常数Inc_refine的倍数，bRefine＝Yes；否则 bRefine＝No。
从原理上说，我们在可靠点数目不足或姿态变化过大时需要重做特征检测和匹配；否则 就做特征跟踪。求精则每隔固定的时间间隔(Inc_refine)就做一次。
现在令t＝t+1，准备进入对下一时刻的处理。
1.1.6姿态求精
bRefine为Yes表示在当前时刻做姿态求精，这时系统还需要设置求精运算的起始时刻 tStart，实际我们取其为t时刻之前第二个设置bRefine为Yes的时刻，如果这样得到的tStart 小于t0就令其为t0。
1.1.4节的姿态参数初始化只利用了两个时刻(t和t+1)的三维对应。为得到更精确的姿 态参数，我们选择那些在连续Num＝t-tStart+1个时刻始终正确的三维对应来求精所有相关 的姿态参数{RtStart+τ，TtStart+τ}，τ＝0，...，t-tStart-1。具体地说，在一次求精运算中，我们先计算 {RatStart+τ，TatStart+τ，τ＝1，...，t-tStart；MtStart(kpp)，pp＝1，...，Kr(t)}，使得式(1.1.10)中的f_obj达到其 最小值。
其中RatStart+τ，TatStart+τ是从tStart时刻到tStart+τ时刻的刚体变换系数； {MtStart(kpp)}，pp＝1，...，Kr(t)，是在这Num个时刻中都可靠的三维点数组，包含Kr(t)个点，pp 表示三维点的标号。m为L或R，表示某一个摄像机。是二维检测或跟踪的结果， 在m号摄像机tStart+τ时刻的图象中，即
而则是pp号三维点在tStart+τ时刻，投影到m摄像机上的结果，即
和 表示从L摄像机坐标系到m摄像机坐标系的刚体变换。即 另外有
从物理意义上看，式(1.1.10)表示的是所有图象中二维重投影的总误差。优化该目标能得 到姿态参数的精确值。其中参数{RatStart+τ，TatStart+r，τ＝1，...，t-tStart；MtStart(kpp)，pp＝1，...，Kr(t)}的 求解可调用通用的Levenberg-Marquardt算法(利用matlab工具包)来完成。
得到了{RatStart+τ，TatStart+τ}，τ＝1，...，t-tStart后按下面的公式计算出 {RtStart+τ，TtStart+τ}，τ＝0，...，t-tStart-1，从而达到求精目的。
RtStart＝RatStart+1，TtStart＝TatStart+1，
上面公式的物理意义是刚体运动的复合，例如：{RtStart+τ-1，TtStart+τ-1}是tStart+τ-1时刻到 tStart+τ时刻的刚体变换，它等于两个刚体变换的复合，即先作从tStart+τ-1时刻到tStart时 刻的刚体变换(是{RatStart+τ-1，TatStart+τ-1}的逆变换)，再作从tStart时刻到tStart+τ时刻的刚体 变换(即{RatStart+τ，TatStart+τ})。
1.2模型和姿态的初始化
按通用模型特定化的思路，首先需要从Faceg得到一个初始模型Faceinit；另外，在上面 的姿态估计步骤中得到了{Rt，Tt，t＝1，...，N-1}，而为了获得模型和投影图象间的对应，我们还需 要计算从generic坐标系(下段介绍，是一个固定在人头上的物体坐标系)到各帧L摄像机坐 标系的刚体变换Rg，t，Tg，t，t＝1，...，N。这两个任务分别由模型初始化和姿态初始化完成。
先解释一下generic坐标系。我们取坐标系原点为左眼内角LE、右眼内角RE、左外嘴唇 LM、右外嘴唇RM四个点的重心Middle；取X轴方向为从右眼RE指向左眼LE的方向，取 Y轴方向从下頦指向额头的方向(竖直向上的方向)，取Z轴方向为从后脑指向鼻子的方向。 三方向两两正交，并构成右手坐标系。上面提到的各特征点参见图5。
模型初始化的过程实际是对原始通用模型施加一个尺度变换。即按下面的公式(1.2.1)完 成：
其中Mg和Minit分别是Faceg和Faceinit上的任一组对应点。三个比例因子sx，sy，sz利用人 脸上的四个主要特征点的距离配置确定(包括左内眼角、右内眼角、左嘴角和右嘴角)。具 体地：sx是建模对象和Faceg中Seg_E长度(两内眼角的距离)的比值，sy则是Seg_EM长 度(两眼中心到嘴中心的距离)的比值，sz就取为sy。这些线段也在图5中标出。
确定了这三个比例系数后对Faceg中的所有点按公式(1.2.1)作尺度变换，结果即为Faceinit。
姿态初始化的任务是要计算Rg，1，Tg，1，即从generic坐标系到第1帧时L摄像机坐标系的刚 体变换，即要求解(1.2.2)中的参数
M1＝Rg，1Minit+Tg，1 (1.2.2)
其中M1和Minit分别是时刻1时人脸和Faceinit这两个模型间的任一组对应点。
为求解Rg，1，Tg，1，我们选择这两个模型在LE(左眼)、RE(右眼)和Mid_M(嘴中心) 这三个位置处的对应点，按前面1.1.4节中姿态初始化的方法得到结果。
获得了Rg，1，Tg，1后，系统根据刚体运动的复合，按下面的公式计算其他各时刻的姿态参数 Rg，t，Tg，t，t＝2，...，N：
Rg，t＝Rt-1Rg，t-1，Tg，t＝Rt-1*Tg，t-1+Tt-1，t＝2，...，N
现在系统会自动判断对应于侧面视图的时刻side，即选择side为使Rg，t的旋转角最接近 90度的那个t。
1.3形状修正A
这一步要把Faceinit变形为Face1。如前所述，我们保持Face的拓扑结构，而只是修改其 中三维点的坐标。具体地，我们先确定三维结点列表P中一个子集Subset的新位置，再通过 一个已知的径向基函数算法来确定P中所有结点的坐标。
先简单介绍径向基函数算法的基本原理。设要将Facestart变形为Facenew，已知一个Subset 中的结点在Facenew中的坐标为{New1，…，News}，记其在Facestart中的对应点为{Sta1，…，StaS}； 这里S为Subset中点的数目。现在对于模型中的任一位置pt，其在Facenew中的坐标可 由下面的已知公式求出：
其中是pt点在Facestart中的坐标，而{C1，…，Cs}可按下面的公式求出
该方法实际是根据Facestart中各点到一些标定点(即Subset中的点)的距离关系来确定该 点在Facenew中的位置。φ是已知的径向基函数(定义域和值域都为非负实数)，它是自变量 r(表示距离)的递减函数，即距离远的标定点对该点的影响较小。
这一步在利用上面的径向基函数算法时，我们有Facestart＝Faceinit，Facenew＝Face1。而 Subset中的点包括两类：一类是三维显著特征点，共4个，由二维显著特征点(指人脸上的 左眼内角LE、右眼内角RE、左外嘴唇LM、右外嘴唇RM四个点在图象平面上的投影)通 过前面的立体重建算法得到；另一类是侧面特征点(指那些在人脸侧面投影视图上易于辨认 的点，如下頦、嘴、鼻尖、鼻梁上方、前额等。典型的人脸侧面投影以及我们定义的一些侧 面特征点在图6中用小的圆圈表示，它的一个实例参见图16)，实际使用的有24个，它们 的三维位置按下面的步骤从手工标出的二维点中获得。
设手工标出的侧面特征点为 kk=1,...,NumSF,其中NumSF＝24 为侧面特征点的数目，side是前面已经得到的侧面帧标号。它们在generic坐标系中的三维坐 标为 则对任一kk，我们有
其中 是姿态参数的分量表示。(1.3.1)中有两个线性方 程，因而可直接用线性方程组的求解方法得到中的两个未知数和
该方法实际是利用了侧面特征点的特殊位置(即在generic坐标系中坐标的X分量为0)， 因而可直接从单幅图象中求出三维点。
1.4形状修正B
这一步同样利用前面的径向基函数算法，Facestart＝Face1，Facenew＝Face2，Subset中除 了包括前面“形状修正A”里已经确定了的点之外，还有一些新添加的点，即根据手工标出 的正面人脸轮廓(指正面姿态下图像上的人脸轮廓，用闭合的多边形轮廓表示，一个例子如 图7所示。该图将多边形轮廓叠加到了人脸图象上，看得更清楚)恢复出来的三维点。得到 这些三维点的基本方法是将当前人脸模型向图象平面作透视投影，并计算三维模型透视投影 的二维轮廓点，最后依赖于这些轮廓点和图象上的人脸轮廓的匹配来计算其三维坐标。具体 的步骤如下所述：
首先需要计算出Face1在t＝1时刻(是手工选择的正面视图)向图象平面投影的轮廓 其中是P(即模型的三维结点列表)中点的标号，显然 有 而nNum1是Cont1中的顶点数目。Cont1中的点就是这一步形状修正中 新添加到Subset中的点，它们是通过判断投影线和模型的交点来选择的，即在某一轮廓点上， 投影线与模型不能有其他的交点，详细的计算方法如下所述：
记三维模型为Faceproj(这里就是Face1)，透视投影中心为 即L摄像机 光心在generic坐标系中的坐标，则我们对模型中的所有顶点计算过该点的投影线和模型的交 点，据此判断该顶点是否在投影轮廓中。具体地，Mnow为Faceproj中任一点，Planejp为Faceproj 中任一面片构成的平面，计算过Center和Mnow的直线(即投影线)与Planejp的交点为Cnow； 如果对某一个满足1≤jp≤nMesh的jp(如前所述，nMesh是Faceproj中的面片总数)求出的 Cnow在线段CMnow的内部，则Mnow不在投影轮廓中；否则就应属于投影轮廓。
接着要计算Cont1中的任一点在Face2的三维坐标。设图象上的人脸轮廓是 它是由二维点列表构成的多边形， 而nNumCI1是点的数目。则我们按下面的方法确定在Face2的三维坐标：
设 点在Faceg中的法向为vn，在Face1中的坐标是 而 经过该点的投影方向为 (这三者都是已知量)，则pti点在Face2中的坐标 满足
其中
v＝(vx vy vz)T＝vpn×(vn×vpn) (1.4.2)
参数tline可根据二维直线 与闭合多边形Cont_img1的交点求 出。具体步骤如下：
设Cont_img1中的任一线段可表示为seg＝(x0+sdx y0+sdy)T，0≤s≤1，其中s为线段参 数，其他都是已知量。
我们先按下面的公式计算
对所有的seg求出这两个数值si和tline。对于很多seg上面的公式无解，即公式里的矩阵 不可逆。有解且得到的si满足0≤si≤1的seg有一个或两个，取与点 距离近的 那一个seg，它得到的tline即为所求。
这里v是点的变形方向，对它的计算利用了通用人脸的局部形状。实际上，该方向 是Face1上该点处法向垂直于投影线的分量，这样一方面能尽量保持局部曲率；另一方面能尽 量减小点的变形距离。
1.5形状修正C
这一步同样利用前面的径向基函数算法，Facestart＝Face2，Facenew＝Face3，Subset中 除了包括前面“形状修正A、B”里已经确定了的点之外，还有一些新添加的点，即根据手工 标出的正面人脸特征轮廓点(指正面姿态下手工标出的眼睛、鼻孔和嘴的轮廓点，是一些离 散点，如图8所示，同样是叠加显示)恢复出来的三维点。具体的恢复步骤如下所述：
设某特征轮廓点的二维坐标是 该点在Face2中的坐标为 我们假设该点的Z坐标已经计算正确，即它在Face3中的Z坐标 等于则它在Face3中的坐标为
即这些特征轮廓点在正面姿态下的投影等于实际图象上的特征点。
1.6形状修正D
这一步同样利用前面的径向基函数算法，Facestart＝Face3，Facenew＝Face4，Subset中 除了包括前面“形状修正A、B、C”里已经确定了的点之外，还有一些新添加的点，即从手 工标出的中间视图1(该视图的选择由系统自动判断，实际是使Rg，t的旋转角最接近25度的 时刻t。所谓“中间视图”是指介于正面和侧面之间的中间投影姿态)人脸轮廓(同样用闭合 多边形表示，典型的例子是图9)恢复出来的三维点，具体的恢复步骤如下所述：
设中间视图1对应于时刻int1，与前面1.4节对正面图人脸轮廓的处理类似，也要先算出 Face3在int1时的投影轮廓 其中是P(即人脸模 型的三维结点列表)中点的标号，而nNum_int1是投影轮廓中的顶点数目。投影轮廓的计算 方法如前所述(1.4节)。
Contint1中的点在Face4中的坐标按下面的与1.4节类似的步骤确定：对于任一 定义 其中 是该点在Face3 中的三维坐标，而 如前所述是L摄像机光心在generic坐标系中的坐标。则 该点在Face4中的坐标满足
其中参数tline2根据二维直线 与闭合多边形 Cont_imgint1的交点求出，方法如1.4节所述。
1.7形状修正E
这一步同样利用前面的径向基函数算法，Facestart＝Face4，Facenew＝Face5，Subset中 除了包括前面“形状修正A、B、C、D”里已经确定了的点之外，还有一些新添加的点，即 从手工标出的中间视图2(该视图的选择也由系统自动判断，实际是使Rg，t的旋转角最接近 40度的时刻t)人脸轮廓(典型的例子是图10)恢复出来的三维点，具体的恢复步骤和1.6 节中的处理完全一样。
1.8纹理映射
纹理映射是要为创建的模型Faces生成柱面纹理图Textures，使得最终的建模结果看起来 有相片一样的真实感。目前的纹理是通过两幅图象生成，即L摄像机拍摄的正面和侧面视图 (IL，1和IL，side)，这部分的基本原理是把采集图象转换到一个统一的柱面坐标系中，并进行融 合，其流程如图3所示，具体步骤如下：
1.8.1柱面展开图的生成
系统首先对Faces中的任一三维点 作公式(1.8.1)所示的柱面映射：
映射的结果是一个二维平面，其上的点为{Cynpt}＝{(lgpt ltpt)}。这里lgpt是柱面上的经 度位置，用弧度值表示；ltpt则是柱面旋转轴向的坐标位置。
1.8.2正面纹理图
将所有在正面姿态参数Rg，1，Tg，1下按公式(1.8.2)进行投影：
其中 是旋转矩阵和平移向量的分量表示。
这时 是Faces中三维点在正面视图IL，1上的二维投影。现在我们得到 一个IL，1平面上点和Cyn平面上点Cynpt＝(lgpt ltpt)T的二维对应关系，该关系是由于它 们都和三维点相对应而得到的(分别通过公式1.8.2和1.8.1)。我们基于这种对应关系 将IL，1中的各像素映射到Cyn平面上，即对Faces中的任一面片facet＝(pt1pt2pt3)T，设 mproj＝Δ(p1，p2，p3)表示这三点在IL，1平面上构成的三角形，mcyn＝Δ(p1′，p2′，p3′)是Cyn平面 上的对应三角形，即 iT1＝1，2，3。对三角形mcyn内的 点[lg lt]按下面的公式(1.8.4)计算[xy]，并令
即将正面纹理图映射到柱面展开图上。
其中的六个未知量aiA，iA＝1，...，6按公式(1.8.5)求出。
正面纹理映射的结果记为Texture1。
1.8.3侧面纹理图
与正面图的处理非常类似，将所有在侧面姿态参数Rg，side，Tg，side下按公式(1.8.6)进行 投影：
这里同样有分量表示
接着对Faces中的任一面片facet＝(pt1 pt2 pt3)T，设msproj＝Δ(ps1，ps2，ps3)，其中 iT2＝1，2，3。对三角形mcyn内的点[lg lt]按后面的公式(1.8.8)计算 [xs ys]，并令
即将侧面纹理图映射到柱面展开图上。
其中的六个未知量asiAS，iAS＝1，...，6按公式(1.8.9)求出。
侧面纹理映射的结果记为Textureside。
1.8.4侧面纹理图的反射
由于人脸只向一侧旋转，所以只能得到单侧的侧面纹理，另一侧的纹理可通过对 Textureside作反射处理得到。具体地说，由于Faces的拓扑结构是完全对称的，设有直接获得 纹理的那一侧脸上的任一三角面片mcyn＝Δ(p1，p2，p3)和它在另一侧脸上的对称面片 mcyn′＝Δ(p1′，p2′，p3′)，其中 iT3＝1，2，3，则对mcyn中 的任一点p＝(lg lt)T按下面的公式(1.8.9)计算[lg′lt′]，并令
即从直接获得纹理的一侧反射生成另一侧。
其中的六个未知量rsiRS，iRS＝1，...，6按公式(1.8.10)求出。
现在我们得到了完整的Textureside。
1.8.5两纹理图的融合
我们融合Texture1和Textureside来得到最终的纹理图Textures。具体地，设定阈值lgmin和 lgmax，Textures在任一位置Cyn＝(lg lt)T处的取值由下面的式(1.8.11)决定
我们在一个普通的桌面PC系统上实现了上面的建模方法。该计算机CPU为Intel PIII933， 配有256M的Rambus内存，安装了Microsoft Windows 2000 Professional操作系统。为进行数 据采集，我们使用了两个相同的SANYO VCC-5972P型摄像机，以及两块相同的Matrox MeteorII采集卡来采集两路同步视频，图象分辨率为768*576，真彩色图象，帧速率15帧/ 秒。数据采集时立体摄像机按图11所示的上下配置(一方面，这种配置使两摄像机的公共视 场较大；另一方面，这时外极线的方向基本沿竖直方向，与人脸上灰度变化多数沿水平方向 的实际情况相适应)，一个实际采集的立体视频序列中的典型帧如图12所示。实际采集的视 频共有50帧(即N＝50)，其中第一帧是手工选择的正面姿态图象。
我们使用的通用模型Faceg来自MPEG4的参考软件MOMUSYS的人脸动画部分，称为 ISTFACE，由葡萄牙的Instituto Superior Tecnico开发。它包含1039(nVertex)个结点，组成 1704(nMesh)个三角形面片，形成一个完整的人头，不但包括人脸表面(例如眼睛、鼻子、嘴、 耳朵等器官和脸颊、额头、后脑等区域)，还包括眼珠、牙齿、舌头这些内部器官和颈部区 域。它可由MPEG4人脸动画参数流驱动，通过三维网格点的变形实现动画。图13是该模型 的中性状态和动画状态(嘴唇的运动)。
下面详细介绍该实施实例中的各个步骤：
第一步：初始化：使用图11所示的摄像机配置拍摄图12所示的立体视频；并设置：
●时间初值t＝t0＝1。其中t表示时间，对应于视频序列中的图象序号，表示当前处理的图 象；t0同样表示时间，表示当前跟踪着的特征被检测出来的时刻；
●计算基础矩阵 (其中各已知量的意义参见1.1.2节中 的解释)。实际数值是 其中 fL＝1330.3，fR＝1330.6，u0L＝384，v0L＝288，u0R＝384，v0R＝288；
●width，height：图象的宽度和高度，在具体实现中，width＝768，height＝576；
●blockx，blocky：决定了特征检测、跟踪以及匹配时局部窗口中包含的网格点数目。实际 上，该窗口以待测点(i，j)为中心，网格点总数为(2*blockx+1)*(2*blocky+1)。在具体实现中， blockx＝3，blocky＝3；
●quality，dis_min：特征点检测时的阈值，前者控制特征点的质量，后者控制特征点的间 距。在具体实现中，quality＝0.01，dis_min＝10；
●Num_subset，ratio：姿态参数初始化(1.1.4节)时使用的常数，前者是随机选择的子 集的数目，后者用来决定如何更新{bReliablek}。实际取值为Num_subset＝100，ratio＝0.9；
●θ_bound，K_bound，Inc_refine：是我们为计算bReDetect和bRefine而使用的阈值。它 们的取值是：θ_bound＝20，单位是度；K_bound＝10；Inc_refine＝4；
●lgmin和lgmax：用来确定Texture1和Textureside之间的分界点，是外眼角和耳朵上沿之间 的一个位置。具体地，lgmin对应于右侧脸上的该点，在Faces中标号为306；lgmax对应于左侧 脸上的该点，在Faces中标号为120；
●φ：即形状修正时已知的径向基函数，我们选择 其中表示函数取值为负时就置其为0。
第二步：在图象IL，t上计算梯度向量
以及
其中1≤i≤width，1≤j≤height表示图象上任一点，blockx，blocky是与窗口尺寸相关的常 量。如果计算时用到了超出图象边界的点，直接令 接着利用已有的数学工 具软件matlab对GL，t(i，j)作奇异值分解，并取其最小的奇异值组合得到 {min_val(i，j)}，1≤i≤width，1≤j≤height，这是与原图象尺寸相同的浮点数矩阵。
举一个实际例子。在图12所示立体视频序列的IL，1上，以某点(i＝346，j＝297)为中心的 7*7(blockx＝blocky＝3)局域图象上每一点的灰度值为
而7*7的{gx}和{gy}矩阵为
按前述公式计算出的2*2的矩阵 其两个奇异值为31561.2和 22449.8，故min_val(346,297)＝22449.8(这里min_val是两个奇异值中较小的一个，函数的 自变量是图象上点的二维坐标)。
第三步：在矩阵{min_val(i，j)}中寻找满足下面几个条件的位置作为检测出的特征点。(1) 该点处的min_val取值较大；具体地，设整个矩阵的最大值为max_min_val，则特征点处的取 值应大于quality*max_min_val，其中quality为一预先设定的阈值，实际取为0.01；(2)该点 处的min_val是其3*3邻域内的最大值；(3)任意两特征点的距离应大于dis_min：如果两特 征点的距离小于或等于该值，则min_val取值较小的一个位置不被认为是特征点；dis_min实 际取为10。按这三条规则我们检测出IL，t上的特征点序列，并记为{pL，t(k)}，1≤k≤K(t)。检测 的结果参见图14。 第四步：对任一pL，t(k)，计算其对应的外极线并计算
其中(i，j)是直线l(其向量表示是l＝[l1 l2 l3]，即直线的方程是l1i+l2j+l3＝0)上的 任一点，即满足 则pL，t(k)的匹配点即为 pR，t(k)＝(imin jmin)T。对所有检测出的特征点进行相同的处理，我们得到IR，t上的匹配点序列， 并记为{pR，t(k)}，1≤k≤K(t)。
第五步：根据任一对匹配pL，t(k)和pR，t(k)，计算 并最终得到三维点 对 所有成对匹配进行相同处理，我们得到t时刻的三维特征点序列{Mt(k)|k＝1，...，K(t)}。同时初 始化布尔变量序列{bReliablek}，k＝1，...，K(t)全为Yes，表示所有点都是可靠的。
第六步：在任一pL，t(k)＝(i，j)处计算 以及 则该点在L视频序列中t+1时刻的对应点为 对所有的特征点进行相同的处理，我们得到IL，t+1上的跟踪结果， 并记为{pL，t+1(k)}，1≤k≤K(t)。
第七步：用类似的方法处理R视频序列。具体地，先按前述公式在R，t图上计算 和{GR，t(iR，jR)}，接着对任一pR，t(k)＝(iR，jR)计算 最终有 则该点在R视频序列中t+1时刻的对应点为 对所有特征点进行相同的处理，我们得到IR，t+1上的跟踪结果， 并记为{pR，t+1(k)}，1≤k≤K(t)。
第八步：用与第五步类似的方法得到t+1时刻的三维点序列。具体地，根据任一对匹配 pL，t+1(k)和pR，t+1(k)，计算 并最 终得到三维点 对所有成对的跟踪结果 进行相同处理，我们得到t+1时刻的三维特征点序列{Mt+1(k)|k＝1，...，K(t)}。
第九步：设序列{bReliablek}，k＝1，...，K(t)中取值为Yes的项的标号为(ks}，s＝1，...，Kr(t)， 其中Kr(t)是取值为Yes的总数目。随机选取Num_subset(实际取为100)个三元子集 {Setn＝{N1，n，N2，n，N3，n}}，n＝1，...，Num_subset，其中1≤N1，n，N2，n，N3，n≤Kr(t)。根据 和 的对应关系，按1.1.4节的方法 估计出刚体运动参数Rt，n，Tt，n，然后对所有的可靠三维对应{Mt(ks)，Mt+1(ks)}，s＝1，...，Kr(t)计 算εn，s＝‖Mt+1(ks)-Rt，nMt(ks)-Tt，n‖2，n＝1，...，Num_subset；s＝1，...，Kr(t)，并按下面的公式得到 {Rt，Tt}：
{Rt，Tt}＝{Rn_m，Tn_m}，其中
第十步：计算序列 并更新 {b Reliableks}，s＝1，...，Kr(t)。具体地，我们求出Errt中前K_True个较小的εs对应的标号 {su}，u＝1，...，K_True，然后把中对应于这些标号的位置处置为Yes，其他的置为 No。K_True的具体计算步骤如前所述。
第十一步：按1.1.5节的四步骤方法设置两个布尔型的标志量bReDetect和bRefine，并令 t＝t+1。
第十二步：如果bRefine＝No，转第十三步；否则，令tStart为t时刻之前第二个设置bRefine 为Yes的帧，如果这样得到的tStart小于t0就令其为t0。然后调用通用的Levenberg-Marquardt 算法(利用matlab工具包)求解式(1.1.10)，得到 {RatStart+τ，TatStart+τ，τ＝1，...，t-tStart；MtStart(kpp)，pp＝1，...，Kr(t)}，再按下面的公式计算出 {RtStart+τ，TtStart+τ}，τ＝0，...，t-tStart-1：
RtStart＝RatStart+1，TtStart＝TatStart+1，
第十三步：如果bReDetect＝Yes，则令t0＝t，转第二步；否则，令t＝t+1，如果这时t＝N， 则转第十四步，否则转第六步；
上面各步骤完成了姿态估计(1.1节的内容)。在对图12的输入进行处理时，特征检测 和匹配在第1，15，29，38帧处(即t0的不同取值)，实际结果如图14所示。每次检测和匹配 上的全部特征数目分别是89、97、107和61(即各个K(t0))，在t0保持不变的所有时刻(即 两次检测之间的特征跟踪)里一直可靠的特征点(即序列{bReliablek}中取值为Yes点的数目) 分别有12、18、36、14个。
第十四步：在IL，1和IR，1上手工选择左内眼角、右内眼角、左嘴角和右嘴角这四个二维显 著特征点，实际结果如图15所示。根据这四个位置处的两幅图间的匹配，可按前述的立体重 建算法得到四个三维显著特征点，即左内眼角右内眼角左嘴角右嘴角 然后按公式(1.2.1)将原始通用模型Faceg变形为Faceinit。
其中三个比例因子sx，sy，sz的计算方法如下：令 即建模对象和 Faceg中Seg_E长度(两内眼角的距离)的比值；令 即Seg_EM长度(两眼中心到嘴中心的距离)的比值；令sz＝sy。其中 是在Faceg上这四个主要特征点的三维坐标。
第十五步：接着我们求解generic坐标系到L摄像机坐标系的刚体变换Rg，1，Tg，1。我们选 择 和 这三个三维对应点，按前述的 姿态参数初始化一节的方法求解。
获得了Rg，1，Tg，1后，系统按下面的公式计算Rg，t，Tg，t，t＝2，...，N：
Rg，t＝Rt-1Rg，t-1，Tg，t＝Rt-1*Tg，t-1+Tt-1，t＝2，...，N
接着系统会自动判断对应于侧面视图的帧side，即选择side为使Rg，t的旋转角最接近90 度的时刻t。我们实际得到的side＝50，对应的旋转角度是86.9度。
第十六步：在图象IL，side上手工拾取侧面特征点，实验结果如图16所示，记为 kk＝1，...，NumSF，其中NumSF为侧面特征点的数目，实际为24。
第十七步：按前述方法求解式(1.3.1)，得到三维侧面特征点
第十八步：按前述的1.3节中的径向基函数算法把Faceinit变形为Face1。Subset中的点包 括两类：一类是三维显著特征点，共4个，即另一类是侧面特征点， 即十七步求出的kk＝1，...，NumSF。
第十九步：按前述1.4节中计算投影轮廓的方法得到Face1在时刻1时投影到图象平面的 轮廓点 Face1就是那里的Faceproj。
第二十步：手工标出IL，1上的人脸轮廓，如图7所示。这是一个闭合的多边形，记为 是连续的二维顶点列表，nNumCI1是 二维顶点的数目。
第二十一步：按照前述的公式(1.4.1)确定Cont1中点在Face2中的三维坐标
第二十二步：按前述1.3节中的径向基函数算法把Face1变形为Face2。Subset中新添加 的点为第二十一步求出的
第二十三步：手工标出IL，1图上的正面人脸特征轮廓点(实验结果为图8)，结果为 其中nNumCF是点的总数目。
第二十四步：按前述的公式(1.5.1)求解第二十三步标出的特征轮廓点的三维坐标
第二十五步：按前述1.3节中的径向基函数算法把Face2变形为Face3。Subset中新添加 的点为第二十四步求出的
第二十六步：系统自动判断对应于中间视图1的帧int1，即选择int1为使Rg，t的旋转角最 接近25度的时刻t。接着在IL，int1图上的手工标出人脸轮廓，如图9所示。同样用闭合多边形 表示，记为 nNumCI int1是顶 点数。
第二十七步：按前述1.4节中计算投影轮廓的方法得到Face3在时刻int1时投影到图象平 面的轮廓点 Face3就是那里的Faceproj。
第二十八步：根据Cont_imgint1确定Contint1中点在Face4中的三维坐标。具体地，对于任 定义 其中 是 该点在Face3中的三维坐标。则该点在Face4中的坐标满足 其中 参数tline2根据二维直线 与闭合多边形Cont_imgint1的交点求 出，方法如1.4节所述。对Contint1中所有点按上面的步骤计算，得到
第二十九步：按前述1.3节中的径向基函数算法把Face3变形为Face4。Subset中新添加 的点为第二十八步求出的
第三十步：系统自动判断对应于中间视图2的帧int2，即选择int2为使Rg，t的旋转角最接 近40度的时刻t。接着在IL，int2图上的手工标出人脸轮廓，如图10所示。同样用闭合多边形 表示，记为 nNumCIint2是 顶点数。
第三十一步：按前述1.4节中计算投影轮廓的方法得到Face4在时刻int2时投影到图象平 面的轮廓点 Face4就是那里的Faceproj。
第三十二步：根据Cont_imgint 2确定Contint 2中点在Face5中的三维坐标。具体地，对于任 一 定义 其中 是该点在Face4中的三维坐标。则该点在Face5中的坐标满足 其 中参数tline3根据二维直线 与闭合多边形Cont_imgint2的交点 求出，方法如1.4节所述。对Contint2中所有点按上面的步骤计算，得到
第三十三步：按前述1.3节中的径向基函数算法把Face4变形为Face5。Subset中新添加 的点为第三十二步求出的 Face5即为最终的形状模型Faces。
这样我们一共进行了五次形状修正，对应于1.3至1.7节的形状修正A—E，实际上各步 使用的Subset中的点数目如下面的表1所示。
不同模型 新加入的Subset点及其数目 Subset点总数 Face1 4个正面显著特征(两内眼角、两嘴角) +24个侧面特征 28 Face2 正面视图的投影轮廓点(41个) 69 Face3 正面视图中两眼、两鼻孔和嘴的特征轮廓 点(44个) 113 Face4 中间视图1时的投影轮廓点(37个) 150 Faces 中间视图2时的投影轮廓点(44个) 194
表1 形状修正中每一步新加入的Subset点的位置及数目
第三十四步：对Faces中的任一三维点 按公式 作柱面映射，结果为{Cynpt}＝{(lgpt ltpt)}。
第三十五步：将所有在正面姿态参数Rg，1，Tg，1下按公式 进行投 影，得到 然后按公式(1.8.3)得到正面纹理图Texture1，实际例子如图17 所示。
第三十五步：与第三十四步类似，获得Faces的侧面纹理图Textureside。具体地，将所有 在侧面姿态参数Rg，side，Tg，side下按公式 进行投影，得到 然后按公式(1.8.7)得到侧面纹理图Textureside，实际例子如图18所示。
第三十六步：在柱面图{(lg，lt)}上进行纹理反射。实际实验中直接获得纹理的那一侧脸是 左侧，因此这一步根据左侧的纹理信息按公式(1.8.8)生成右侧脸。这时完整的Textureside如图 19所示。
第三十七步：我们按公式融合Texture1和 Textureside来得到最终的纹理图Textures，实际结果如图20所示。
我们最终的、映射了Textures的建模结果Faces如图21所示。
| 标题 | 发布/更新时间 | 阅读量 |
|---|---|---|
| 一种三维京剧脸谱自动化妆方法 | 2020-05-20 | 179 |
| 一种利用多目3D视觉传感技术防控儿童近视的系统及方法 | 2020-05-21 | 624 |
| 一种三维人脸的人脸级联回归模型训练方法及重建方法 | 2020-05-21 | 937 |
| 一种基于大姿态对准的人脸识别方法及系统 | 2020-05-12 | 419 |
| 一种高压氧舱人脸识别交互方法及系统 | 2020-05-15 | 388 |
| 一种基于单张图片的三维头发重建方法 | 2020-05-08 | 448 |
| 一种基于3D打印技术精准修复遗体面部残缺的方法 | 2020-05-14 | 586 |
| 一种实时人脸表情重建方法 | 2020-05-15 | 723 |
| 用于生成信息的方法和装置 | 2020-05-08 | 508 |
| 三维表情基的生成方法、装置、设备及存储介质 | 2020-05-11 | 696 |
高效检索全球专利专利汇是专利免费检索,专利查询,专利分析-国家发明专利查询检索分析平台，是提供专利分析，专利查询，专利检索等数据服务功能的知识产权数据服务商。
我们的产品包含105个国家的1.26亿组数据，免费查、免费专利分析。
专利汇分析报告产品可以对行业情报数据进行梳理分析，涉及维度包括行业专利基本状况分析、地域分析、技术分析、发明人分析、申请人分析、专利权人分析、失效分析、核心专利分析、法律分析、研发重点分析、企业专利处境分析、技术处境分析、专利寿命分析、企业定位分析、引证分析等超过60个分析角度，系统通过AI智能系统对图表进行解读，只需1分钟，一键生成行业专利分析报告。
