基于金字塔的LK特征点跟踪算法

采用金字塔执行L-K特征点跟踪算法的描述西北工业大学覃勋辉译1、问题描述 假设有两幅大小相同的2D灰度图像I和J，它们在点[]Txy=x处的灰度值分别为()(,)IIxy=x和()(,)JJxy=x，通常称I为第一幅图像，J为第二幅图像。实际中，图像I和J均是离散的函数，图像左上顶点的坐标值为[]00T。令xn和yn分别是图像的宽度和高度，则图片右下顶点的坐标值为11Txynn⎡⎤−−⎣⎦。 若已知第一幅图像I中有特征点[]Txyuu=u，在第二幅图像J中存在特征点[]Txyuu=+=vud使得()Iu和()Jv近似，则特征点跟踪的任务是计算J中的对应特征点v。 向量[]Txydd=d叫做在点x处的图像速率，也称为在点x处的光流。由于摄像机针孔成像模型带来透视关系使得完全的对应点并不存在，所以有必要在特征点的2D图像邻域中定义相似性这一概念方便描述特征点的对应关系。相似性定义如下：令xw和yw为两个正整数，存在光流d使得残差ε最小，此时图像特征点间相似度最高。可表示为： ()()2(),(,)(,)yyxxxxyyuwuwxyxyxuwyuwddIxyJxdydεε+−=−=−==−++∑∑d        （1） 由上式可知，相似性函数可由图像特征点的(21)(21)xxww+∗+邻域求得，这个邻域范围称为综合窗口，其中xw和yw通常取2~7pixs。 2、描述跟踪算法 精确性和鲁棒性是判断特征点跟踪算法优劣的关键因素。跟踪算法的局部亚像素定位能力决定了算法的精确性。很明显小的综合窗口（对应小的xw和yw值）有助于保持图像中的局部细节，当图像块间存在较大速率时（形变严重）常用小的综合窗口。 鲁棒性是指算法对光照变化、图像区域移动时稳健的跟踪能力。特别是为了使算法能够处理大的图像区域移动，需要设置一个大的综合窗口。的确，如果仅考虑公式（1），算法已经默认为xxdw≤和yydw≤（除非有先验的匹配信息可用否则只能选择大的综合窗口）。因此，难以选择一个适当的综合窗口使其同时满足算法的精确性和鲁棒性要求。为了解决这个问题，本文引入迭代的高斯金字塔实现经典的L‐K特征点跟踪算法，使得L‐K光流计算达到理想的局部跟踪定位精度。 2.1 建立图像金字塔 下面建立大小为xynn∗图像I的金字塔。令0II=是第0层的图像，它是金字塔图像中分辨率最高的图像，图像的宽度和高度分别定义为0xxnn=和0yynn=。以一种递归的方式建立金字塔：从0I中计算1I，从1I中计算2I，···。令1,2,...L=代表金字塔的层数，1LI−是第1L−层的图像，1Lxn−和1Lyn−分别是图像1LI−的宽度和高度。图像LI可按如下方式由1LI−求得： ()112,24LLIIxy−=+ ()()()()()1111121,221,22,212,218LLLLIxyIxyIxyIxy−−−−−+++−+++ ()()()()()1111121,2121,2121,2121,2116LLLLIxyIxyIxyIxy−−−−−−++−+−++++ （2） 为了满足公式（2），必须在图像1LI−（101Lxxn−≤≤−，101Lyyn−≤≤−）的周围定义一些虚拟的像素（实际并不存在）： 111111111111111111(1,)(0,)(,1)(,0)(,)(1,)(,)(,1)(,)(1,1)LLLLLLLLxxLLLLyyLLLLLLxyxyIyIyIxIxInyInyIxnIxnInnInn−−−−−−−−−−−−−−−−−−−=−==−=−=−− 此时，由公式（2）定义的图像LI坐标(),xy需满足1021Lxxn−≤≤−和1021Lyyn−≤≤−。因此，图像LI的宽度Lxn和高度Lyn的最大整数满足下列两个条件： 112LLxxnn−+≤                                （3） 112LyLynn−+≤                                （4） 由方程（2）（3）（4）可以从下到上的建立图像I的高斯金字塔：{}0,...,LILLm=和图像J的高斯金字塔{}0,...,LJLLm=，其中mL是指金字塔的高度，通常取2,3,4。对于大小为640480×的图片，图片123,,III和4I的大小分别为320240,160120,8060,4030××××，在大多数情况下没有必要建立更高层数的金字塔。因为高斯金字塔图像中心的移动，使其可以处理光流大于综合窗口大小的特征点跟踪问题，由此可以根据估计的最大光流选择适当的图像金字塔的高度。在下一部分将详细的描述跟踪算法的细节以进一步的说明本文算法。注意：用公式（2）对图像进行下采样时中用到了低通滤波器[][]1/41/21/41/41/21/4T×防止图像失真。在实际工程中也可以使用平滑效果更强的低通滤波器[][]1/161/43/81/41/161/161/43/81/41/16T×建立图像金字塔，两者效果基本相同。 2.2 金字塔式的特征点跟踪 特征点跟踪的目的是给定图像I中的点u，计算图像J中对应的特征点[]Txyuu=+=vud，或则计算u点的光流d。对于0,...,mLL=，定义LLLxyuu⎡⎤=⎣⎦u是图像金字塔中图像LI的特征点u的坐标值。根据图像金字塔计算公式（2,3,4）可以计算Lu： 2LL=uu                                  （5） 公式（5）中的除法运算被应用在独立的两个坐标系中，同样相应的乘法运算将被应用在后续的公式中。特别的，0=uu。 金字塔跟踪算法总体流程如下：首先，计算金字塔最上层图像mLI上的光流。然后，根据最上层光流的计算结果估计次上层1mL−光流的初始值，再计算次上层图像上光流的精确值。最后，将次上层光流的计算结果估计下一层2mL−光流的初始值，计算其精确值后再带入下一层，直至最底层原始图像。 现用数学符号描述两个基本的图像层L和1L+间的递归操作。假设已经从第1L+层图像计算了第L层的光流初始估计值TLLLxygg⎡⎤=⎣⎦g。为了计算第L层的光流，有必要引入残差位移矢量TLLLxydd⎡⎤=⎣⎦d使得图像间的匹配误差函数Lε达到最小： ()()()2,(,)(,)LLxxxxLLxxxxuwuwLLLLLLLLLxyxxxyxuwyuwddIxyJxgdygdεε++=−=−==−++++∑∑d（6） 由公式（6）可以看出，每一层图像上综合窗口的大小都保持(21)(21)xyww+×+不变。注意应该在第二幅图像中以初始光流Lg值预先平移综合窗口，这样操作后计算的残差位移矢量TLLLxydd⎡⎤=⎣⎦d就非常小，则可通过标准的L‐K算法计算。 详细的计算残差位移矢量Ld的方法将在2.3中阐述。现在假设已经计算出残差位移矢量Ld，将其带入下一层1L−作为该层新的初始光流估计值1L−g。 12()LLL−=+ggd                          （7） 通过同样计算方法可以计算1L−层的残差位移矢量1L−d，此过程一直进行到图像最底层0L=（原始图像），此时图像和位移矢量都为原始的分辨率。假设算法中图像的最高层没有可靠的初始估计值，将其置为0矢量： []00mTL=g                           （8） 在最底层图像的位移矢量达到最佳分辨率时，可以计算公式（1）中光流的最终解d： 00=+dgd                              （9） 可以观测到光流矢量的解也可表示为下列展开式： 02mLLLL==∑dd                            （10） 由公式（10）可以看出，用金字塔实现L‐K算法的一个明显优点是：在计算一个很大的特征点平移矢量（光流矢量）d时，每一层的残差位移矢量Ld始终保持一个非常小的值。假设每一层元素的残差位移矢量都达到了maxd，那么金字塔算法可以处理的像素运动差和为1max_max(21)mLfinaldd+=−。比如，综合窗口大小为11×，高度为3mL=的金字塔，金字塔算法可以处理的最大像素差是15pixels。这意味着用相对小的综合窗口就能够处理大的像素平移。 2.3光流迭代计算法（迭代L­K） 现在描述如何计算光流。在图像金字塔的每一层L，目的是找到当层图像中的位移残差矢量Ld使得公式（6）中匹配误差函数()LLεd最小。因为每一层图像上的计算方式大致相同，故以第L层为例详细描述该层上位移残差矢量Ld的计算过程。定义如下两图像区域A和B： [](,)1,11,1xxxxyyyyxypwpwpwpw⎡⎤∀∈−−++×−−++⎣⎦, ()(),,LAxyIxy                          （11） [](,),,xxxxyyyyxypwpwpwpw⎡⎤∀∈−+×−+⎣⎦, ()(),,LLLxyBxyJxgyg++                    （12） 注意定义的(),Axy和(),Bxy区域有细微的差别。定义的(),Axy窗口大小是(23)(23)xyww+×+而不是(21)(21)xyww+×+。当使用中心差分原理计算(),Axy的空间微分（梯度）（公式（19,20））时，能够更加明显的体现两个区域范围不同及其不同的原因。令位移矢量TLxyvv⎡⎤==⎣⎦vd，特征点位置矢量TLxypp⎡⎤==⎣⎦pu，特征点跟踪的目的是计算位移矢量TLxyvv⎡⎤==⎣⎦vd使得匹配误差函数最小： ()()()2,(,)(,)yxxxxxyxpwpwxyxyxpwypwvvAxyBxvyvεε++=−=−==−++∑∑v         （13） 应用标准的迭代L‐K算法。当达到最优值时，ε关于v的一阶偏导数为0： ()[]00optε=∂=∂vvvv                       （14） 将倒数式展开得： ()()2(,)(,).yxxxxxyxpwpwxyxpwypwBBAxyBxvyvxyε++=−=−∂⎡⎤∂∂=−−++⎢⎥∂∂∂⎣⎦∑∑vv     （15） 现在将(,)xyBxvyv++用一阶泰勒展开式在点[]00T=v处展开（因为金字塔策略使所估计的平移矢量很小所以可以有效的近似）。 ()2(,)(,).yxxxxxyxpwpwxpwypwBBBBAxyBxyxyxyε++=−=−∂⎛⎞⎡⎤⎡⎤∂∂∂∂≈−−−⎜⎟⎢⎥⎢⎥∂∂∂∂∂⎣⎦⎣⎦⎝⎠∑∑vvv   （16） 注意到(,)(,)AxyBxy−的值可以被描述为在点[]Txy处图像关于时间的倒数： [](,),,xxxxyyyyxypwpwpwpw⎡⎤∀∈−+×−+⎣⎦ (),(,)(,)IxyAxyBxyδ−                    （17） 矩阵BBxy⎡⎤∂∂⎢⎥∂∂⎣⎦是图像的梯度矢量，重新定义表示符号： TxyIBBIIxy⎡⎤⎡⎤∂∂∇=⎢⎥⎢⎥∂∂⎣⎦⎣⎦                       （18） 注意图像梯度xI和yI可直接由图像(,)Axy在特征点p的(21)(21)xyww+×+邻域中计算，而与第二幅图像(,)Bxy无关。（在描述光流迭代算法时会显现这一处理的重要性。）用中心差分方法计算图像的梯度可表示为： [](,),,xxxxyyyyxypwpwpwpw⎡⎤∀∈−+×−+⎣⎦ (,)(1,)(1,)(,)2xAxyAxyAxyIxyx∂+−−==∂              （19） (,)(,1)(,1)(,)2yAxyAxyAxyIxyy∂+−−==∂              （20） 实际上也可用sharr方法计算图像的梯度（效果和中心差法接近）。公式（16）用新的符号表示为： ()()12yyxxxxyypwpwTTxpwypwIIIεδ++=−=−∂≈∇−∇∂∑∑vvv                  （21） ()22.1.2yyxxxxyypwpwxxxyxpwypwyxyyI