[发明专利]噪声环境下光流场快速稳健估计方法

摘要

本发明涉及一种噪声环境下光流场快速稳健估计方法。所提算法基于噪声环境下光流场估计方法（ML法），引入惩罚因子以增强光流计算稳健性，并在光流计算迭代公式中加入动量因子缩短光流计算收敛时间以加快光流场计算，而后基于变分方法极小化光流能量函数求解欧拉‑拉格朗日方程，最后通过迭代方法求得速度场。仿真结果表明，对视频中连续两帧图片加入不同高斯噪声后，与M算法及ML方法相比，所提算法可显著增强光流场计算稳健性，缩短光流计算收敛时间，加快光流场计算。

主权项

噪声环境下光流场快速稳健估计方法，其特征在于，该方法包括以下步骤：步骤1：光流约束方程假设图像上一点(x,y)在t时刻的亮度为I(x,y,t)，在Δt时间后该像素点亮度变为I(x+Δx,y+Δy,t+Δt)，当Δt趋于无穷小时认为该点亮度不变，得到等式：∂I∂xΔxΔt+∂I∂yΔyΔt+∂I∂x=0---(1)]]>将改写为Ix,Iy,It，并令则重写上式，得到光流计算的基本等式：Ixu+Iyv+It＝0     (2)u,v表示速度场矢量的两个分量；步骤2：构建能量函数方程①光流约束方程修正设物体面元由位置P运动到位置P'，对应成像投影在像平面上由p移动到p'，假设物体表面为朗伯面，则p和p'的灰度分别为：Ip=ρdIq(N·L)=ρdIqcosθIp′=ρdIq(N′·L)=ρdIqcosθ′---(3)]]>其中，ρd为物体表面漫反射系数，Iq为入射光强度，设为常数；N和N'分别为物体面元在位置P和P'的单位法向量；L为入射光的单位方向向量；θ为N和L的夹角；θ'为N'和L的夹角；光流基本方程的修正项为：q=I(x+dx,y+dy,t+dt)-I(x,y,t)=ρdIq[(N′-N)·L]=ρdIq(ΔN·L)=cI(x,y,t)---(4)]]>其中：ΔN表示曲面法线运动变化；从而将基本方程修正为：cI+Ixu+Iyv+It＝0      (5)定义光流约束因子为：e≈(cI+Ixu+Iyv+It)2     (6)其中，c为参量，连同u,v采用全部平滑约束；全部平滑约束用梯度模的平方和来测度；定义光流场的全局平滑约束因子为：e1＝||▽u||2+||▽v||2＝(ux2+vx2+uy2+vy2)     (7)定义参量c的全局平滑约束因子为:e2＝||▽c||2＝(cx2+cy2)       (8)设I'(x,y,z)为I(x,y,z)受加性噪声n(x,y,z)污染后图像点(x,y)的实际灰度值，则：I(x,y,z)＝I(x,y,z)‑n(x,y,z)     (9)假设式(9)中的噪声n(x,y,z)是均值为0，方差为δ2的高斯白噪声，则：∫∫ΩI(x,y,t)dxdy=∫∫Ω[I′(x,y,t)-n(x,y,t)]dxdy=∫∫ΩI′(x,y,t)dxdy-∫∫Ωn(x,y,t)dxdy---(10)]]>∫∫Ω[I(x,y,t)-I′(x,y,t)]2dxdy=∫∫Ωn(x,y,t)dxdy=δ2---(11)]]>I(x,y,z)((x,y)∈Ω)未受噪声污染，故应满足全局平滑约束，定义灰度全局平滑约束因子为：e3＝||▽I||2＝(Ix2+Iy2)     (12)②引入惩罚项引入二次项β2(u2+v2+c2)作为惩罚项：e4＝β2(u2+v2+c2)      (13)噪声环境下光流场估计问题转化为式(11)约束下极小化：E=∫∫Ω(e+α1e1+α2e2+α3e3+e4)dxdy---(14)]]>其中，α1,α2,α3为调和4类约束因子的权参数；步骤3：求解光流场①光流求解迭代公式根据拉格朗日乘数法，式(11)约束下极小化式(14)即是极小化：∫∫Ω([e+α1e1+α2e2+α3e3+e4+λ(I-I′)-λδ2])dxdy---(15)]]>令：F(x,y,I,u,v,c,Ix,Iy,ux,uy,vx,vy,cx,cy)＝e+α1e1+α2e2+α3e3+e4+λ(I‑I')‑λδ2    (16)由变分法，式(15)取极小值的必要条件为：∂F∂I-∂∂x(∂F∂Ix)-∂∂y(∂F∂Iy)=0∂F∂u-∂∂x(∂F∂ux)-∂∂y(∂F∂uy)=0∂F∂v-∂∂x(∂F∂vx)-∂∂y(∂F∂vy)=0∂F∂c-∂∂x(∂F∂cx)-∂∂y(∂F∂cy)=0---(17)]]>将式(16)带入式(17)，化简后得：c(Ixu+Iyv+It+cI)=α3(Ixx+Iyy)+∂∂x[u(Ixu+Iyv+It+cI)]+∂∂y[v(Ixu+Iyv+It+cI)]-λ(I-I′)................(a)Ix(Ixu+Iyv+It+cI)=α1(uxx+uyy)+β2u...........(b)Iy(Ixu+Iyv+It+cI)=α1(vxx+vyy)+β2v............(c)I(Ixu+Iyv+It+cI)=α2(cxx+cyy)+β2c............(d)---(18)]]>其中，sxx+syy＝▽2s(s＝I,u,v,c)为s分量的拉普拉斯算子，采用下式估算：▿2s≈K(s‾i,j,k-si,j,k)---(19)]]>和si,j,k分别为tk时刻(xi,yi)处s分量的取值及其邻域平均值，K＝1；定义为：s‾i,j,k=16(si-1,j,k+si,j-1,k+si,j+1,k+si+1,j,k)+112(si-1,j-1,k+si-1,j+1,k+si+1,j-1,k+si+1,j+1,k)---(20)]]>由式(18)式中(d)可得：(Ixu+Iyv+It+cI)＝α2(cxx+cyy)/I     (21)化简得：(B2+B5+B6+B7)(I-I‾)=-[c(A1I‾+A2+A6)+A4I‾+A5I‾]-D(B1+B2+B3+B4)(u-u‾)=-[Ix(A1I‾+A2+A3)+A4I‾+A5I‾](B1+B2+B3+B4)(v-v‾)=-[Iy(A1I‾+A2+A3)+A4I‾+A5I‾](B1+B2+B3+B4)(u-u‾)=-[I(A1I‾+A2+A6)+A4I‾+A5I‾]---(22)]]>其中：D=λ(I-I′)+α2(Ix2+Iy2+α1)[(I2-II‾)(ux+uy)+(I‾Ix-II‾x)u+(I‾Iy-II‾y)v]/I2---(23)]]>式(23)两端乘以(I‑I')并在图像平面Ω上积分得：∫∫ΩD(I-I′)2dxdy=λ∫∫Ω(I-I′)2dxdy+∫∫Ωα2(Ix2+Iy2+α1)(I-I′)[(I2-II‾)(ux+uy)+(I‾Ix-II‾x)u+(I‾Iy-II‾y)v]/I2dxdy---(24)]]>将式(11)带入式(24)可得:λ=1δ2∫∫Ω[D-α2(Ix2+Iy2+α1)[(I2-II‾)(ux+uy)+(I‾Ix-II‾x)u+(I‾Iy-II‾y)v]/I2](I-I′)dxdy---(25)]]>由式(22)和(25)进一步可得：λn+1=1δ2∫∫Ω{Dn-α2((Ixn)2+(Iyn)2+α1)[((In)2-InI‾n)(ux+uy)+(I‾nIxn-InI‾xn)u+(I‾nIyn-InI‾yn)v]/(In)2}(In-I′)dxdy---(26)]]>In+1=I‾n-cn(A1I‾n+A2+A6)+A4I‾n+A5I‾n+DnB2+B5+B6+B7un+1=u‾n-Ixn(A1u‾n+A2+A3)+A4u‾n+A5u‾nB1+B2+B3+B4vn+1=v‾n-Iyn(A1v‾n+A2+A3)+A4v‾n+A5v‾nB1+B2+B3+B4cn+1=c‾n-In(A1c‾n+A2+A6)+A4c‾n+A5c‾nB1+B2+B3+B4---(27)]]>②添加动量因子后的光流求解迭代公式在光流迭代公式中添加动量因子以加快光流场计算，即：In+1=I‾n-cn(A1I‾n+A2+A6)+A4I‾n+A5I‾n+DnB2+B5+B6+B7un+1=u‾n-Ixn(A1u‾n+A2+A3)+A4u‾n+A5u‾nB1+B2+B3+B4+μn(uin-uin-1)vn+1=v‾n-Iyn(A1v‾n+A2+A3)+A4v‾n+A5v‾nB1+B2+B3+B4+μn(vin-vin-1)cn+1=c‾n-In(A1c‾n+A2+A6)+A4c‾n+A5c‾nB1+B2+B3+B4---(28)]]>其中：A1=Ixn(β2α2+α1-β4+α1α2)A2=Itn(β4+α1α2-α1β2-α2β2)A3=c‾nα2In(α1-β2)+α2v‾nIyn(α2-β2)A4=β2[α2(Iyn)2-β2(Iyn)2-α1+α1α2]A5=β4(1-α1-α2+β2)A6=c‾nα1In(α1-β2)+α1v‾nIyn(α2-β2)---(29)]]>B1=-α1α2[(Ixn)2+(Iyn)2+(In)2]B2=α1β2[(Ixn)2+(Iyn)2+(In)2]B3=α12(1-α2+β2)B4=α1β2(α2-β2-1)B5=-α1α3[(Ixn)2+(Iyn)2+(In)2]B6=α12(1-α3+β2)B7=α1β2(α3-β2-1)---(30)]]>其中，在对图像求x方向、y方向和t方向的导数时，采用的公式如下：sx=14{[s(xi,yj+1,tk-1)-s(xi,yj,tk-1)]+[s(xi+1,yj+1,tk-1)-s(xi+1,yj,tk-1)]+[s(xi,yj+1,tk)-s(xi,yj,tk)]+[s(xi+1,yj+1,tk)-s(xi+1,yj,tk)]}(s=I,u,v)---(31)]]>sy=14{[s(xi+1,yj,tk-1)-s(xi,yj,tk-1)]+[s(xi+1,yj+1,tk-1)-s(xi,yj+1,tk-1)]+[s(xi+1,yj,tk)-s(xi,yj,tk)]+[s(xi+1,yj+1,tk)-s(xi,yj+1,tk)]}(s=I,u,v)---(32)]]>st=14{[s(xi,yj,tk)-s(xi,yj,tk-1)]+[s(xi+1,yj,tk)-s(xi+1,yj,tk-1)]+[s(xi,yj+1,tk)-s(xi,yj+1,tk-1)]+[s(xi+1,yj+1,tk)-s(xi+1,yj+1,tk-1)]}(s=I)---(33)]]>步骤4：迭代求解根据式(28)，用Gauss‑Seidel方法进行迭代，并计算两次迭代光流值的误差，如果小于0.005，或者迭代次数超过200次，则此时所得即为所求光流值。