【目标跟踪】光流法运动视频跟踪【含Matlab源码 1357期】

⛄一、光流场简介

1 案例背景
运动视觉研究的内容是如何从变化场景中的一系列不同时刻的图像中提取有关场景中物体的形状、位置和运动的信息。根据研究的方法，它可以分为两类：基于特征的方法和基于光流场的方法。基于特征的方法抽取特征点，是离散的；光流场属于运动数据研究范畴，是基于连续的图像序列，并直接对其进行运动估计，可以求得图像中每一像素处所对应物体的运动信息。
当物体运动时，在图像上对应物体的亮度模式也在运动。光流(Optical Flow) 是指图像中亮度模式运动的速度，光流场是一种二维瞬时速度场，它是景物中可见点的三维速度矢量在成像表面的投影。光流不仅包含了被观察物体的运动信息，而且携带着有关场景的三维结构信息。本案例基于Computer Vision System Toolbox，使用光流场算法对交通视频中汽车的运动进行检测和估计。

2 理论基础
2.1光流法检测运动原理
光流场是指图像灰度模式的表观运动，它是一种像素级的运动。光流法检测运动物体的基本原理是：根据各个像素点的速度矢量特征，可以对图像进行动态分析。如果图像中没有运动物体，则光流矢量在整个图像区域是连续变化的；当图像中有运动物体时，由于目标和图像背景存在相对运动，所以运动物体所形成的速度矢量必然和邻域背景速度矢量不同，从而检测出运动物体及其位置。但是光流法的优点在于，光流不仅携带了运动物体的运动信息，还携带了有关三维结构的丰富信息，它能够在不知道场景任何信息的情况下，检测出运动的图像。基于光流场的运动检测的步骤如图1所示。在理想情况下，光流场和二维运动场互相吻合，但这一命题不总是对的。如图2所示，一个均匀球体在某一光源照射下，亮度呈现一定的明暗模式。当球体绕中心轴旋转时，明暗模式并不随着表面运动，所以图像也没有变化，此时光流在任意地方都等于零，然而运动场却不等于零。如果球体不动而光源运动，则明暗模式将开始随着光源运动，此时光流不等于零但运动场为零。
在这里插入图片描述
图1基于光流场的运动检测的步骤图2光流场和运动场的联系和区别
光流场的方法能够较好地用于二维运动估计，它也可以同时给出全局点的运动估计，但是光流场并不等价于运动场，因此其本身必然存在着一些问题：遮挡问题、孔径问题、光照问题等。

2.2 光流的主要计算方法
光流场的计算主要有基于梯度的方法、基于匹配的方法、基于能量的方法和基于相位的方法。另外，近几年神经网络动力学也颇受学者重视。基于梯度的方法利用图像灰度的梯度来计算光流，是研究得最多的方法，比如Horn-S chun ck算法、Lucas-Kanade算法和Nagel算法。基于梯度的方法以运动前后图像灰度保持不变作为先决条件，导出光流约束方程。由于光流约束方程并不能唯一地确定光流，因此需要导入其他约束。根据引入的约束不同，基于梯度的方法又可以分为全局约束方法和局部约束方法。全局约束的方法假定光流在整个图像范围内满足一定的约束条件；而局部约束的方法假定在给定点周围的一个小区域内，光流满足一定的约束条件。
基于匹配的方法包括基于特征匹配和基于区域匹配两种。基于区域匹配技术在视频编码中得到了广泛应用，它通过对图像序列中相邻两帧图像间的子块匹配来进行运动估值。在区域匹配算法中，图像被分割为子块，子块中所有像素的运动被认为是相同的，由于复杂的运动可以被近似地分解为一组平移运动之和，所以区域匹配算法采用的运动模型假设图像中的运动物体由做平移运动的刚体组成，且假设在图像场景中没有大的遮挡物。
基于能量的方法首先要对输入图像序列进行时空滤波处理，这是一种时间和空间整合。对于均匀的流场，要获得正确的速度估计，这种时空整合是非常必要的。然而，这样做会降低光流估计的空间和时间分辨率。尤其是当时空整合区域包含几个运动成分(如运动边缘)时，估计精度将会恶化。此外，基于能量的光流技术涉及大量的滤波器，因此存在高计算负荷的问题。然而可以预期，随着相应硬件的发展，在不久的将来，滤波将不再是一个严重的限制因素，所有这些技术都可以在帧速下加以实现。
基于相位的方法由Fleet和Jepson提出，该算法根据带通滤波器输出的相位特性来确定光流。通过与带通速度调谐滤波器输出中的等相位轮廓相垂直的瞬时运动来定义分速度。带通滤波器按照尺度、速度和定向来分离输入信号。基于相位的光流技术的综合性能比较优秀，光流估计比较精确且具有较高的空间分辨率，对图像序列的适用范围也比较宽。
对于光流计算来讲，如果说前面的基于能量或相位的模型有一定的生物合理性的话，那么近几年出现的利用神经网络建立的视觉运动感知的神经动力学模型则是对生物视觉系统功能与结构的更为直接的模拟。尽管现有的神经动力学模型还不成熟，然而这些方法及其结论为其进一步研究打下了良好的基础，是将神经机制引入运动计算方面所做的极有意义的尝试。
目前，对光流的研究方兴未艾，新的计算方法还在不断涌现。这里对光流技术的发展趋势与方向提出以下看法。
(1)现有技术有各自的优点与缺陷，方法之间相互结合，优势互补，建立光流计算的多阶段或分层模型是光流技术发展的一个趋势。
(2)通过深入的研究发现，现有光流方法之间有许多共通之处。如微分法和匹配法的前提假设极为相似；某些基于能量的方法等效于区域匹配技术；而相位方法则将相位梯度用于法向速度的计算。
(3)尽管光流计算的神经动力学方法还很不成熟，然而对它的研究却具有极其深远的意义。随着生物视觉研究的不断深入，神经网络方法无疑会不断完善，也许光流计算乃至计算机视觉的根本出路就在于神经机制的引入。

2.3 梯度光流场约束方程
假定像素点(x，y)在1时刻的灰度值为I(x，y，t)，在1+dr时刻，该像素点运动到新的位置(x+dx，y+dy)，此时对应的灰度值为I(x+dx，y+dy，t+dr)。根据图像的一致性假设，当dt→0时，图像沿着运动轨迹的亮度保持不变，即：I(x，y，t)=I(x+dx，y+dy，t+dt)(23.1)
如果图像灰度随(x，y，t)缓慢变换，则将(23.1)式进行泰勒级数展开：
（1）在这里插入图片描述
如果图像灰度随(x，y，t)缓慢变换，则将(1)式进行泰勒级数展开：

于是

式中，VI-[1…1，]表示梯度方向；U=[u，可表示光流。(23.4)式称为光流约束方程，是
所有基于梯度的光流计算方法的基础。
考虑由u和v组成的二维空间，那么(4)式定义了一条直线，所有满足约束方程的U=[u，都在该直线上，如图3所示，该直线和梯度VI-[1…1，]垂直。由于光流约束方程包含u和v两个未知量，显然由一个方程并不能唯一确定，为了求解光流场，必须引入新的约束条件。
在这里插入图片描述
图3光流基本方程所确定的约束线
根据约束条件的不同，梯度光流法又分为全局约束方法和局部约束方法。全局约束方法假定光流在整个图像范围内满足一定的约束条件，而局部约束的方法假定在给定像素点周围的一个小区域内，光流满足于一定的约束条件。常用的基于梯度的光流计算方法如下。
（1）运动场平滑
Horn-S chun ck假设光流在整个图像上光滑变化，即运动场既满足光流约束方程又满足全局平滑性。将光滑性测度同加权微分约束测度组合起来，其中加权参数控制图像流约束微分和光滑性微分之间的平衡。
（2）预测校正
Lucas-kanade假设在一个小的空间邻域上运动矢量保持恒定，然后使用加权最小二乘的思想来估计光流，它是一种基于像素递归的光流算法，就是预测校正型的位移估算器。预测值可以作为前一个像素位置的运动估算值，或作为当前像素邻域内的运动估算线性组合。依据该像素上的位移帧差的梯度最小值，对预测做进一步的修正。
（3）平滑约束
与Hom-S chun ck算法一样， Nagel也使用了全局平滑约束来建立光流误差测度函数，
但是Nagel提出的一种面向平滑的约束，并不是强加在亮度梯度变化最剧烈的方向(比如
边缘方向)上的，这样做的目的是为了处理遮挡。

2.4 Horn-S chun ck光流算法
Horn-S chun ck算法是一种全局约束的方法，其提出了光流的平滑性约束条件，即图像
上任一点的光流并不是独立的，光流在整个图像范围内平滑变化。所谓平滑，就是在给定
邻域内其速度分量平方和积分最小：
在这里插入图片描述
在实际情况下，(5)式可以使用下面的表达式代替：

式中，u和▼分别表示u邻域和v邻域中的均值。
根据光流基本方程(4) 式考虑光流误差， Horn-S chun ck算法将光流求解归结为如下极值问题：

式中，入控制平滑度，它的取值要考虑图中的噪声情况，如果噪声较强，则说明图像数据
本身的置信度较低，需要更多地依赖光流约束，所以1可以取较大的值；反之，可以取较
小的值。

⛄二、部分源代码

function varargout = guidemo(varargin)
% GUIDEMO M-file for guidemo.fig
% GUIDEMO, by itself, creates a new GUIDEMO or raises the existing
% singleton*.
%
% H = GUIDEMO returns the handle to a new GUIDEMO or the handle to
% the existing singleton*.
%
% GUIDEMO(‘CALLBACK’,hObject,eventData,handles,…) calls the local
% function named CALLBACK in GUIDEMO.M with the given input arguments.
%
% GUIDEMO(‘Property’,‘Value’,…) creates a new GUIDEMO or raises the
% existing singleton*. Starting from the left, property value pairs are
% applied to the GUI before guidemo_OpeningFcn gets called. An
% unrecognized property name or invalid value makes property application
% stop. All inputs are passed to guidemo_OpeningFcn via varargin.
%
% *See GUI Options on GUIDE’s Tools menu. Choose “GUI allows only one
% instance to run (singleton)”.
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help guidemo

% Last Modified by GUIDE v2.5 14-Sep-2021 10:00:52

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct(‘gui_Name’, mfilename, …
‘gui_Singleton’, gui_Singleton, …
‘gui_OpeningFcn’, @guidemo_OpeningFcn, …
‘gui_OutputFcn’, @guidemo_OutputFcn, …
‘gui_LayoutFcn’, [] , …
‘gui_Callback’, []);
if nargin && ischar(varargin{1})
gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% — Executes just before guidemo is made visible.
function guidemo_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% varargin command line arguments to guidemo (see VARARGIN)

% Choose default command line output for guidemo
handles.output = hObject;

a=ones(320,650);
imshow(a);

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes guidemo wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% — Outputs from this function are returned to the command line.
function varargout = guidemo_OutputFcn(hObject, eventdata, handles)
% varargout cell array for returning output args (see VARARGOUT);
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% — Executes on button press in Browse.
function Browse_Callback(hObject, eventdata, handles)
% hObject handle to Browse (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

[filename, pathname] = uigetfile(‘*.avi’, ‘Pick an VIDEO’);
if isequal(filename,0) || isequal(pathname,0)
disp(‘User pressed cancel’)
else
a=aviread(filename);

  handles.a=a;
  
end

% Update handles structure

guidata(hObject, handles);

% — Executes on button press in PlayMovie.
function PlayMovie_Callback(hObject, eventdata, handles)
% hObject handle to PlayMovie (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

B = ‘.png’;
C = ‘yjf-00_1-’;
k=1;

for kk=1:50

    d=num2str(kk);
    e2=strcat(C,d);
    filename3 =strcat(e2,B);
    b=imread( filename3);
    imshow(b);
    pause(0.5);

end
% — Executes on button press in framesepration.
function framesepration_Callback(hObject, eventdata, handles)
% hObject handle to framesepration (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

function BinaryPattern_Callback(hObject, eventdata, handles)
% hObject handle to BinaryPattern (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

% — Executes on button press in Tracking.
function Tracking_Callback(hObject, eventdata, handles)
% hObject handle to Tracking (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

k=1;

for kk=1:50

    d=num2str(kk);
    e2=strcat(C,d);
    filename3 =strcat(e2,B);
    a=imread('yjf-00_1-bk.png');
    b=imread( filename3);
    a1=a(:,:,1);
    a2=a(:,:,2);
    a3=a(:,:,3);
    [r c]=size(a1);
    b1=b(:,:,1);
    b2=b(:,:,2);
    b3=b(:,:,3);
    z1= imsubtract(a1,b1);
    z2= imsubtract(a2,b2);
    z3= imsubtract(a3,b3);

   
    fig(:,:,2)=z2;
    fig(:,:,3)=z3;


    I2=rgb2gray(fig);
    I3 = imadjust(I2, stretchlim(I2), [0 1]);
    level = graythresh(I3);
   
    K6 = medfilt2(K,[5,5]);

% k6=double(K6);
% K6=edge(K6,‘canny’);
% imshow(K6,[])

    k=k+1;

    Image  = K6;
    % get size of image as H,W
    imshow(b);
    [H,W] = size(Image);
    [r,c] = find( bwperim(Image,4) == 1 );
    [tr,tc] = boundarytrack(r,c,H,W,0);
    X1=tr;
  
    YB =tc;
    Xmin = min(XB);
    Ymin = min(YB);
    Xmax = max(XB);
    Ymax = max(YB);
    X = [Xmin Xmax];
    Y = [Ymin Ymin];
    line(Y,X) ;
    X = [Xmin Xmin];
    Y = [Ymin Ymax];
    line(Y,X) ;
    X = [Xmax Xmax];
    Y = [Ymin Ymax];
    line(Y,X) ;
    X = [Xmin Xmax];
    Y = [Ymax Ymax];
    line(Y,X) ;
 
    pause(0.5);

end