基于Lucas-Kanade算法的双目图像光流提取matlab仿真

1.算法仿真效果

matlab2022a仿真结果如下：

2.算法涉及理论知识概要

        1950年，Gibson首先提出了光流的概念，所谓光流就是指图像表现运动的速度。物体在运动的时候之所以能被人眼发现，就是因为当物体运动时，会在人的视网膜上形成一系列的连续变化的图像，这些变化信息在不同时间，不断的流过眼睛视网膜，就好像一种光流过一样，故称之为光流。光流法检测运动物体的原理：首先给图像中每个像素点赋予一个速度矢量（光流），这样就形成了光流场。如果图像中没有运动物体，光流场连续均匀，如果有运动物体，运动物体的光流和图像的光流不同，光流场不再连续均匀。从而可以检测出运动物体及位置。

应用背景：

       根据图像前景和背景的运动，检测视频的变化，空间运动物体在观察成像平面上的像素运动的瞬时速度，是利用图像序列中像素在时间域上的变化以及相邻帧之间的相关性来找到上一帧跟当前帧之间存在的对应关系，从而计算出相邻帧之间物体的运动信息的一种方法。可以用来检测运动抖动物体

关键技术：

       当人的眼睛观察运动物体时，物体的景象在人眼的视网膜上形成一系列连续变化的图像，这一系列连续变化的信息不断“流过”视网膜（即图像平面），好像一种光的“流”，故称之为光流（optical flow）。

       光流法(Optical flow or optic flow)是关于视域中的物体运动检测中的概念。用来描述相对于观察者的运动所造成的观测目标、表面或边缘的运动。光流法在样型识别、计算机视觉以及其他影像处理领域中非常有用，可用于运动检测、物件切割、碰撞时间与物体膨胀的计算、运动补偿编码，或者通过物体表面与边缘进行立体的测量等等。

光流场，它是指图像中所有像素点构成的一种二维(2D)瞬时速度场，其中的二维速度矢量是景物中可见点的三维速度矢量在成像表面的投影。

所以光流不仅包含了被观察物体的运动信息，而且还包含有关景物三维结构的丰富信息"

光流的概念:(Optical flow or optic flow)

       它是一种运动模式，这种运动模式指的是一个物体、表面、边缘在一个视角下由一个观察者（比如眼睛、摄像头等）和背景之间形成的明显移动。光流技术，如运动检测和图像分割，时间碰撞，运动补偿编码，三维立体视差，都是利用了这种边缘或表面运动的技术。

二维图像的移动相对于观察者而言是三维物体移动的在图像平面的投影。

有序的图像可以估计出二维图像的瞬时图像速率或离散图像转移。

光流算法：

      它评估了两幅图像的之间的变形，它的基本假设是体素和图像像素守恒。它假设一个物体的颜色在前后两帧没有巨大而明显的变化。基于这个思路，我们可以得到图像约束方程。不同的光流算法解决了假定了不同附加条件的光流问题。

Lucas–Kanade算法：

      这个算法是最常见，最流行的。它计算两帧在时间t 到t + δt之间每个每个像素点位置的移动。 由于它是基于图像信号的泰勒级数，这种方法称为差分，这就是对于空间和时间坐标使用偏导数。

图像约束方程可以写为I (x ,y ,z ,t ) = I (x + δx ,y + δy ,z + δz ,t + δt )

I(x, y,z, t) 为在（x,y,z）位置的体素。

我们假设移动足够的小，那么对图像约束方程使用泰勒公式，我们可以得到：

3.MATLAB核心程序

if size(size(im1),2)==3
    im1=rgb2gray(im1);
end
if size(size(im2),2)==3
    im2=rgb2gray(im2);
end
 
im1=double(im1);
im2=double(im2);
 
im1 = imgaussfilt(im1,2);
im2 = imgaussfilt(im2,2);
 
 
 
u = zeros(size(im1)); % create u and v matrices
v = zeros(size(im2));
 
 
IxMask = [-1 1; -1 1]; % Create the X derivative filter
IyMask = [-1 -1; 1 1]; % Create the Y derivative filter
 
ItMask_1 = [-1 -1; -1 -1]; %Create both the filters for the time derivative
ItMask_2 = [1 1; 1 1];
 
im1x = conv2(im1, IxMask,'same');
im2x = conv2(im2, IxMask,'same');
im1y = conv2(im1, IyMask,'same');
im2y = conv2(im2, IyMask,'same');
 
im1t = conv2(im1, ItMask_1,'same');
im2t = conv2(im2, ItMask_2,'same');
 
Ix = (im1x + im2x) ./ 2;
Iy = (im1y + im2y) ./ 2;
It = (im1t + im2t) ./ 2;
 
newnumb1 = size(im1,1) - mod(size(im1,1),windowSize);
newnumb2 = size(im1,2) - mod(size(im1,2),windowSize);
 
for i = 1 : newnumb1 - windowSize
    for j = 1 : newnumb2 - windowSize
        
        Ixwindow=Ix(i  : i + windowSize , j : j + windowSize);
        Iywindow=Iy(i  : i + windowSize , j : j + windowSize);
        Itwindow=It(i  : i + windowSize , j : j + windowSize);
        
        Ixwindow = reshape(Ixwindow , [1 , size(Ixwindow,1)^2 ])';
        Iywindow = reshape(Iywindow , [1 , size(Iywindow,1)^2 ])';
        Itwindow = reshape(Itwindow , [1 , size(Itwindow,1)^2 ])';
        
        A=[Ixwindow Iywindow];
        B=Itwindow;
        [outpute]=inv(A'*A)*A'*B;
        
        u(i  , j) = outpute(1);
        v(i  , j) = outpute(2);
        
        
    end
end
u(isnan(u))=0;
v(isnan(v))=0;