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