1 function [us,vs] = HSoptflow(Xrgb,n)
 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 3 % Author: Gregory Power [email protected]
 4 % This MATLAB code shows a Motion Estimation map created by
 5 % using a Horn and Schunck motion estimation technique on two
 6 % consecutive frames.  Input requires.
 7 %     Xrgb(h,w,d,N) where X is a frame sequence of a certain
 8 %                height(h), width (w), depth (d=3 for color), 
 9 %                and number of frames (N). 
10 %     n= is the starting frame number which is less than N 
11 %     V= the output variable which is a 2D velocity array
12 %
13 % Sample Call: V=HSoptflow(X,3);
14 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
15 [h,w,d,N]=size(Xrgb);
16 if n>N-1
17    error(1,'requested file greater than frame number required');
18 end; 
19 %get two image frames
20 if d==1
21     Xn=double(Xrgb(:,:,1,n));
22     Xnp1=double(Xrgb(:,:,1,n+1));
23 elseif d==3
24     Xn=double(Xrgb(:,:,1,n)*0.299+Xrgb(:,:,2,n)*0.587+Xrgb(:,:,3,n)*0.114);
25     Xnp1=double(Xrgb(:,:,1,n+1)*0.299+Xrgb(:,:,2,n+1)*0.587+Xrgb(:,:,3,n+1)*0.114);
26 else
27     error(2,'not an RGB or Monochrome image file');
28 end;
29 
30 %get image size and adjust for border
31 size(Xn);
32 hm5=h-5; wm5=w-5;
33 z=zeros(h,w); v1=z; v2=z;
34 
35 %initialize
36 dsx2=v1; dsx1=v1; dst=v1;
37 alpha2=625;
38 imax=20;
39 
40 %Calculate gradients
41 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;
42 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;
43 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;
44 
45 for i=1:imax
46    delta=(dsx1.*v1+dsx2.*v2+dst)./(alpha2+dsx1.^2+dsx2.^2);
47    v1=v1-dsx1.*delta;
48    v2=v2-dsx2.*delta;
49 end;
50 u=z; u(5:hm5,5:wm5)=v1(5:hm5,5:wm5);
51 v=z; v(5:hm5,5:wm5)=v2(5:hm5,5:wm5);
52 
53 xskip=round(h/64);
54 [hs,ws]=size(u(1:xskip:h,1:xskip:w))
55 us=zeros(hs,ws); vs=us;
56 
57 N=xskip^2;
58 for i=1:hs-1
59   for j=1:ws-1
60      hk=i*xskip-xskip+1;
61      hl=i*xskip;
62      wk=j*xskip-xskip+1;
63      wl=j*xskip;
64      us(i,j)=sum(sum(u(hk:hl,wk:wl)))/N;
65      vs(i,j)=sum(sum(v(hk:hl,wk:wl)))/N;
66    end;
67 end;
68 
69 figure(1);
70 quiver(us,vs);
71 colormap('default');
72 axis ij;
73 axis tight;
74 axis equal;

 1 close all;
 2 clear all;
 3 clc;
 4 n=40;
 5 
 6 raw=zeros(200,256,1,n);
 7 for i=1:n
 8     filename=strcat('E:\matlab_work\lianxi\data1\1 (',int2str(i),').bmp');
 9     raw(:,:,1,i)=imread(filename);
10 end
11 
12 for i=1:n-1
13     V=HSoptflow(raw,i);
14 end