1.算法描述
1950年，Gibson首先提出了光流的概念，所谓光流就是指图像表现运动的速度。物体在运动的时候之所以能被人眼发现，就是因为当物体运动时，会在人的视网膜上形成一系列的连续变化的图像，这些变化信息在不同时间，不断的流过眼睛视网膜，就好像一种光流过一样，故称之为光流。光流法检测运动物体的原理：首先给图像中每个像素点赋予一个速度矢量（光流），这样就形成了光流场。如果图像中没有运动物体，光流场连续均匀，如果有运动物体，运动物体的光流和图像的光流不同，光流场不再连续均匀。从而可以检测出运动物体及位置。
应用背景：
根据图像前景和背景的运动，检测视频的变化，空间运动物体在观察成像平面上的像素运动的瞬时速度，是利用图像序列中像素在时间域上的变化以及相邻帧之间的相关性来找到上一帧跟当前帧之间存在的对应关系，从而计算出相邻帧之间物体的运动信息的一种方法。可以用来检测运动抖动物体
关键技术：
当人的眼睛观察运动物体时，物体的景象在人眼的视网膜上形成一系列连续变化的图像，这一系列连续变化的信息不断“流过”视网膜（即图像平面），好像一种光的“流”，故称之为光流（optical flow）。
光流法(Optical flow or optic flow)是关于视域中的物体运动检测中的概念。用来描述相对于观察者的运动所造成的观测目标、表面或边缘的运动。光流法在样型识别、计算机视觉以及其他影像处理领域中非常有用，可用于运动检测、物件切割、碰撞时间与物体膨胀的计算、运动补偿编码，或者通过物体表面与边缘进行立体的测量等等。
光流场，它是指图像中所有像素点构成的一种二维(2D)瞬时速度场，其中的二维速度矢量是景物中可见点的三维速度矢量在成像表面的投影。
所以光流不仅包含了被观察物体的运动信息，而且还包含有关景物三维结构的丰富信息"
下面是针孔相机模型,随着3D点在空间中运动,相应的图像点也在移动.运动场由图像中所有图像点的运动矢量组成
这样,我们回想SLAM14讲的内容:Puv=KPc
假设相机描绘的是动态场景,现在将上式对时间求导,可以得到:
这里
就是我们说的运动场,向量u取决于图像上的2d坐标和时间t.
这里
是相应的3D运动,其与运动场的关系是: 其中,M是一个2*3的矩阵.
这就代表着,对于一个特定的图像点上, 运动场 相对于 位于M的零空间中的3D运动 是不变的.
运动场是理想的构造,描述了2D-3D之间的运动关系.
但实际上,只能基于对图像数据的测量来近似真实的运动场.
问题在于，在大多数情况下，每个图像点都有一个单独的运动，因此必须通过对图像数据的邻域操作来局部测量。结果，无法为某些类型的邻域确定正确的运动场，而是通常被称为光流的近似值
总之,不能正确测量所有像点的运动场，故光流是运动场的近似值。
有几种不同的方法可以根据应如何进行光学估算的不同标准来计算光流。
2.仿真效果预览
matlab2022a仿真结果如下：
3.MATLAB核心程序
fr_f40=rgb2gray(f40);
Xn=double(fr_f1);
Xnp1=double(fr_f40);
%get image size and adjust for border 获取图像对边界进行调整
[h,w]=size(fr_f1);
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;