optiflow 光流法程序

光流法程序，利用传统光流法思想，亲测有效，效果不错

function [usvs] = HSoptflow（Xrgbn）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author: Gregory Power gregory.power@wpafb.af.mil
% This MATLAB code shows a Motion Estimation map created by
% using a Horn and Schunck motion estimation technique on two
% consecutive frames.  Input requires.
%     Xrgb（hwdN） where X is a frame sequence of a certain
%                height（h） width （w） depth （d=3 for color） 
%                and number of frames （N）. 
%     n= is the starting frame number which is less than N 
%     V= the output variable which is a 2D velocity array
%
% Sample Call: V=HSoptflow（X3）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[hwdN]=size（Xrgb）;
if n>N-1
   error（1‘requested file greater than frame number required‘）;
end
%get two image frames
if d==1
    Xn=double（Xrgb（::1n））;
    Xnp1=double（Xrgb（::1n+1））;
elseif d==3
    Xn=double（Xrgb（::1n）*0.299+Xrgb（::2n）*0.587+Xrgb（::3n）*0.114）;
    Xnp1=double（Xrgb（::1n+1）*0.299+Xrgb（::2n+1）*0.587+Xrgb（::3n+1）*0.114）;
else
    error（2‘not an RGB or Monochrome image file‘）;
end

%get image size and adjust for border
% size（Xn）
hm5=h-5; wm5=w-5;
z=zeros（hw）; v1=z; v2=z;

%initialize
dsx2=v1; dsx1=v1; dst=v1;
alpha2=625;%%%%平衡两个误差的常数
imax=20;%%%%%迭代最大次数

%Calculate gradients
dst（5:hm55:wm5） = （ Xnp1（6:hm5+16:wm5+1）-Xn（6:hm5+16:wm5+1） + Xnp1（6:hm5+15:wm5）-Xn（6:hm5+15:wm5）+ Xnp1（5:hm56:wm5+1）-Xn（5:hm56:wm5+1） +Xnp1（5:hm55:wm5）-Xn（5:hm55:wm5））/4;%%%对应点梯度，即随时间变化的梯度
dsx2（5:hm55:wm5） = （ Xnp1（6:hm5+16:wm5+1）-Xnp1（5:hm56:wm5+1） + Xnp1（6:hm5+15:wm5）-Xnp1（5:hm55:wm5）+ Xn（6:hm5+16:wm5+1）-Xn（5:hm56:wm5+1） +Xn（6:hm5+15:wm5）-Xn（5:hm55:wm5））/4;%%%纵向梯度
dsx1（5:hm55:wm5） = （ Xnp1（6:hm5+16:wm5+1）-Xnp1（6:hm5+15:wm5） + Xnp1（5:hm56:wm5+1）-Xnp1（5:hm55:wm5）+ Xn（6:hm5+16:wm5+1）-Xn（6:hm5+15:wm5） +Xn（5:hm56:wm5+1）-Xn（5:hm55:wm5））/4;%%%横向梯度


for i=1:imax
   delta=（dsx1.*v1+dsx2.*v2+dst）./（alpha2+dsx1.^2+dsx2.^2）;
   v1=v1-dsx1.*delta;
   v2=v2-dsx2.*delta;
end
u=z; u（5:hm55:wm5）=v1（5:hm55:wm5）;
v=z; v（5:hm55:wm5）=v2（5:hm55:wm5）;

xskip=round（h/32）;
[hsws]=size（u（1:xskip:h1:xskip:w））;
us=zeros（hsws）; vs=us;

N=xskip^2;

for i=1:hs-1
  for j=1:ws-1
     hk=i*xskip-xskip+1;
     hl=i*xskip;
     wk=j*xskip-xskip+1;
     wl=j*xskip;
     us（ij）=sum（sum（u（hk:hlwk:wl）））/N;
     vs（ij）=sum（sum（v（hk:hlwk:wl）））/N;
  end
end

figure（1）;
quiver（usvs）;
colormap（‘default‘）;
axis ij;
axis tight;
axis equal;

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件      51271  2011-07-21 14:17  optiflow\1.jpg

     文件      51379  2011-07-21 14:17  optiflow\10.jpg

     文件      51210  2011-07-21 14:17  optiflow\11.jpg

     文件      51337  2011-07-21 14:17  optiflow\12.jpg

     文件      51169  2011-07-21 14:17  optiflow\13.jpg

     文件      51281  2011-07-21 14:17  optiflow\14.jpg

     文件      51367  2011-07-21 14:17  optiflow\15.jpg

     文件      51297  2011-07-21 14:17  optiflow\16.jpg

     文件      51305  2011-07-21 14:17  optiflow\17.jpg

     文件      51303  2011-07-21 14:17  optiflow\2.jpg

     文件      51393  2011-07-21 14:17  optiflow\3.jpg

     文件      51406  2011-07-21 14:17  optiflow\4.jpg

     文件      51413  2011-07-21 14:17  optiflow\5.jpg

     文件      51313  2011-07-21 14:17  optiflow\6.jpg

     文件      51406  2011-07-21 14:17  optiflow\7.jpg

     文件      51294  2011-07-21 14:17  optiflow\8.jpg

     文件      51331  2011-07-21 14:17  optiflow\9.jpg

     文件       2589  2012-07-31 12:30  optiflow\HSoptflow.m

     文件        508  2012-07-18 13:42  optiflow\main.asv

     文件        612  2012-07-18 15:57  optiflow\main.m

    ..A.SH.     67584  2012-07-18 14:06  optiflow\Thumbs.db

     目录          0  2012-07-17 23:15  optiflow

----------- ---------  ---------- -----  ----

               943768                    22