运动目标检测光流法

1、沈阳理工大学数字图像处理课程设计摘要运动目标检测方法是研究如何完成对视频图像序列中感兴趣的运动目标区 域的“准确定位”问题。光流场指图像灰度模式的表面运动，它可以反映视频相 邻帧之间的运动信息，因而可以用于运动目标的检测。MATLAB这种语言可移植性好、可扩展性强，再加上其中有丰富的图像处理函数，所以利用MATLAB软件来用光流法对运动目标的检测中具有很大的优势。本设计主要可以借助 matlab软件编写程序，运用 Horn-Schunek算法对图像前后两帧进行处理，画出 图像的光流场。而图像的光流场每个像素都有一个运动矢量，因此可以反映相邻 帧之间的运动，分析图像的光流场就可以得出图像中的运动

2、目标的运动情况。关键字：光流法；Horn-Schunek算法；matlab2目录1光流法的设计目的1.2光流法的原理1.2.1光流法的介绍1.2.1.1光流与光流场的概念 12.1光流法检测运动目标的原理 22.1.1光流场计算的基本原理 22.2.2基于梯度的光流场算法22.2.3Hor n-Schu nek 算法3.2.2.4光流法检测运动目标物体的基本原理概述 53光流法的程序具体实现 6.3.1源代码.6.3.1.1求解光流场函数.6.3.1.2求导函数9.3.1.3高斯滤波函数 .9.3.1.4平滑性约束条件函数 1.03.1.5画图函数104仿真图及分析1.2结论1.3参考文献14

3、沈阳理工大学数字图像处理课程设计1光流法的设计目的数字图像处理，就是用数字计算机及其他有关数字技术，对图像进行处理， 以达到预期的目的。随着计算机的发展，图像处理技术在许多领域得到了广泛应用，数字图像处理已成为电子信息、通信、计算机、自动化、信号处理 等专 业的重要课程。数字图像处理课程设计是在学习完数字图像处理的相关理论后，进行的综合性训练课程，其目的是：使学生进一步巩固数字图像处理的基本概念、理论、分 析方法和实现方法；增强学生应用Matlab编写数字图像处理的应用程序及分析、 解决实际问题的能力；尝试所学的内容解决实际工程问题， 培养学生的工程实践 能力。运动目标检测是数字图像处理技术的

4、一个主要部分，近些年来，随着多媒体技术的迅猛发展和计算机性能的不断提高， 动态图像处理技术日益受到人们的青 睞，并且取得了丰硕的成果，广泛应用于交通管理、军事目标跟踪、生物医学等 领域。因此，基于光流法，实现运动目标的检测是本文的研究对象。 结合图书馆书 籍、网上资料以及现有期刊杂志，初步建立起运动目标检测的整体思路和方法。2光流法的原理2.1光流法的介绍2.1.1光流与光流场的概念光流是指空间运动物体在观测成像面上的像素运动的瞬时速度， 它利用图像 序列像素强度数据的时域变化和相关性来确定各自像素位置的“运动” ，即反映 图像灰度在时间上的变化与景物中物体结构及其运动的关系。 将二维图像平面

5、特 定坐标点上的灰度瞬时变化率定义为光流矢量。视觉心理学认为人与被观察物体 发生相对运动时，被观察物体表面带光学特征的部位的移动给人们提供了运动和 结构的信息。当相机与场景目标间有相对运动时所观察到的亮度模式运动称之为 光流(optical flow)，或者说物体带光学特征部位的移动投影到视网膜平面(也即图 像平面)上就形成了光流。光流场是指图像灰度模式的表观运动， 它是一个二维矢量场，所包含的信息 就是各个像素点的瞬时运动速度矢量信息。光流场每个像素都有一个运动矢量，因此可以反映相邻帧之间的运动。2.1光流法检测运动目标的原理2.1.1光流场计算的基本原理一般情况下，光流由相机运动，场景中目

6、标运动或两者的共同运动产生。 光 流场的计算大致可分为三类：基于梯度的方法(Horn-Schunek)；基于匹配的方法； 基于能量的方法。基于梯度的方法利用图像灰度的梯度来计算光流，是研究最多的方法。基于梯度的方法根据运动前后图像灰度保持不变这个基本假设，导出光流约束方程。 由于光流约束方程并不能唯一的确定光流，因此需要导入其他的约束。根据引入的约束不同，基于梯度的方法又可以分为全局约束方法和局部约束方法。全局约束的方法假定光流在整个图像范围内满足一定的约束条件；而局部约束的方法假定在给定点周围的一个小区域内，光流满足一定的约束条件。基于匹配的方法，这类方法是将速度vm定义为视差d=(dx,d

7、y)t,使得两个时 刻的图像区域的匹配最佳。为了找到最佳匹配，我们可以对定义在d上的相似度 量，如规一化的互相关系数，进行最大化，也可以对某一距离度量，如光强度差 的平方和，进行最小化。2.2.2基于梯度的光流场算法梯度光流法又分为全局约束方法和局部约束方法。全局约束方法假定光流在 整个图像范围内满足一定的约束条件，而局部约束的方法假定在给定点周围的一 个小区域内，光流满足于一定的约束条件。下面先导出光流约束方程。然后给出 两种比较典型的基于梯度的方法。假定图像上点m =(x,y)在时刻t的灰度值为I (x, y,t),经过时间间隔dt 后，对应点的灰度为I (x dx, y dy, t dt

8、)3，当时，可以认为两点的灰度不变， 也就是：I ( x dx, y dy ,t dt)= I ( x, y,t)(2.1)如果图像灰度随x,y,t缓慢变化，可以将(1)式左边泰勒级数展开：I (x dx, y dy,t dt)I (x, y,t)亠亠：(2.2)其中&代表二阶无穷小项。由于dt-0,忽略&，可以得到：dx ldyVdt=0(2.3)x:x:t令u=, v=代表x,y方向上的光流，Ix=, Iy=, It=分别代表图像灰度相对 于x,y,t的偏导，式(2.3)可以写成：IxU lyV h =0(2.4)此式即光流场的基本方程。写成向量形式，即：Vm It =0(

9、2.5)其中=(Ix, Iy)是图像在点m处的梯度，Vm(U,V)是点m的光流。上式称为光 流约束方程，是所有基于梯度的光流计算方法的基础。2.2.3 Horn-Schu nek 算法Horn-Schunek算法提出了光流的平滑性约束。即：图像上任一点的光流并不 是独立的，光流在整个图像范围内平滑变化。因此 Horn-Schu nek算法是一种全 局约束的方法5。设平滑性约束项为极小化：Es = "(u：+v； +v：)dxdy(2.6)由基本等式，显然要求极小化：2E . (Ixu Iyy It) dxd y(2.7)于是，由(2.6)和(2.7)式可知，最后求得光流应满足(2.8

10、)式：152 2 2 2 2min 口丸(Ux +Uy +Vx 7)十(* +lyV+IJ dxdy(2.8)这里入的取值要考虑图中的噪声情况，如果噪声较强，说明图像数据本身的置信度较低，需要更多的依赖光流约束，所以可以取较大的值；反之，取较小的值。为了满足(2.8)，可将该式对u和v分别求导，并取导数为0。这样就得到5:l；u IxlyV2iu lxlt( 2.9)I；V IxlyU 二-円、v-lylt(2.10)以上两式也称为Euler方程。如果令U和V分别表示u邻域和v邻域中的均值(可用图像平滑算子求得)，并令' u=u-u和v=vv，则式(2.9)和(2.10 )改写成:h

11、= ( I i,j,k 1 ' l i 1, j,k 1' l i, j 1,k -1 ' l i 1, j 1 ,k 1)_( l i,j,k Ti 1,j,k ' l i,j 1,k ' l i 1, j 1,k)44(2.11)(I： +丸2) u + lxlyv = /u Ixlt(2.12)(I： +A2) v + lxlyU =fv_lylt从上式解得：(2.13)Ix(IxU lyv It)I, I,(2.14)ly(lxu lyvIt)2 l2 i：(2.15)式(2.14)和 (2.15)提供了用迭代法求解u和v的基础。实际中，常用松弛

12、迭方程进行求解:,(k 1)一(k)u-u- lx(k 1)一(k)v=7"I：(k)Ixu+一(k)lyvIt2I2 I 2xy(k)(k)Ixu+lyvIt2-| 2 | 21y(2.16)(2.17)(0) (0) 其中k是循环数，u和v是初始值，可以取为0。u和v是局部平均，' 为权重系数，根据导数求取的精确度确定。在实际求解过程中，需要估计亮度的时间和空间微分。 这可在图像点的一个 2X 2X 2立方邻域中估计，如果下标i，j，k分别对应x，y，t，那么3个一阶偏导 分别是：11lx ( 1 i 1,j,k 1 i 1,j 1,k 1 i 1,j,k 1 1 i 1

13、,j -1,k 1)- ( Ii,j,k 1 i,j 1,k 1 i,j,k d 1 i,j 1,k 1)44(2.18)111 y = ( 1 i,j 1,k ' 1 i 1,j 1,k' 1 i,j 1,k 1 ' 1 i 1,j -1,k 1)一 二(1 i,j,k 1 i -1,j,k ' 1 i,j,k ' 1 i 1,j,k 1)44(2.19)It 二(1 i,j ,k 11 i 1,j,k 1 ' 1 i,j 1,k 1 ' 1 i 1,j 1,k 1)- ( I i, j,k T i 1,j,k T i,j 1,k T

14、i 1,j 1,k)44(2.20)也就是用一阶差分来替代灰度I关于x,y,t轴的偏导。上述算法的实现相对简单，计算复杂性较低。但是这种技术存在着严重缺陷。 首先，图像灰度保持假设对于许多自然图像序列来讲都是不合适的，尤其是在图像的遮合边缘处和(或)当运动速度较高时，基于灰度保持假设的约束存在较大误 差。其次，在图像的遮合区域，速度场是突变的，而总体平滑约束则迫使所估计 的光流场平滑地穿过这一区域，此过程平滑掉了有关物体形状的非常重要的信 息。第二，微分技术的一个要求是l(x，y，t)必须是可微的，这暗示着需对图像 数据进行时空预平滑9，以避免混叠效应;而且数值微分的求取具有病态性，如果 处理

15、不当将对最终的速度估计产生显著影响。2.2.4光流法检测运动目标物体的基本原理概述给图像中的每一个像素点赋予一个速度矢量，这就形成了一个图像运动场， 在运动的一个特定时刻，图像上的点与三维物体上的点一一对应，这种对应关系 可由投影关系得到，根据各个像素点的速度矢量特征，可以对图像进行动态分析。 如果图像中没有运动物体，则光流矢量在整个图像区域是连续变化的。 当图像中 有运动物体时，目标与图像背景存在相对运动，运动物体所形成的速度矢量必然与邻域背景速度矢量不同，从而检测出运动物体及位置。光流法的优点：光流不仅携带了运动物体的运动信息， 而且携带了有关景物三维结构的丰富信息，它能够在不知道场景的任

16、何信息的情况下， 检测出运动对10象 。3光流法的程序具体实现下面，以光流场计算的典型方法Horn-Schunek算法位例，介绍光流法的Matlab程序实现对图像运动目标的检测。3.1源代码3.1.1求解光流场函数 fun cti on u, v = HS(im1, im2, alpha, ite, uln itial, vI nitial, displayFlow, displaylm)%输入：%im1-输入图像1%im2-输入图像2%alpha-反映HS光流算法的平滑性约束条件的参数%ita-( 3.15.9)式中的迭代次数%uInitial-光流横向分量初始值%vInitial-光流纵向

17、分量初始值%displayFlow-光流场显示参数，其值为1时显示，为0时不显示%displayImg-显示光流场的指定图像，如果为空矩阵，则无指定图像输出%输出：%u-横向光流矢量%v-纵向光流矢量%初始化参数g)%功能：求解光流场 if nargin<1 | nargin<2im1=imread('yos9.tif); im2=imread('yos10.tif);endif nargin<3 alpha=1;endif nargin<4ite=100;endif nargin<5 | nargin<6uln itial = zeros(

18、size(im1(:,:,1);vIn itial = zeros(size(im2(:,:,1);elseif size( uln itial,1) =0 | size( vln itial,1)=0uln itial = zeros(size(im1(:,:,1);vln itial = zeros(size(im2(:,:,1);endif nargin<7 displayFlow=1;endif nargin<8displaylmg=im1;end% 将RGB图像转化为灰度图像if size(size(im1),2)=3 im1=rgb2gray(im1);endif si

19、ze(size(im2),2)=3im2=rgb2gray(im2);endim仁double(iml);im2=double(im2);%调用平滑性约束函数对图像进行平滑im1=smoothlmg(im1,1);im2=smoothlmg(im2,1);tic;%为光流矢量设置初始值u = ul nitial;v = vl nitial;fx, fy, ft = computeDerivatives(im1, im2); %调用求导函数对时间分量和空间分量进行求导kernel_仁1/12 1/6 1/12;1/6 0 1/6;1/12 1/6 1/12; % 均值模板%迭代次数为100for

20、 i=1:ite %计算光流矢量的局部均值uAvg=c onv 2(u,kernel_1,'same');vAvg=conv2(v,kernel_1,'same'); % 用迭代法求解光流矢量u= uAvg - ( fx .* ( ( fx .* uAvg ) + ( fy .* vAvg ) + ft ) ) ./ ( alphaA2 + fx.A2 + fy.A2);v= vAvg - ( fy .* ( ( fx .* uAvg ) + ( fy .* vAvg ) + ft ) ) ./ ( alphaA2 + fx42 + fy.A2);endu(is

21、 nan(u)=0;v(is nan(v)=0;%画图if displayFlow=1plotFlow(u, v, displaylmg, 5, 5);% 调用画图函数end3.1.2求导函数 function fx, fy, ft = computeDerivatives(im1, im2)%功能：求输入图像参考像素点的像素值沿三轴方向的偏导数%输入：%im1-输入图像1%im2-输入图像2%输出：%fx-参考像素点的灰度值沿x方向的偏导数%fy-参考像素点的灰度值沿y方向的偏导数%fz-参考像素点的灰度值沿z方向的偏导数if size(im2,1)=0im2=zeros(size(im1)

22、;end%利用标准模板求得式(3.15.5)中的偏导数lx, ly, Itfx = conv2(im1,0.25* -1 1; -1 1,'same') + conv2(im2, 0.25*-1 1; -1 1,'same');fy = conv2(im1,0.25*-1 -1; 1 1, 'same') + conv2(im2, 0.25*-1 -1; 1 1, 'same');ft = con v2(im1,0.25*o nes(2),'same') + con v2(im2, -0.25*o nes(2),&

23、#39;same');3.1.3高斯滤波函数 fun cti on G=gaussFilter(segma,kSize)%功能：实现高斯滤波%输入：%sigma-高斯分布的概率密度函数的方差%kSize-高斯向量的模板尺寸大小%输出：%G-方差为segma,大小为kSize的一维高斯向量模板 if nargin<1segma=1;endif nargin<2kSize=2*(segma*3);endx=-(kSize/2):(1+1/kSize):(kSize/2);%利用均值为0,方差为segma高斯分布概率密度函数求解一维高斯向量模板G=(1/(sqrt(2*pi)*s

24、egma) * exp (-(x.A2)/(2*segmaA2);3.1.4平滑性约束条件函数fun cti on smoothedlmg=smoothlmg(img,segma)%功能：实现平滑性约束条件%输入：%img-数字图像%sigma-高斯分布的方差%输出：%smoothedlmg经高斯滤波的图像矩阵if nargin<2segma=1;end%调用高斯滤波函数G=gaussFilter(segma);%对图像进行平滑约束smoothedlmg=c onv 2(img,G'same');%二次平滑smoothedlmg=c onv 2(smoothedlmg,G

25、','same');3.1.5画图函数fun cti on plotFlow(u, v, imgOrigi nal, rSize, scale)%功能：绘制光流场图%输入：u-横向光流矢量v-纵向光流矢量imgOriginal-光流场显示的图像%rSize-可见光流矢量区域的尺寸scale-光流场规模figure();if nargin>2if sum(sum(imgOrigi nal)=0imshow(imgOriginal,0 255);hold on;endendif nargin<4rSize=5;endif nargin<5scale=3;endfor i=1:size(u,1)for j=1:size(u,2)if floor(i/rSize)=i/rSize | floor(j/rSize)=j/rSizeu(i,j)=0;v(i,j)=0;endendendquiver(u, v, scale, 'color', 'g', 'linewidth', 2);set(gca,'YD,'reverse');4仿真图及分析将上述程序保存为matlab.