融合多视角、多线索二维信息的人脸三维模型的建立方法属于计算机人脸图像重建技术领域。

现有的根据图象创建三维人脸模型的方法缺乏从多种二维线索中提取三维信息并进行融合的策略。Kang，Sing Bing等人在“image method and apparatus for semi-automatically mappinga face on to a wireframe topology”(United States Patent 6,031,539)的专利技术中提出了一种根据单幅正面图象创建人脸模型的半自动方法。该方法的原理是：自动地检测出人脸图象上的多个特征点，进而根据特征点定位面部的界标点(landmark points)，然后据此得到人脸线框模型中所有网格点的位移，创建出与该图象相对应的人脸模型。该方法的主要缺点是没有充分利用各视角的信息，这限制了模型的精确程度。
在“利用桌面摄像机克隆人脸”(《第八届计算机视觉国际会议论文集》，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)：G(i,j)=Σii=-blockxblockxΣjj=-blockyblockyg→(i+ii,j+jj)g→T(i+ii,j+jj),]]>为待测点p处图像的梯度向量：g→(i,j)=gxgyT=[I(i+1,j)-I(i-1,j)2I(i,j+1)-I(i,j-1)2]T,]]>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(x)yL，t(k))T，它的匹配点pR，t(k)＝(xR，t(k)yR，t(k))T应位于一条由pL，t(k)决定的直线上，在这个一维空间中，选择使一个局部窗口W内总灰度差异最小的点作为pL，t(k)的匹配点pR，t(k)，总灰度差异的表达式：Diff(i,j)=Σii=-blockxblockxΣjj=-blockyblocky(IL,t(xL,t(k)+ii,yL,t(k)+jj)-IR,t(i+ii,j+jj))2]]>(i，f)是直线 F上的任一点；记时刻t的匹配结果为{pR，t(k)}，总数目为K(t)个，每点都与pL，t(k)相对应，在下述公式表示的直线上搜索，使Diff(i，j)最小的位置：pL,t(k)1TFpR,t(k)1=0]]>，其中：F=fL0uOL0fLvOL001-T0-Tc,zTc,yTc,z0-Tc,x-Tc,yTc,x0RcfR0uOR0fRvOR001-1]]>，fL，fR，uOL，vOL，uOR，vOR为摄像机内参数，已知量；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分别是同一点在这两个摄像机坐标系中的坐标；以及摄像机透视投影公式：smAm-1pm,t1=Mm,t]]>，sm(即sL或sR)是待定的比例因子(由稍后的公式求出)，得到三维点在摄像机L、R中的坐标ML，t和MR，t：ML,t=12(sLAL-1pL,t1+sRRcAR-1pR,t1+Tc);]]>AL=fL0uOL0fLvOL001,]]>AR=fR0uOR0fRvOR001;]]>SLSR=(BTB)-1BTTc,]]>而B=[AL-1pL,t1-RcAR-1pR,t1:]]>MR,t=Rc-1(ML,t-Tc);]]>
在姿态估计的余下部分，三维点都用L摄像机坐标系中的坐标来表示，因此省略下标L，即将ML，t简记为Mt；(2.2)，后继时刻：对t+1时刻的立体图对用KLT算法进行二维特征跟踪，再同样通过立体重建得到各后继时刻的三维点；跟踪是指同一摄像机拍摄的不同时刻图像间的二维对应：在二维特征跟踪时，设检测的特征点在图像Im，t(p)(m为L或R)上，p(x＝i，y＝j)是任一特征点，则其在图像Im，t+1上的跟踪结果为 ，其中d→=G-1e→]]>，G已如上述，而e→(i,j)=Σii=-blockxblockxΣjj=-blockyblockyg→(i+ii,j+jj)-(Im,t+1(i+ii,j+jj)-Im,t(i+ii,j+jj))]]>记跟踪得到的特征点为{pm，t+1(k)}，m表示L或R，t+1为时刻：对于摄像机L为：pL,t+1(k)=pL,t(k)-d→L,t(pL,t(k)),1≤k≤K(t);]]>对于摄像机R为：pR,t+1(k)=pR,t(k)-d→R,t(pR,t(k)),1≤k≤K(t);]]>根据任一对匹配的跟踪点计算三维点{Mt+1(k)|k＝1，…K(t)}(如前所述，这里省略了下标L)；Mt+1(k)=12(sLAL-1pL,t+1(k)1+sRRcAR-1pR,t+1(k)1+Tc);]]>其中SLSR=(BTB)-1BTTc,]]>B=[AL-1pL,t+1(k)1-RcAR-1pR,t+1(k)1]:]]>(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))，每个三元子集中三维对应点为：{Mt(kN1,n),Mt(kN2,n),Mt(kN3,n)}]]>和{Mt+1(kN1,n),Mt(kN2,n),Mt+1(kN3,n)}]]>，可简记为{(Mim，Mim’)|im＝0，1，2}，它表示三对不共线的三维对应点，Rt，Tt满足Mim′＝RtMim+Tt|im＝0，1，2；然后从下列公式中求出R和T(略去下标t)：R=I+sin||r||||r||[r]x+1-cos||r||||r||2[r]x2,]]>其中r=θ*r^]]>,θ=arccos(r^×dir1)·(r^×dir1′)||(r^×dir1)·(r^×dir1′)||,]]>r^=(dir1′-dir1)×(dir2′-dir2)||(dir1′-dir1)×(dir2′-dir2),]]>是旋转的中心轴，称旋转轴，3*1的向量；θ是绕 的旋转角；dirim=Mim-M0||Mim-M0||,im=1,2;]]>dirim′=Mim′-M0′||Mim′-M0′||,im=1,2;]]>[r]x是根据3*1向量r＝[rxryrz]T定义的反对称矩阵：[r]x=rxryrzx=0-rzryrz0-rx-ryrx0,]]>而I=100010001]]>是单位矩阵；平移分量T由下式估计出：T=13Σim=02(Mim′-RMim);]]>对每个子集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}，n_m=argminnmedsn{ϵn,s|s=1,…,Kr(t)},]]>表示针对下标s取集合{εn，s|s＝1，…，Kr(t)}的中值，它对不同的n构成一个新的数组{medsn|n=1,…,Num_subset}]]>；而argmin则是在这个中值数组中选择取最小值的那个n；(2.4)，在t＝t+1时刻，设置：是否需要重新调用特征检测和匹配模块的布尔标志量bReDetect；是否需要进行姿态求精的布尔标志量bRefine；(2.5)，当t＝N时，姿态估计结束；(3).模型初始化和姿态初始化：模型初始化是通用模型特定化方法的初始步骤，它从通用模型Faceg得到一个初始化模型Faceinit，即对通用模型进行尺度变换，可用下式进行：Minit=sx000sy0OOszMg]]>，Mg和Minit分别是Faceg和Faceinit上的任一组对应点；sx是建模对象和Faceg的两内眼角距离Seg_E长度的比值，sy是两眼中心到嘴中心距离Seg_EM长度的比值；sz＝sy；再根据此对Faceg中所有点按上式变换；姿态初始化是为了获得模型和图象间的对应，即计算从人脸generic坐标系到L摄像机坐标系的刚体变换；即从下式求解Rg，l，Tg，l：Ml＝Rg，lMinit+Tg，l，Ml和Minit分别是时刻t＝1时人脸和Faceinit间的任一组对应点；本发明取上述两个模型在LE(左眼)、RE(右眼)和Mid_M(嘴中心)这三个位置处的对应点，按上述r、θ、R、[r]x、T各式依次求得Rg，l和Tg，l；再按下式从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)，手工拾取侧面特征点：设：手工标出的侧面特征点为{mL,sidekk,=xL,sidekkyL,sidekkT},]]>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中的坐标 可由下面的已知公式求出：Mptnew=ΣjRBF=1SCjRBFφ(||Mptstart-StajRBF||);]]>是pt点在Facestart中的坐标，而{C1，…，CS}可按下面的公式求出C1T...CST=φ(||Sta1-Sta1||)...φ(||Sta1-StaS||).........φ(||StaS-Sta1||)...φ(||StaS-StaS||)-1New1T...NewST]]>；φ是已知的径向基函数；“形状修正A”在利用上面的径向基函数算法时，设置Facestart＝Faceinit，Facenew＝Face1。而Subset中的点包括两类：一类是人脸上的左眼内角LE、右眼内角RE、左外嘴唇LM、右外嘴唇RM四个三维显著特征点，由二维显著特征点通过上述立体重建方法得出；另一类是上述24个侧面特征点，它们的三维位置按下面的步骤从手工标出的二维点中获得：设：侧面特征点在人脸generic坐标系中的三维坐标为{Msgkk=0YgkkZgkkT},]]>kk＝1，…，NumSF；则对于任一kk，有其中Rg,side=Rg,side,1TRg,side,2TRg,side,3T,]]>Tg,side=Tg,side,1Tg,side,2Tg,side,3]]>是姿态参数的分量表示，为已知量；可直接用线性方程组的求解方法得到 中的两个未知数 和 (4.3)，形状修正B：按前述的径向基函数算法把Face1变形为Face2；设Facestart＝Face1，Facenew＝Face2，Subset中除了包括“形状修正A”里已经确定了的点之外，还有一些新添加的点，即从正面姿态下手工标出的人脸轮廓恢复出来的三维点，得到这些三维点的基本方法是把当前人脸模型Face1向正面姿态下的图象平面作透视投影并计算三维模型投影结果的二维轮廓点，最后依赖于这些轮廓点和上述图象上的人脸轮廓即闭合多边形的匹配来计算其三维坐标：首先，计算出Face1在t＝1时刻向手工选择的正面视图的图象平面上的投影轮廓Cont1={pt1iC1|iC1=1,…,nNum1},]]>是模型三维结点列表P中点的标号，1≤pt1iC1≤nVertex;]]>而nNum1是Cont1中顶点的数目，它们即形状修正B中向Subset新添加的点；它们是通过判断投影线与模型的交点来选择的，即在某一轮廓点上，投影线与模型不能有其他交点，其算法如下：Face1为三维模型，透视投影中心为Center=-Rg,1-1Tg,1]]>已知量，则对模型中所有顶点计算过Center点的投影线和模型的交点，据此判断该顶点是否在投影轮廓上，若Mnow为Face1中任一点，Planejp为Face1中任一面片构成的平面，计算过Center点和Mnow的直线，即投影线，与Planejp的交点为Cnow；如果对某一面片1≤ip≤nMesh而言，点Cnow在线段CMnow的内部，则Mnow不在投影轮廓上；否则即为投影轮廓点；再计算Cont1中的任一点 在Face2的三维坐标：设上述图象上的人脸轮廓是Cont_img1={xCI1iCimgyCI1iCimg|iCimg=1,...,nNumCI1}]]>，它是由二维点列表构成的多边形，nNumCI1是点的数目；则pti=pt1iC1]]>在Face2的三维坐标 M2pti=M1pti+tlinev;]]>v＝(vxvyvz)T＝vpn×(vn×vpn)；其中，M1pti=X1ptiY1ptiZ1ptiT]]>是pti点在Face1中的坐标；经 点的投影方向为vpn为已知量；vn为pti点在Faceg中的法向；参数tline可根据二维直线X1ptitlinevxZ1pti+tlinevzY1pti+tlinevyZ1pti+tlinevzT]]>与闭合多边形Cont_img1的交点求出：设Cont_img1中的任一线段为seg＝(x0+sdxy0+sdy)T，0≤s≤1，s为线段参数，其他都是已知量；先按下面的公式计算sitline=dx-(X1pti-vx*Z1pti/vz)dy-(Y1pti-vy*Z1pti/vz)-1x0-vx/vzy0-vy/vz]]>对所有的seg求出这两个数值si和tline；对于很多seg上面的公式无解，即公式里的矩阵不可逆；有解且得到的si满足0≤si≤1的seg有一个或两个，取与点X1ptiZ1ptiY1ptiZ1ptiT]]>距离近的那一个seg，它得到的tline即为所求；(4.4)，形状修正C：按前述的径向基函数算法把Face2变形为Face3；设Facestart＝Face2，Facenew＝Face3，Subset中除了包括“形状修正A、B”里已经确定了的点之外，还有一些新添加的点，即从手工标出的正面人脸特征轮廓点：眼睛、鼻孔和嘴的轮廓点恢复出来的三维点：设某特征轮廓点的二维坐标是xCF1iCFyCF1iCF]]>，该点在Face2中的坐标为M2iCF=X2iCFY2iCFZ2iCFT]]>，设该点的Z坐标已经计算正确，即它在Face3中的Z坐标等于，则它在Face3中的坐标为M3iCF=xCF1iCF*Z2iCFyCF1iCF*Z2iCFZ2iCFT]]>，即这些特征轮廓点在正面姿态下的投影等于实际图象上的特征点；(4.5)，形状修正D：按前述的径向基函数算法把Face3变形为Face4；设Facestart＝Face3，Facenew＝Face4，Subset中除了包括“形状修正A、B、C”里已经确定了的点之外，还有一些新添加的点，即从手工标出的中间视图1即使Rg，t的旋转角最接近25度的时刻t的人脸轮廓中恢复出来的三维点，具体步骤如下：设该视图对应于时刻int1，与形状修正B所述的相同，先算出Face3在int1时的投影轮廓Contint1={ptint1iC2|iC2=1,…,nNum_int1},]]>是P中点的标号，nNum_int1是投影轮廓中的顶点数目；对于任一点ptli=ptint1iC2]]>，它在Face3中的三维坐标M3ptli=X3pt1iY3pt1iZ3pt1iT]]>，而Center=-Rg,1-1Tg,1]]>如前所述是L摄像机光心在generic坐标系中的坐标，则该点在Face4中的坐标应满足M4ptli=M3ptli+tline2v2,]]>v2=v2xv2yv2zT=M3ptli-Center,]]>参数tline2同样根据直线X3pt1i+tline2v2xZ3pt1i+tline2v2zY3pt1i+tline2v2yZ3pt1i+tline2v2zT]]>与闭合多边形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中的任一三维点Mspt=XsptYsptZsptT]]>按下式作柱面映射：，映射的结果是一个二维平面，其上的点为{Cynpt}＝{(lgptltpt)}；这里lgpt是柱面上的经度位置，用弧度值表示；ltpt则是柱面旋转轴向的坐标位置；(5.2)，生成正面纹理图：把所有点 在已知的正面姿态参数Rg，l，Tg，l下按下式在正面视图IL，1上投影，p→1pt=x1pty1ptT]]>是Faces中三维点 在IL，1上的二维投影：，其中Rg,1=R1g,1TR2g,1TR3g,1T,]]>Tg,1=T1g,1T2g,1T3g,1]]>是旋转矩阵和平移向量的分量表示；接着，把IL，1中的各像素映射到Cyn平面上，即对Faces中的任一面片facet＝(pt1pt2pt3)T，mproj＝Δ(p1，p2，p3)表示这三点在IL，1平面上构成的三角形，mcyn＝Δ(p1′，p2′，p3′)是Cyn平面上的对应三角形，即piT1=x1ptiT1y1ptiT1T,]]>piT1′=lgptiT1ltptiT1T]]>iT1＝1，2，3；再对三角形mcyn内的点[lglt]按下式计算[xy]：xy1=a1a2a3a4a5a6001-1lglt1;]]>a1～a6由下式求出：再令正面纹理图ITexture1(lg,lt)=IL,1(x,y);]]>(5.3)，生成侧面纹理图：把所有点 在已知的侧面姿态参数Rg，side，Tg，side下按下式在侧面视图IL，side上投影，p→sidept=xsideptysideptT]]>是Faces中三维点 在IL，side上的二维投影：Rg,side=R1g,sideTR2g,sideTR3g,sideT,]]>Tg,side=T1g,sideT2g,sideT3g,side;]]>再对Faces中的任一面片facet＝(pt1pt2pt3)T，设msproj＝Δ(ps1，ps2，ps3)，psiT2=xsideptiT2ysideptiT2T]]>，iT2＝1，2，3；再对三角形mcyn内的点[lg lt]按下式计算[xs ys]：xsys1=as1as2as3as4as5as6001-1lglt1;]]>as1～as6由下式求出：再令侧面纹理图ITextureside(1g,lt)=IL,side(xs,ys);]]>(5.4)，反射侧面纹理图以得到另一侧的纹理图，这是因为Faces的拓扑结构是完全对称的；设：直接获得纹理的那一侧脸上的任一三角面片mcyn＝Δ(p1，p2，p3)，它在另一侧脸上的对称面片是mcyn′＝Δ(p1′，p2′，p3′)，其中piT3=lgptiT3ltptiT3T,]]>piT3′=lg′ptiT3lt′ptiT3T,]]>iT3＝1，2，3；对mcyn中的任一点p＝(lg lt)T按下式计算[lg′lt′]，lg′lt′1=rs1rs2rs3rs4rs5rs6001-1lglt1;]]>rs1～rs6按下式求出：再从直接有纹理的一侧反射生成另一侧：ITextureside(lg′,lt′)=ITextureside(lg,lt);]]>(5.5)，正面的和侧面的纹理图的融合，得到最终纹理图Textures；设定阈值lgmin和lgmax，Textures在任一位置Cyn＝(lg lt)T处的取值按下式给出：2.所述的姿态估计步骤中的更新可靠度量序列{bReliablek}的方法，它依次含有如下步骤：(1).根据对序列Errt的统计分析计算正确对应的数目K_True；(1.1)，设：K_True的初始值为Kr(t)*ratio，ratio为预先设定的常数；Errt＝{εs＝‖Mt+1(ks)-RtMt(ks)-Tt‖2＝1，…，Kr(t)}中按从小到大排序的前K_True个元素组成的子序列为E；(1.2)，计算E的均值ϵ‾=1K_TrueΣϵs∈Eϵs]]>和标准差σ=1K_True-1Σϵs∈E(ϵs-ϵ‾)2;]]>(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，把{bReliablek}中对应于这些标号的位置赋值为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)(即{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.所述的姿态估计步骤中，bRefine为Yes表示在当前时刻作姿态求精，即得出更精确的姿态参数：设选择在连续Num＝t-tStart+1个时刻始终正确的三维对应来求精所有相关的姿态参数{RtStart+τ，TtStart+τ}，τ＝0，…，t-tStart-1，tStart为求精运算的起始时刻，现取t时刻之前的第二个设置bRefine为Yes的时刻，若tStart＜t0，取其为t0；在一次求精运算中，先计算{RatStart+τ，TatStart+τ，τ＝1，…，t-tStart；MtStart(kpp)，pp＝1，…，Kr(t)}，使所有图像中二维重投影的总误差f_obj达到最小值：f_obj=ΣmΣτ=0t-tStartΣpp=1Kr(t)[(app,mτ-xpp,mτ)2+(bpp,mτ-ypp,mτ)2];]]>其中RatStart+τ，TatStart+τ是从tStart时刻到tStart+τ时刻的刚体变换系数；{MtStart(kpp)}，pp＝1，…，Kr(t)，是在这Num个时刻中都可靠的三维点数组，共Kr(t)个，pp为三维点的标号；m为L或R，表示摄像机； 是在m号摄像机、tStart+τ时刻图象中的二维检测或跟踪的结果，即app,mτbpp,mτ1=Am-1pm,tStart+τ(kpp)1,]]>Am为AL或AR； 是pp号三维点在tStart+τ时刻，投影到m摄像机上的结果，即，其中Rc,m=R1c,mTR2c,mTR3c,mT]]>和Tc,m=T1c,mT2c,mT3c,m]]>是从L摄像机坐标系到m摄像机坐标系的刚体变换，即Rc,L=100010001,]]>Tc,L=000,]]>Rc,R=Rc-1,Tc,R=-Rc-1Tc,R]]>；另外，RatStart=100010001,]]>TatStart=000;]]>
优化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，，τ＝2，…，t-tStart试验证明：它达到了预期目的。
附图说明
图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(从不同视角查看)具体实施方式我们在一个通用模型特定化的框架下创建人脸模型，即：从通用模型Faceg开始，根据获得的建模对象的二维信息逐步进行修改，最终得到特定人模型Faces。在这一过程中，修改的只是P中的三维点坐标，而F则保持不变，即拓扑结构恒定(P和F共同表示了人脸的形状模型，即Face＝{P，F}。具体地，P＝{Mpti＝(XptiYptiZpti)T|pti＝1，…，nVertex}是人脸上的三维结点列表，而F＝{facetfj＝(pt1，fjpt2，fjpt3，fj)T|fj＝1，…，nMesh}是人脸上的三角形面片列表，其中每个facetfj＝(pt1，fj，pt2，fjpt3，fj)T表示由P中三个点Mpt1,fj,Mpt2,fj,Mpt3,fj]]>组成的一个三角形。ptpti，fj，pti＝1，…，3，fj＝1，…，nMesh，是P中各点的序号，满足1≤ptpti，fj≤nVertex。nVertex和nMesh分别是三维点和三角面片的总数目。总的说来，P表示模型形状，F则表示模型的拓扑结构)。沿这条技术路线可以保持模型的拓扑结构，因而能方便地实现人脸动画；同时，该通用模型也表示了我们对人脸这种非常相似的物体类的一般知识，能为具体的处理过程提供附加约束。
下面我们将详细说明建模步骤，并给出一个具体实施实例。
人脸建模的总体流程如图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两个图象序列中分别进行跟踪，接着同样通过立体重建得到。有了不同时刻的三维点后我们就可以估计姿态。从未知量的个数来考虑，只需要三对两个时刻的对应三维点就可完成，但要求重建出的三维点是可靠的，而实际实验中很难保证所有的三维点都可靠。为了评价三维点的可靠性，我们定义并不断更新数组{bReliablek}，这是一个布尔数组，取值为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)表示图象上的任一点；该点处图象的梯度向量为g→=▿Im,t(p)=(∂Im,t∂x(p)∂Im,t∂y(p))T]]>，计算2*2的矩阵G=∫Wg→g→Tdp]]>，其中W是以p为中心的一个矩形窗口。我们求出G的最小奇异值sing_min(利用已有的matlab工具包)，并选择那些sing_min很大的p点作为检测出的特征点。从物理意义上看，G实际是对图象上局部灰度变化丰富程度的表示：如果局部的灰度恒定，这时梯度向量 为0，G为0，所有的奇异值都是0；如果局部的灰度变化有固定的方向(例如图象上比较强的单一边缘处)，梯度向量 基本都垂直于该方向，这时G最小的奇异值接近于0；只有在沿各个方向的灰度变化都比较显著的局部(例如两条边的交点)才满足G的最小奇异值很大。因此，我们实际是选择那些局部灰度变化丰富的点作为特征点，这样的点在跟踪和匹配时可靠程度较高。
计算机实现时我们需要对上面的公式作离散化处理，具体的计算公式如下：g→(i,j)=gxgyT=[Im,t(i+1,j)-Im,t(i-1,j)2Im,t(i,j+1)-Im,t(i,j-1)2]T---(1.1.1)]]>G(i,j)=Σii=-blockxblockxΣjj=-blockyblockyg→(i+ii,j+jj)g→T(i+ii,j+jj)---(1.1.2)]]>这里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)后，令F=fL0uOL0fLvOL001-T0-Tc,zTc,yTc,z0-Tc,x-Tc,yTc,x0RcfR0uOR0fRvOR001-1]]>(该矩阵被称为基础矩阵)，其中fL，fR，uOL，vOL，uOR，vOR是摄像机的内参数，都是已知量；其矩阵表示是AL=fL0uOL0fLvOL001,]]>AR=fR0uOR0fRvOR001.]]>则对于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)：pL,t(t)1TFpR,t(k)1=0---(1.1.4)]]>(1.1.4)式的物理意义是pR，t(k)位于由pL，t(k)决定的一条直线上，这样为寻找该点的搜索空间从整个图象降为该直线，即从二维降到了一维。在这个一维空间里，选择使一个局部窗口Wm内总灰度差异最小的点作为匹配结果，即要使下面的(1.1.5)式最小：Diff=∫Wm[IR,t(pR,t(k))-IL,t(pL,t(k))]2dp---(1.1.5)]]>离散化后，实际计算时的公式是：Diff(i,j)=Σii=-blockxblockxΣjj=-blockyblocky(IL,t(xL,t(k)+ii,yL,t(k)+jj)-IR,t(i+ii,j+jj))2---(1.1.6)]]>其中(i，j)是直线 F上的任一点。
图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)}。立体重建的过程对每一对匹配点是独立的，原理如下：首先，我们有摄像机的透视投影原理，即下面的已知公式，smAm-1pm,t1=Mm,t]]>，其中m为L或R，表示两个摄像机中的任一个；Mm，t为t时刻三维点在m摄像机坐标系中的坐标；sm是待求的一个比例因子。我们省略了(k)来表示任一点。
再结合两摄像机间的配置关系(1.1.3)可求解出三维点，ML，t的表达式为ML,t=12(sLAL-1pL,t1+sRRcAR-1pR,t1+Tc)]]>其中SLSR=(BTB)-1BTTc]]>，而B=[AL-1pL,t1-RcAR-1pR,t1].]]>而MR，t则可按下面的公式求出MR,t=Rc-1(ML,t-Tc)]]>如果不特别说明，在本文中提到的立体重建出的三维点都是指在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上跟踪该点。基于局部纹理信息的匹配，令e→=∫W(Im,t+1(p)-Im,t(p))g→dp,]]>d→=G-1e→]]>(G的计算如前所述)，则 即为图象Im，t+1上对p的跟踪结果。
实际计算时同样需要对 的求解作离散化处理：e→(i,j)=Σii=-blockxblockxΣjj=-blockyblockyg→(i+ii,j+jj)(Im,t+1(i+ii,j+jj)-Im,t(i+ii,j+jj))---(1.1.7)]]>跟踪得到的特征点在图2中表示为(pm，t+1(k)}。
1.1.4姿态参数初始化对于跟踪的结果用和匹配时相同的立体重建方法可获得t+l时刻的三维点序列，这样我们就获得了这两个时刻(t和t+1)的一组三维对应点{M，(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的向量)和绕该轴的旋转角度θ。定义单位向量dirim=Mim-M0||Mim-M0||,]]>im＝1，2和dirim′=Mim′-M0′||Mim′-M0′||,]]>im＝1，2。则r^=(dir1′-dir1)×(dir2′-dir2)||(dir1′-dir1)×(dir2′-dir2)||]]>在得到旋转轴 后，θ可以按下式进行求解：θ=arccos(r^×dir1)·(r^×dir1′)||(r^×dir1)·(r^×dir1′)||]]>得到了 和θ之后，旋转矩阵Rt可由下式给出：Rt=I+sin||r||||r||[r]x+1-cos||r||||r||2[r]x2]]>其中r=θ*r^]]>，[r]x是根据3*1向量r＝[rxryrz]T定义的反对称矩阵，即[r]x=rxryrzx=0-rzryrz0-rx-ryrx0,]]>而I=100010001]]>是单位矩阵。估计出旋转分量Rt之后，平移分量Tt可由下式得到：Tt=13Σim=02(Mim′-RtMim)]]>上面的算法虽然在理论上能得到准确的结果，但是数值稳定性很差，即在其输入{(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}，其中n_m=argminnmedsn{ϵn,s|s=1,…,Kr(t)}---(1.1.9)]]>其中 表示针对下标s取集合{εn，s|s＝1，…，Kr(t)}的中值，它对不同的n构成一个新的数组{medsn|n=1,…,Num_subset}]]>；而 则是在这个中值数组中选择取最小值的那个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的均值ϵ‾=1K_TrueΣϵs∈Eϵs]]>和标准差σ=1K_True-1Σϵs∈E(ϵs--ϵ-)2;]]>(3)在Errt中寻找满足‖εs- ε‖＞3*σ的点，设数目为InValid_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(r)}，使得式(1.1.10)中的f_obj达到其最小值。f_obj=ΣmΣτ=0t-tStartΣpp=1Kr(t)[(app,mτ-xpp,mτ)2+(bpp,mτ-ypp,mτ)2]---(1.1.10)]]>其中RatStart+τ，TatStart+τ是从tStart时刻到tStart+τ时刻的刚体变换系数；{MtStart(kpp)}，pp＝1，…，Kr(t)，是在这Num个时刻中都可靠的三维点数组，包含Kr(t)个点，pp表示三维点的标号。m为L或R，表示某一个摄像机。 是二维检测或跟踪的结果，在m号摄像机tStart+τ时刻的图象中，即app,mτbpp,mτ1=Am-1pm,tStart+τ(kpp)1---(1.1.11)]]>而 则是pp号三维点在tStart+τ时刻，投影到m摄像机上的结果，即Rc,m=R1c,mTR2c,mTR3c,mT]]>和Tc,m=T1c,mT2c,mT3c,m]]>表示从L摄像机坐标系到m摄像机坐标系的刚体变换。即Rc,L=100010001,]]>Tc,L=000,]]>Rc,R=Rc-1,Tc,R=-Rc-1Tc,R]]>。另外有RatStart=100010001,]]>TatStart=000.]]>从物理意义上看，式(1.1.10)表示的是所有图象中二维重投影的总误差。优化该目标能得到姿态参数的精确值。其中参数{RatStart+τ，TatStart+τ，τ＝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，，τ＝2，…，t-tStart上面公式的物理意义是刚体运动的复合，例如：{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)完成：Minit=sx000sy000szMg---(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，l，Tg，l，即从generic坐标系到第1帧时L摄像机坐标系的刚体变换，即要求解(1.2.2)中的参数M1＝Rg，lMinit+Tg，l(1.2.2)其中Ml和Minit分别是时刻1时人脸和Faceinit这两个模型间的任一组对应点。
为求解Rg，l，Tg，l，我们选择这两个模型在LE(左眼)、RE(右眼)和Mid_M(嘴中心)这三个位置处的对应点，按前面1.1.4节中姿态初始化的方法得到结果。
获得了Rg，l，Tg，l后，系统根据刚体运动的复合，按下面的公式计算其他各时刻的姿态参数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中的坐标 可由下面的已知公式求出：Mptnew=ΣjRBF=1SCjRBFφ(||Mptstart-StajRBF||)]]>其中 是pt点在Facestart中的坐标，而{C1，…，CS}可按下面的公式求出C1T...CST=φ(||Sta1-Sta1||)...φ(||Sta1-StaS||).........φ(||StaS-Sta1||)...φ(||StaS-StaS||)-1New1T...NewST]]>
该方法实际是根据Facestart中各点到一些标定点(即Subset中的点)的距离关系来确定该点在Facenew中的位置。φ是已知的径向基函数(定义域和值域都为非负实数)，它是自变量r(表示距离)的递减函数，即距离远的标定点对该点的影响较小。
这一步在利用上面的径向基函数算法时，我们有Facestart＝Faceinit，Facenew＝Face1。而Subset中的点包括两类：一类是三维显著特征点，共4个，由二维显著特征点(指人脸上的左眼内角LE、右眼内角RE、左外嘴唇LM、右外嘴唇RM四个点在图象平面上的投影)通过前面的立体重建算法得到；另一类是侧面特征点(指那些在人脸侧面投影视图上易于辨认的点，如下頦、嘴、鼻尖、鼻梁上方、前额等。典型的人脸侧面投影以及我们定义的一些侧面特征点在图6中用小的圆圈表示，它的一个实例参见图16)，实际使用的有24个，它们的三维位置按下面的步骤从手工标出的二维点中获得。
设手工标出的侧面特征点为{mL,sidekk=xL,sidekkyL,sidekkT},]]>kk＝1，…，NumSF，其中NumSF＝24为侧面特征点的数目，side是前面已经得到的侧面帧标号。它们在generic坐标系中的三维坐标为{Msgkk=0YgkkZgkkT},]]>kk＝1，…，NumSF，则对任一kk，我们有其中Rg,side=Rg,side,1TRg,side,2TRg,side,3T,]]>Tg,side=Tg,side,1Tg,side,2Tg,side,3]]>是姿态参数的分量表示。(1.3.1)中有两个线性方程，因而可直接用线性方程组的求解方法得到 中的两个未知数 和 该方法实际是利用了侧面特征点的特殊位置(即在generic坐标系中坐标的X分量为0)，因而可直接从单幅图象中求出三维点。
1.4形状修正B这一步同样利用前面的径向基函数算法，Facestart＝Face1，Facenew＝Face2，Subset中除了包括前面“形状修正A”里已经确定了的点之外，还有一些新添加的点，即根据手工标出的正面人脸轮廓(指正面姿态下图像上的人脸轮廓，用闭合的多边形轮廓表示，一个例子如图7所示。该图将多边形轮廓叠加到了人脸图象上，看得更清楚)恢复出来的三维点。得到这些三维点的基本方法是将当前人脸模型向图象平面作透视投影，并计算三维模型透视投影的二维轮廓点，最后依赖于这些轮廓点和图象上的人脸轮廓的匹配来计算其三维坐标。具体的步骤如下所述：首先需要计算出Face1在t＝1时刻(是手工选择的正面视图)向图象平面投影的轮廓Cont1={pt1iC1|iC1=1,…,nNum1}]]>，其中 是P(即模型的三维结点列表)中点的标号，显然有1≤p1iC1≤nVertex]]>；而nNum1是Cont1中的顶点数目。Cont1中的点就是这一步形状修正中新添加到Subset中的点，它们是通过判断投影线和模型的交点来选择的，即在某一轮廓点上，投影线与模型不能有其他的交点，详细的计算方法如下所述：记三维模型为Faceproj(这里就是Face1)，透视投影中心为Center=-Rg,1-1Tg,1]]>即L摄像机光心在generic坐标系中的坐标，则我们对模型中的所有顶点计算过该点的投影线和模型的交点，据此判断该顶点是否在投影轮廓中。具体地，Mnow为Faceproj中任一点，Planejp为Faceproj中任一面片构成的平面，计算过Center和Mnow的直线(即投影线)与Planejp的交点为Cnow；如果对某一个满足1≤ip≤nMesh的ip(如前所述，nMesh是Faceproj中的面片总数)求出的Cnow在线段CMnow的内部，则Mnow不在投影轮廓中；否则就应属于投影轮廓。
接着要计算Cont1中的任一点 在Face2的三维坐标。设图象上的人脸轮廓是Cont_img1={xCI1iCimgyCI1iCimg|iCimg=1,...,nNumCI1}]]>，它是由二维点列表构成的多边形，而nNumCI1是点的数目。则我们按下面的方法确定 在Face2的三维坐标：设pti=pt1iC1]]>点在Faceg中的法向为vn，在Face1中的坐标是M1pti=X1ptiY1ptiZ1ptiT]]>，而经过该点的投影方向为vpn=M1pti+Rg,1-1Tg,1]]>(这三者都是已知量)，则pti点在Face2中的坐标满足M2pti=M1pti+tlinev---(1.4.1)]]>其中v＝(vxvyvz)T＝vpn×(vn×vpn)            (1.4.2)参数tline可根据二维直线X1pti+tlinevxZ1pti+tlinevzY1pti+tlinevyZ1pti+tlinevzT]]>与闭合多边形Cont_img1的交点求出。具体步骤如下：设Cont_tmg1中的任一线段可表示为seg＝(x0+sdxy0+sdy)T，0≤s≤1，其中s为线段参数，其他都是已知量。
我们先按下面的公式计算sitline=dx-(X1pti-vx*Z1pti/vz)dy-(Y1pti-vy*Z1pti/vz)-1x0-vx/vzy0-vy/vz]]>对所有的seg求出这两个数值si和tline。对于很多seg上面的公式无解，即公式里的矩阵不可逆。有解且得到的si满足0≤si≤1的seg有一个或两个，取与点X1ptiZ1ptiY1ptiZ1ptiT]]>距离近的那一个seg，它得到的tline即为所求。
这里v是 点的变形方向，对它的计算利用了通用人脸的局部形状。实际上，该方向是Face1上该点处法向垂直于投影线的分量，这样一方面能尽量保持局部曲率；另一方面能尽量减小点的变形距离。
1.5形状修正C这一步同样利用前面的径向基函数算法，Facestart＝Face2，Facenew＝Face3，Subset中除了包括前面“形状修正A、B”里已经确定了的点之外，还有一些新添加的点，即根据手工标出的正面人脸特征轮廓点(指正面姿态下手工标出的眼睛、鼻孔和嘴的轮廓点，是一些离散点，如图8所示，同样是叠加显示)恢复出来的三维点。具体的恢复步骤如下所述：设某特征轮廓点的二维坐标是xCF1iCFyCF1iCF]]>，该点在Face2中的坐标为M2iCF=X2iCFY2iCFZ2iCFT]]>，我们假设该点的Z坐标已经计算正确，即它在Face3中的Z坐标等于 ，则它在Face3中的坐标为M3iCF=xCF1iCF*Z2iCFyCF1iCF*Z2iCFZ2iCFT---(1.5.1)]]>即这些特征轮廓点在正面姿态下的投影等于实际图象上的特征点。
1.6形状修正D这一步同样利用前面的径向基函数算法，Facestart＝Face3，Facenew＝Face4，Subset中除了包括前面“形状修正A、B、C”里已经确定了的点之外，还有一些新添加的点，即从手工标出的中间视图1(该视图的选择由系统自动判断，实际是使Rg，t的旋转角最接近25度的时刻t。所谓“中间视图”是指介于正面和侧面之间的中间投影姿态)人脸轮廓(同样用闭合多边形表示，典型的例子是图9)恢复出来的三维点，具体的恢复步骤如下所述：设中间视图1对应于时刻int1，与前面1.4节对正面图人脸轮廓的处理类似，也要先算出Face3在int1时的投影轮廓Contint1={ptint1iC2|iC2=1,…,nNum_int1}]]>，其中 是P(即人脸模型的三维结点列表)中点的标号，而nNum_int1是投影轮廓中的顶点数目。投影轮廓的计算方法如前所述(1.4节)。
Contint1中的点在Face4中的坐标按下面的与1.4节类似的步骤确定：对于任一ptli=ptint1iC2,]]>定义v2=v2xv2yv2zT=M3ptli-Center]]>，其中M3pt1i=X3pt1iY3pt1iZ3pt1iT]]>是该点在Face3中的三维坐标，而Center=-Rg,1-1Tg,1]]>如前所述是L摄像机光心在generic坐标系中的坐标。则该点在Face4中的坐标 满足M4ptli=M3ptli+tline2v2---(1.6.1)]]>其中参数tline2根据二维直线X3pt1i+tline2v2xZ3pt1i+tline2v2zY3pt1i+tline2v2yZ3pt1i+tline2v2zT]]>与闭合多边形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中的任一三维点Mspt=XsptYsptZsptT]]>作公式(1.8.1)所示的柱面映射：映射的结果是一个二维平面，其上的点为{Cynpt}={lgptltpt}]]>。这里lgpt是柱面上的经度位置，用弧度值表示；ltpt则是柱面旋转轴向的坐标位置。
1.8.2正面纹理图将所有 在正面姿态参数Rg，l，Tg，l下按公式(1.8.2)进行投影：其中Rg,1=R1g,1TR2g,2TR3g,3T,]]>Tg,1=T1g,1T2g,1T3g,1]]>是旋转矩阵和平移向量的分量表示。
这时p→1pt=x1pty1ptT]]>是Faces中三维点 在正面视图IL，1上的二维投影。现在我们得到一个IL，1平面上点 和Cyn平面上点Cynpt＝(lgptltpt)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平面上的对应三角形，即piT1=x1ptiT1y1ptiT1T]]>piT1′=lgptiT1ltptiT1T,]]>iT1＝1，2，3。对三角形mcyn内的点[lg lt]按下面的公式(1.8.4)计算[xy]，并令ITexture1(lg,lt)=IL,1(x,y)---(1.8.3)]]>即将正面纹理图映射到柱面展开图上。xy1=a1a2a3a4a5a6001-1lglt1---(1.8.4)]]>其中的六个未知量aiA，iA＝1，…，6按公式(1.8.5)求出。正面纹理映射的结果记为Texture1。
1.8.3侧面纹理图与正面图的处理非常类似，将所有 在侧面姿态参数Rg，side，Tg，side下按公式(1.8.6)进行投影：这里同样有分量表示Rg,side=R1g,sideTR2g,sideTR3g,sideT,]]>Tg,side=T1g,sideT2g,sideT3g,side.]]>接着对Faces中的任一面片facet＝(pt1pt2pt3)T，设msproj＝Δ(ps1，ps2，ps3)，其中psiT2=xsideptiT2ysideptiT2T]]>，iT2＝1，2，3。对三角形mcyn内的点[lg lt]按后面的公式(1.8.8)计算[xs ys]，并令ITextureside(lg,lt)=IL,side(xs,ys)---(1.8.7)]]>即将侧面纹理图映射到柱面展开图上。xsys1=as1as2as3as4as5as6001-1lglt1---(1.8.8)]]>其中的六个未知量asiAS，iAS＝1，…，6按公式(1.8.9)求出。侧面纹理映射的结果记为Textureside。
1.8.4侧面纹理图的反射由于人脸只向一侧旋转，所以只能得到单侧的侧面纹理，另一侧的纹理可通过对Textureside作反射处理得到。具体地说，由于Faces的拓扑结构是完全对称的，设有直接获得纹理的那一侧脸上的任一三角面片mcyn＝Δ(p1，p2，p3)和它在另一侧脸上的对称面片mcyn′＝Δ(p1′，p2′，p3′)，其中piT3=lgptiT3ltptiT3T,]]>piT3′=lg′ptiT3lt′ptiT3T,]]>iT3＝1，2，3，则对mcyn中的任一点p＝(lg lt)T按下面的公式(1.8.9)计算[lg′lt′]，并令ITextureside(lg′,lt′)=ITextureside(lg,lt)---(1.8.8)]]>即从直接获得纹理的一侧反射生成另一侧。lg′lt′1=rs1rs2rs3rs4rs5rs6001-1lglt1---(1.8.9)]]>其中的六个未知量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为IntelPIII933，配有256M的Rambus内存，安装了Microsoft Windows 2000 Professional操作系统。为进行数据采集，我们使用了两个相同的SANYO VCC-5972P型摄像机，以及两块相同的MatroxMeteorII采集卡来采集两路同步视频，图象分辨率为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同样表示时间，表示当前跟踪着的特征被检测出来的时刻；●计算基础矩阵F=AL-T0-Tc,zTc,yTc,z0-Tc,x-Tc,yTc,x0RcAR-1]]>(其中各已知量的意义参见1.1.2节中的解释)。实际数值是Rc=0.999952-0.005045+0.0084570.0065790.981257-0.192594-0.0073270.1926400.981242,]]>Tc=0.0124130.6112400.088092]]>AL=fL0uOL0fLvOL001,]]>AR=fR0uOR0fRvOR001,]]>其中fL＝1330.3，fR＝1330.6，uOL＝384，vOL＝288，uOR＝384，vOR＝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；●φ：即形状修正时已知的径向基函数，我们选择φ=(1-r)+6(6+36r+82r2+72r3+30r4+5r5)]]>，其中 表示函数取值为负时就置其为0。
第二步：在图象IL，t上计算梯度向量g→L,t(i,j)=gxgyT=[IL,t(i+1,j)-IL,t(i-1,j)2IL,t(i,j+1)-IL,t(i,j-1)2]T]]>以及G(i,j)=Σii=-blockxblockxΣjj=-blockyblockyg→L,t(i+ii,j+jj)g→L,tT(i+ii,j+jj),]]>其中1≤i≤width，1≤j≤height表示图象上任一点，blockx，blocky是与窗口尺寸相关的常量。如果计算时用到了超出图象边界的点，直接令GL,t(i,j)=0000]]>。接着利用已有的数学工具软件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)局域图象上每一点的灰度值为48638197102907257781091351371219657981461751771591306110215318017716213359951341501461301005480107121113957248576169665650]]>而7*7的{gx}和{gy}矩阵为按前述公式计算出的2*2的矩阵G=233112665.52665.530700]]>，其两个奇异值为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)，计算其对应的外极线 F，并计算Diff(i,j)=Σii=-blockxblockxΣjj=-blockyblocky(IL,t(xL,t(k)+ii,yL,t(k)+jj)-IR,t(i+ii,j+jj))2]]>其中(i，j)是直线l(其向量表示是l＝[l1l2l3]，即直线的方程是l1i+l2j+l3＝0)上的任一点，即满足||l1i+l2j+l3||l12+l22<1.]]>(imin,jmin)=argmin(i,j)Diff(i,j).]]>则pL，t(k)的匹配点即为pR，t(k)＝(iminjmin)T。对所有检测出的特征点进行相同的处理，我们得到IR，t上的匹配点序列，并记为{pR，t(k)}，1≤k≤K(t)。
第五步：根据任一对匹配pL，t(k)和pR，T(k)，计算B=[AL-1pL,t1-RcAR-1pR,t1],]]>sLsR=(BTB)-1BTTc]]>，并最终得到三维点Mt(k)=12(sLAL-1pL,t(k)1+sRRcAR-1pR,t(k)1+Tc)]]>。对所有成对匹配进行相同处理，我们得到t时刻的三维特征点序列{Mt(k)|k＝1，…，K(t)}。同时初始化布尔变量序列{bReliablek}，k＝1，…，K(t)全为Yes，表示所有点都是可靠的。
第六步：在任一pL，t(k)＝(i，j)处计算e→L,t(i,j)=Σii=-blockxblockxΣjj=-blockyblockyg→L,t(i+ii,j+jj)(IL,t+1(i+ii,j+jj)-IL,t(i+ii,j+jj))]]>，以及d→L,t(i,j)=GL,t(i,j)-1e→L,t(i,j)]]>。则该点在L视频序列中t+1时刻的对应点为pL,t+1(k)=pL,t(k)-d→L,t(i,j)]]>。对所有的特征点进行相同的处理，我们得到IL，t+1上的跟踪结果，并记为{pL，t+1(k)}，1≤k≤K(k)。
第七步：用类似的方法处理R视频序列。具体地，先按前述公式在R，t图上计算和{GR，t(iR，jR)}，接着对任一pR，t(k)＝(iR，jR)计算e→R,t(iR,jR)=Σii=-blockxblockxΣjj=-blockyblockyg→R,t(iR+ii,jR+jj)(IR,t+1(iR+ii,jR+jj)-IR,t(iR+ii,jR+jj))]]>，最终有d→′R,t(iR,jR)=GR,t(iR,jR)-1e→R,t(iR,jR)]]>，则该点在R视频序列中t+1时刻的对应点为pR,t+1(k)=pR,t(k)-d→′R,t(iR,jR)]]>。对所有特征点进行相同的处理，我们得到IR，t+1上的跟踪结果，并记为{pR，t+1(k)}，1≤k≤K(t)。
第八步：用与第五步类似的方法得到t+1时刻的三维点序列。具体地，根据任一对匹配pL，t+1(k)和pR，t+1(k)，计算B=[AL-1pL,t+1(k)1-RcAR-1pR,t+1(k)1].]]>sLsR=(BTB)-1BTTc,]]>并最终得到三维点Mt+1(k)=12(sLAL-1pL,t+1(k)1+sRRcAR-1pR,t+1(k)1+Tc)]]>。对所有成对的跟踪结果进行相同处理，我们得到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)。根据{Mt(kN1,n),Mt(kN2,n),Mt(kN3,n)}]]>和{Mt+1(kN1,n),Mt+1(kN2,n),Mt+1(kN3,n)}]]>的对应关系，按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}，其中n_m=argminnmedsn{ϵn,s|s=1,…,Kr(t)}]]>第十步：计算序列Errt＝{εs＝‖Mt+1(ks)-RtMt(ks)-Tt‖2|s＝1，…，Kr(t)}，并更新{bReliableks},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和6I(即各个K(t0))，在t0保持不变的所有时刻(即两次检测之间的特征跟踪)里一直可靠的特征点(即序列{bReliablek}中取值为Yes点的数目)分别有12、18、36、14个。
第十四步：在IL，1和IR，1上手工选择左内眼角、右内眼角、左嘴角和右嘴角这四个二维显著特征点，实际结果如图15所示。根据这四个位置处的两幅图间的匹配，可按前述的立体重建算法得到四个三维显著特征点，即左内眼角 ，右内眼角 ，左嘴角 ，右嘴角M 。然后按公式(1.2.1)将原始通用模型Faceg变形为Faceinit。其中三个比例因子sx，sy，sz的计算方法如下：令sx=dist_3D(M1LE,M1RE)dist_3D(MgLE,MgRE)]]>，即建模对象和Faceg中Seg_E长度(两内眼角的距离)的比值；令sy=dist_3D(M1LE+M1RE,M1LM+M1RM)dist_3D(MgLE+MgRE,MgLM+MgRM),]]>即Seg_EM长度(两眼中心到嘴中心的距离)的比值；令sz＝sy。其中MgLE,MgRE,MgLM,MgRM]]>是在Faceg上这四个主要特征点的三维坐标。
第十五步：接着我们求解generic坐标系到L摄像机坐标系的刚体变换Rg，l，Tg，l。我们选择{MinitLE,MinitRE,(MinitLM+MinitRM)/2}]]>和{M1LE,M1RE,(M1LM+M1RM)/2}]]>这三个三维对应点，按前述的姿态参数初始化一节的方法求解。
获得了Rg，l，Tg，l后，系统按下面的公式计算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所示，记为{mL,sidekk=xL,sidekkyL,sidekkT},]]>kk＝1，…，NumSF，其中NumSF为侧面特征点的数目，实际为24。
第十七步：按前述方法求解式(1.3.1)，得到三维侧面特征点{Msgkk=0YgkkZgkkT},]]>kk＝1，…，NumSF。
第十八步：按前述的1.3节中的径向基函数算法把Faceinit变形为Face1。Subset中的点包括两类：一类是三维显著特征点，共4个，即{MinitLE,MinitRE,MinitLM,MinitRM}]]>；另一类是侧面特征点，即十七步求出的 kk＝1，…，NumSF。
第十九步：按前述1.4节中计算投影轮廓的方法得到Face1在时刻1时投影到图象平面的轮廓点Cont1={pt1iC1|iC1=1,…,nNum1}]]>，Face1就是那里的Faceproj。
第二十步：手工标出IL，1上的人脸轮廓，如图7所示。这是一个闭合的多边形，记为Cont_img1={xCI1iCimgyCI1iCimg|iCimg=1,...,nNumCI1}]]>，是连续的二维顶点列表，nNumCI1是二维顶点的数目。
第二十一步：按照前述的公式(1.4.1)确定Cont1中点在Face2中的三维坐标{M2pt1iC1|iC1=1,…,nNum1}.]]>第二十二步：按前述1.3节中的径向基函数算法把Face1变形为Face2。Subset中新添加的点为第二十一步求出的{M2pt1iC1|iC1=1,…,nNum1}.]]>第二十三步：手工标出IL，1图上的正面人脸特征轮廓点(实验结果为图8)，结果为xCF1iCFyCF1iCF|iCF=1,...,nNumCF}]]>，其中nNumCF是点的总数目。
第二十四步：按前述的公式(1.5.1)求解第二十三步标出的特征轮廓点的三维坐标 第二十五步：按前述1.3节中的径向基函数算法把Face2变形为Face3。Subset中新添加的点为第二十四步求出的 第二十六步：系统自动判断对应于中间视图1的帧int1，即选择int1为使Rg，t的旋转角最接近25度的时刻t。接着在IL，int1图上的手工标出人脸轮廓，如图9所示。同样用闭合多边形表示，记为Cont_imgint1={xCIint1iCimg1yCIint1iCimg1|iCimg1=1,...,nNumCIint1}]]>，nNumCIint1是顶点数。
第二十七步：按前述1.4节中计算投影轮廓的方法得到Face3在时刻int1时投影到图象平面的轮廓点Contint1={ptintliC2|iC2=1,…,nNum_int1}]]>，Face3就是那里的Faceproj。
第二十八步：根据Cont_imgint1确定Contint1中点在Face4中的三维坐标。具体地，对于任一pt1i=ptint1iC2]]>，定义v2=v2xv2yv2zT=M3pt1i+Rg,l-1Tg,l]]>，其中M3pt1i=X3pt1iY3pt1iZ3pt1iT]]>是该点在Face3中的三维坐标。则该点在Face4中的坐标 满足M4pt1i=M3pt1i+tline2v2]]>。其中参数tline2根据二维直线X3pt1i+tline2v2xZ3pt1i+tline2v2zY3pt1i+tline2v2yZ3pt1i+tline2v2zT]]>与闭合多边形Cont_imgint1的交点求出，方法如1.4节所述。对Contint1中所有点按上面的步骤计算，得到{M4pt1iC2|iC2=1,…,nNum_int1}.]]>第二十九步：按前述1.3节中的径向基函数算法把Face3变形为Face4。Subset中新添加的点为第二十八步求出的{M4pt1iC2|iC2=1,…,nNum_int1}.]]>第三十步：系统自动判断对应于中间视图2的帧int2，即选择int2为使Rg，t的旋转角最接近40度的时刻t。接着在IL，int2图上的手工标出人脸轮廓，如图10所示。同样用闭合多边形表示，记为Cont_imgint2={xCIint2iCimg2yCIint2iCimg2|iCimg2=1,...,nNumCIint2}]]>，nNumCIint2是顶点数。
第三十一步：按前述1.4节中计算投影轮廓的方法得到Face4在时刻int2时投影到图象平面的轮廓点Contint2={ptint2iC3|iC3=1,…,nNum_int2}]]>，Face4就是那里的Faceproj。第三十二步：根据Cont_imgint2确定Contint2中点在Face5中的三维坐标。具体地，对于任一pt2i=ptint2iC3]]>，定义v3=v3xv3yv3zT=M4pt2i+Rg,1-1Tg,1]]>，其中M4pt2i=[X4pt2iY4pt2iZ4pt2i]T]]>是该点在Face4中的三维坐标。则该点在Face5中的坐标
满足M5pt2i=M4pt2i+tline3v3]]>。其中参数tline3根据二维直线X4pt2i+tline3v3xZ4pt2i+tline3v3zY4-pt2i+tline3v3yZ4pt2i+tline3v3zT]]>与闭合多边形Cont_imgint2的交点求出，方法如1.4节所述。对Contint2中所有点按上面的步骤计算，得到{M5ptint2iC3|iC3=1,…,nNum_int2}.]]>第三十三步：按前述1.3节中的径向基函数算法把Face4变形为Face5。Subset中新添加的点为第三十二步求出的{M5ptint2iC3|iC3=1,…,nNum_int2}]]>。Face5即为最终的形状模型Faces。
这样我们一共进行了五次形状修正，对应于1.3至1.7节的形状修正A-E，实际上各步使用的Subset中的点数目如下面的表1所示。
表1形状修正中每一步新加入的Subset点的位置及数目第三十四步：对Faces中的任一三维点Mspt=XsptYsptZsptT]]>按公式
作柱面映射，结果为{Cynpt}＝{(lgptltpt)}。
第三十五步：将所有 在正面姿态参数Rg，l，Tg，l下按公式 进行投影，得到{p→1pt}={x1pty1ptT}]]>。然后按公式(1.8.3)得到正面纹理图Texture1，实际例子如图17所示。
第三十五步：与第三十四步类似，获得Faces的侧面纹理图Textureside。具体地，将所有 在侧面姿态参数Rg，side，Tg，side下按公式 进行投影，得到{p→sidept}={xsideptysideptT}]]>。然后按公式(1.8.7)得到侧面纹理图Textureside，实际例子如图18所示。
第三十六步：在柱面图{(lg，lt)}上进行纹理反射。实际实验中直接获得纹理的那一侧脸是左侧，因此这一步根据左侧的纹理信息按公式(1.8.8)生成右侧脸。这时完整的Textureside如图19所示。
第三十七步：我们按公式 融合Texture1和Textureside来得到最终的纹理图Textures，实际结果如图20所示。
我们最终的、映射了Textures的建模结果Faces如图21所示。