function [us,vs,hs,ws] = HSoptflow(shot,n)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author: Arron_Tu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
h=shot.height;
w=shot.width;
d=3;
N=size(shot.frames,2);
if n>N-1
   error(1,'requested file greater than frame number required');
end; 
%get two image frames
shot.frames(1,n).cdata=double(shot.frames(1,n).cdata);
shot.frames(1,n+1).cdata=double(shot.frames(1,n+1).cdata);
if d==1
    Xn=shot.frames(1,n).cdata(:,:,1);
    Xnp1=shot.frames(1,n+1).cdata(:,:,1);
elseif d==3
    Xn=shot.frames(1,n).cdata(:,:,1)*0.299+shot.frames(1,n).cdata(:,:,2)*0.587+shot.frames(1,n).cdata(:,:,3)*0.114;
    Xnp1=shot.frames(1,n+1).cdata(:,:,1)*0.299+shot.frames(1,n+1).cdata(:,:,2)*0.587+shot.frames(1,n+1).cdata(:,:,3)*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(h,w); v1=z; v2=z;

%initialize
dsx2=v1; dsx1=v1; dst=v1;
alpha2=625;
imax=20;

%Calculate gradients
dst(5:hm5,5:wm5) = ( Xnp1(6:hm5+1,6:wm5+1)-Xn(6:hm5+1,6:wm5+1) + Xnp1(6:hm5+1,5:wm5)-Xn(6:hm5+1,5:wm5)+... 
                     Xnp1(5:hm5,6:wm5+1)-Xn(5:hm5,6:wm5+1) +Xnp1(5:hm5,5:wm5)-Xn(5:hm5,5:wm5))/4;
dsx2(5:hm5,5:wm5) = ( Xnp1(6:hm5+1,6:wm5+1)-Xnp1(5:hm5,6:wm5+1) + Xnp1(6:hm5+1,5:wm5)-Xnp1(5:hm5,5:wm5)+...
                     Xn(6:hm5+1,6:wm5+1)-Xn(5:hm5,6:wm5+1) +Xn(6:hm5+1,5:wm5)-Xn(5:hm5,5:wm5))/4;
dsx1(5:hm5,5:wm5) = ( Xnp1(6:hm5+1,6:wm5+1)-Xnp1(6:hm5+1,5:wm5) + Xnp1(5:hm5,6:wm5+1)-Xnp1(5:hm5,5:wm5)+...
                     Xn(6:hm5+1,6:wm5+1)-Xn(6:hm5+1,5:wm5) +Xn(5:hm5,6:wm5+1)-Xn(5:hm5,5: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:hm5,5:wm5)=v1(5:hm5,5:wm5);
v=z; v(5:hm5,5:wm5)=v2(5:hm5,5:wm5);

xskip=round(h/32);
[hs,ws]=size(u(1:xskip:h,1:xskip:w));
us=zeros(hs,ws); 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(i,j)=sum(sum(u(hk:hl,wk:wl)))/N;
     vs(i,j)=sum(sum(v(hk:hl,wk:wl)))/N;
   end;
end;

figure(1);
quiver(us,vs);
colormap('default');
axis ij;
axis tight;
axis equal;