运动目标检测光流法

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

2、情况。关键字：光流法；Horn-Schunck算法；matlab目 录1光流法的设计目的12光流法的原理12.1光流法的介绍12.1.1光流与光流场的概念12.1光流法检测运动目标的原理22.1.1光流场计算的基本原理22.2.2基于梯度的光流场算法22.2.3Horn-Schunck算法32.2.4光流法检测运动目标物体的基本原理概述53光流法的程序具体实现63.1源代码63.1.1求解光流场函数63.1.2求导函数83.1.3高斯滤波函数93.1.4平滑性约束条件函数103.1.5画图函数104仿真图及分析12结论13参考文献14I1 光流法的设计目的数字图像处理，就是用数字计算机及其他有

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

4、们的青睞，并且取得了丰硕的成果，广泛应用于交通管理、军事目标跟踪、生物医学等领域。因此，基于光流法，实现运动目标的检测是本文的研究对象。结合图书馆书籍、网上资料以及现有期刊杂志，初步建立起运动目标检测的整体思路和方法。2 光流法的原理2.1 光流法的介绍2.1.1 光流与光流场的概念光流是指空间运动物体在观测成像面上的像素运动的瞬时速度，它利用图像序列像素强度数据的时域变化和相关性来确定各自像素位置的“运动”，即反映图像灰度在时间上的变化与景物中物体结构及其运动的关系。将二维图像平面特定坐标点上的灰度瞬时变化率定义为光流矢量。视觉心理学认为人与被观察物体发生相对运动时，被观察物体表面带光学特征

5、的部位的移动给人们提供了运动和结构的信息。当相机与场景目标间有相对运动时所观察到的亮度模式运动称之为光流(optical flow)，或者说物体带光学特征部位的移动投影到视网膜平面(也即图像平面)上就形成了光流。光流场是指图像灰度模式的表观运动，它是一个二维矢量场，所包含的信息就是各个像素点的瞬时运动速度矢量信息。光流场每个像素都有一个运动矢量，因此可以反映相邻帧之间的运动1。2.1 光流法检测运动目标的原理2.1.1 光流场计算的基本原理一般情况下，光流由相机运动，场景中目标运动或两者的共同运动产生。光流场的计算大致可分为三类：基于梯度的方法(Horn-Schunck)；基于匹配的方法；基于

6、能量的方法。基于梯度的方法利用图像灰度的梯度来计算光流，是研究最多的方法。基于梯度的方法根据运动前后图像灰度保持不变这个基本假设，导出光流约束方程。由于光流约束方程并不能唯一的确定光流，因此需要导入其他的约束。根据引入的约束不同，基于梯度的方法又可以分为全局约束方法和局部约束方法。全局约束的方法假定光流在整个图像范围内满足一定的约束条件；而局部约束的方法假定在给定点周围的一个小区域内，光流满足一定的约束条件。基于匹配的方法，这类方法是将速度vm定义为视差d=(dx,dy)t，使得两个时刻的图像区域的匹配最佳。为了找到最佳匹配，我们可以对定义在d上的相似度量，如规一化的互相关系数，进行最大化，也

7、可以对某一距离度量，如光强度差的平方和，进行最小化2。2.2.2 基于梯度的光流场算法梯度光流法又分为全局约束方法和局部约束方法。全局约束方法假定光流在整个图像范围内满足一定的约束条件，而局部约束的方法假定在给定点周围的一个小区域内，光流满足于一定的约束条件。下面先导出光流约束方程。然后给出两种比较典型的基于梯度的方法。假定图像上点在时刻t的灰度值为，经过时间间隔dt后，对应点的灰度为3，当时，可以认为两点的灰度不变，也就是： （2.1）如果图像灰度随x,y,t缓慢变化，可以将(1)式左边泰勒级数展开： （2.2）其中代表二阶无穷小项。由于dt0，忽略，可以得到： （2.3） 令u=，v=代表

8、x,y方向上的光流，IX=，IY=，IT=分别代表图像灰度相对于x,y,t的偏导4，式(2.3)可以写成： （2.4）此式即光流场的基本方程。写成向量形式，即： （2.5）其中=(IX, Iy)是图像在点m处的梯度，Vm(U,V)是点m的光流。上式称为光流约束方程，是所有基于梯度的光流计算方法的基础。 2.2.3 Horn-Schunck算法Horn-Schunck算法提出了光流的平滑性约束。即：图像上任一点的光流并不是独立的，光流在整个图像范围内平滑变化。因此Horn-Schunck算法是一种全局约束的方法5。设平滑性约束项为极小化： （2.6）由基本等式，显然要求极小化： （2.7）于是，

9、由(2.6)和(2.7)式可知，最后求得光流应满足（2.8）式： （2.8）这里的取值要考虑图中的噪声情况，如果噪声较强，说明图像数据本身的置信度较低，需要更多的依赖光流约束，所以可以取较大的值；反之，取较小的值。为了满足(2.8)，可将该式对u和v分别求导，并取导数为0。这样就得到5： （2.9） （2.10） 以上两式也称为Euler方程。如果令和分别表示u邻域和v邻域中的均值(可用图像平滑算子求得)6，并令和，则式(2.9)和(2.10)改写成： （2.11） （2.12） （2.13） 从上式解得： （2.14） （2.15） 式(2.14)和(2.15)提供了用迭代法求解u和v的基础

10、7。实际中，常用松弛迭方程进行求解： （2.16） （2.17） 其中k是循环数，和是初始值，可以取为0。和是局部平均，为权重系数，根据导数求取的精确度确定。 在实际求解过程中，需要估计亮度的时间和空间微分。这可在图像点的一个222立方邻域中估计8，如果下标分别对应，那么3个一阶偏导分别是： （2.18） （2.19） （2.20）也就是用一阶差分来替代灰度I关于x,y,t轴的偏导。上述算法的实现相对简单，计算复杂性较低。但是这种技术存在着严重缺陷。首先，图像灰度保持假设对于许多自然图像序列来讲都是不合适的，尤其是在图像的遮合边缘处和(或)当运动速度较高时，基于灰度保持假设的约束存在较大误差。

11、其次，在图像的遮合区域，速度场是突变的，而总体平滑约束则迫使所估计的光流场平滑地穿过这一区域，此过程平滑掉了有关物体形状的非常重要的信息。第二，微分技术的一个要求是I(x，y，t)必须是可微的，这暗示着需对图像数据进行时空预平滑9，以避免混叠效应;而且数值微分的求取具有病态性，如果处理不当将对最终的速度估计产生显著影响。2.2.4 光流法检测运动目标物体的基本原理概述给图像中的每一个像素点赋予一个速度矢量，这就形成了一个图像运动场，在运动的一个特定时刻，图像上的点与三维物体上的点一一对应，这种对应关系可由投影关系得到，根据各个像素点的速度矢量特征，可以对图像进行动态分析。如果图像中没有运动物体

12、，则光流矢量在整个图像区域是连续变化的。当图像中有运动物体时，目标与图像背景存在相对运动，运动物体所形成的速度矢量必然与邻域背景速度矢量不同，从而检测出运动物体及位置。光流法的优点：光流不仅携带了运动物体的运动信息，而且携带了有关景物三维结构的丰富信息，它能够在不知道场景的任何信息的情况下，检测出运动对象10。3 光流法的程序具体实现下面，以光流场计算的典型方法Horn-Schunck算法位例，介绍光流法的Matlab程序实现对图像运动目标的检测。3.1 源代码3.1.1 求解光流场函数function u, v = HS(im1, im2, alpha, ite, uInitial, vIn

13、itial, displayFlow, displayIm) % 输入：% im1-输入图像1% im2-输入图像2% alpha-反映HS光流算法的平滑性约束条件的参数% ita-（3.15.9）式中的迭代次数% uInitial-光流横向分量初始值% vInitial-光流纵向分量初始值% displayFlow-光流场显示参数，其值为1时显示，为0时不显示% displayImg-显示光流场的指定图像，如果为空矩阵，则无指定图像输出% 输出：% u-横向光流矢量% v-纵向光流矢量 % 初始化参数g)% 功能：求解光流场 if nargin1 | nargin2 im1=imread(y

14、os9.tif); im2=imread(yos10.tif);endif nargin3 alpha=1;endif nargin4 ite=100;endif nargin5 | nargin6 uInitial = zeros(size(im1(:,:,1); vInitial = zeros(size(im2(:,:,1);elseif size(uInitial,1) =0 | size(vInitial,1)=0 uInitial = zeros(size(im1(:,:,1); vInitial = zeros(size(im2(:,:,1);endif nargin7 disp

15、layFlow=1;endif nargin8 displayImg=im1;end% 将RGB图像转化为灰度图像if size(size(im1),2)=3 im1=rgb2gray(im1);endif size(size(im2),2)=3 im2=rgb2gray(im2);endim1=double(im1);im2=double(im2);% 调用平滑性约束函数对图像进行平滑im1=smoothImg(im1,1);im2=smoothImg(im2,1);tic;% 为光流矢量设置初始值u = uInitial;v = vInitial;fx, fy, ft = computeD

16、erivatives(im1, im2); % 调用求导函数对时间分量和空间分量进行求导kernel_1=1/12 1/6 1/12;1/6 0 1/6;1/12 1/6 1/12; % 均值模板% 迭代次数为100for i=1:ite % 计算光流矢量的局部均值 uAvg=conv2(u,kernel_1,same); vAvg=conv2(v,kernel_1,same); % 用迭代法求解光流矢量 u= uAvg - ( fx .* ( ( fx .* uAvg ) + ( fy .* vAvg ) + ft ) ) ./ ( alpha2 + fx.2 + fy.2); v= vAv

17、g - ( fy .* ( ( fx .* uAvg ) + ( fy .* vAvg ) + ft ) ) ./ ( alpha2 + fx.2 + fy.2);endu(isnan(u)=0;v(isnan(v)=0;% 画图if displayFlow=1 plotFlow(u, v, displayImg, 5, 5); % 调用画图函数end3.1.2 求导函数function fx, fy, ft = computeDerivatives(im1, im2)% 功能：求输入图像参考像素点的像素值沿三轴方向的偏导数% 输入：% im1-输入图像1% im2-输入图像2% 输出：% f

18、x-参考像素点的灰度值沿x方向的偏导数% fy-参考像素点的灰度值沿y方向的偏导数% fz-参考像素点的灰度值沿z方向的偏导数 if size(im2,1)=0 im2=zeros(size(im1);end % 利用标准模板求得式（3.15.5）中的偏导数Ix, Iy, 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 = conv2(im

19、1, 0.25*ones(2),same) + conv2(im2, -0.25*ones(2),same);3.1.3 高斯滤波函数function G=gaussFilter(segma,kSize)% 功能：实现高斯滤波% 输入：% sigma-高斯分布的概率密度函数的方差% kSize-高斯向量的模板尺寸大小% 输出： % G-方差为segma，大小为kSize的一维高斯向量模板 if nargin1 segma=1;endif nargin2 kSize=2*(segma*3);end x=-(kSize/2):(1+1/kSize):(kSize/2);% 利用均值为0，方差为se

20、gma高斯分布概率密度函数求解一维高斯向量模板G=(1/(sqrt(2*pi)*segma) * exp (-(x.2)/(2*segma2);3.1.4 平滑性约束条件函数function smoothedImg=smoothImg(img,segma)% 功能：实现平滑性约束条件% 输入：% img-数字图像% sigma-高斯分布的方差% 输出：% smoothedImg-经高斯滤波的图像矩阵if nargin2 if sum(sum(imgOriginal)=0 imshow(imgOriginal,0 255); hold on; endendif nargin4 rSize=5;e

21、ndif nargin5 scale=3;endfor i=1:size(u,1) for j=1:size(u,2) if floor(i/rSize)=i/rSize | floor(j/rSize)=j/rSize u(i,j)=0; v(i,j)=0; end endendquiver(u, v, scale, color, g, linewidth, 2);set(gca,YDir,reverse);4 仿真图及分析将上述程序保存为matlab.m格式，并输入两帧连续图像，如图4.1所示： 图4.1输入的两帧测试图像由图4.1可知：输入的两帧连续图像为灰度图像，并且有运动的目标即河中水流和空气。 图4.2两帧测试图像的光流场及放大图像由图4.2可知：输入的两帧连续图像，在经过求导函数、高斯滤波函数、平滑条件函数处理之后，用光流法将图像中运动物体的速度矢