光流（运动估计）经典算法

happy

光流（运动估计）经典算法Horn and Schunck算法MATLAB程序

1. function [us,vs] = HSoptflow(Xrgb,n)
   
2. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
3. % Author: Gregory Power gregory.power@wpafb.af.mil
   
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. 
    
46. for i=1:imax
    
47.    delta=(dsx1.*v1+dsx2.*v2+dst)./(alpha2+dsx1.^2+dsx2.^2);
    
48.    v1=v1-dsx1.*delta;
    
49.    v2=v2-dsx2.*delta;
    
50. end;
    
51. u=z; u(5:hm5,5:wm5)=v1(5:hm5,5:wm5);
    
52. v=z; v(5:hm5,5:wm5)=v2(5:hm5,5:wm5);
    
53. 
    
54. xskip=round(h/32);
    
55. [hs,ws]=size(u(1:xskip:h,1:xskip:w))
    
56. us=zeros(hs,ws); vs=us;
    
57. 
    
58. N=xskip^2;
    
59. for i=1:hs-1
    
60.   for j=1:ws-1
    
61.      hk=i*xskip-xskip+1;
    
62.      hl=i*xskip;
    
63.      wk=j*xskip-xskip+1;
    
64.      wl=j*xskip;
    
65.      us(i,j)=sum(sum(u(hk:hl,wk:wl)))/N;
    
66.      vs(i,j)=sum(sum(v(hk:hl,wk:wl)))/N;
    
67.    end;
    
68. end;
    
69. 
    
70. figure(1);
    
71. quiver(us,vs);
    
72. colormap('default');
    
73. axis ij;
    
74. axis tight;
    
75. axis equal;

复制代码