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;

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----

?????文件???????2490??2009-05-08?16:46??optical_flow.m

-----------?---------??----------?-----??----

?????????????????2490????????????????????1