【Matlab】运动目标检测之“光流法”

光流（optical flow）

1950年，Gibson首先提出了光流的概念，所谓光流就是指图像表现运动的速度。物体在运动的时候之所以能被人眼发现，就是因为当物体运动时，会在人的视网膜上形成一系列的连续变化的图像，这些变化信息在不同时间，不断的流过眼睛视网膜，就好像一种光流过一样，故称之为光流。

光流法检测运动物体的原理：首先给图像中每个像素点赋予一个速度矢量（光流），这样就形成了光流场。如果图像中没有运动物体，光流场连续均匀，如果有运动物体，运动物体的光流和图像的光流不同，光流场不再连续均匀。从而可以检测出运动物体及位置。

应用背景：

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

关键技术：

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

编程处理中：

matlab中有现成的！！函数

function [fx, fy, ft] = computeDerivatives(im1, im2)

if size(im2,1)==0
    im2=zeros(size(im1));
end

% Horn-Schunck original method
fx = conv2(im1,0.25* [-1 1; -1 1],\'same\') + conv2(im2, 0.25*[-1 1; -1 1],\'same\');
fy = conv2(im1, 0.25*[-1 -1; 1 1], \'same\') + conv2(im2, 0.25*[-1 -1; 1 1], \'same\');
ft = conv2(im1, 0.25*ones(2),\'same\') + conv2(im2, -0.25*ones(2),\'same\');

% derivatives as in Barron
% fx= conv2(im1,(1/12)*[-1 8 0 -8 1],\'same\');
% fy= conv2(im1,(1/12)*[-1 8 0 -8 1]\',\'same\');
% ft = conv2(im1, 0.25*ones(2),\'same\') + conv2(im2, -0.25*ones(2),\'same\');
% fx=-fx;fy=-fy;

% An alternative way to compute the spatiotemporal derivatives is to use simple finite difference masks.
% fx = conv2(im1,[1 -1]);
% fy = conv2(im1,[1; -1]);
% ft= im2-im1;

也有现成的实例：

Affine optic flow - File Exchange - MATLAB Central https://www.mathworks.com/matlabcentral/fileexchange/27093-affine-optic-flow

Estimate optical flow using Horn-Schunck method - MATLAB http://www.mathworks.com/help/vision/ref/opticalflowhs-class.html

调用系统对象vision.OpticalFlow后产生的混合矩阵数据如何处理 – MATLAB中文论坛 http://www.ilovematlab.cn/thread-292544-1-1.html

Estimate optical flow using Lucas-Kanade method - MATLAB http://www.mathworks.com/help/vision/ref/opticalflowlk-class.html

Lucas-Kanade Tutorial Example 1 - File Exchange - MATLAB Central https://www.mathworks.com/matlabcentral/fileexchange/48744-lucas-kanade-tutorial-example-1?focused=3854179&tab=example

1.首先是假设条件：

（1）亮度恒定，就是同一点随着时间的变化，其亮度不会发生改变。这是基本光流法的假定（所有光流法变种都必须满足），用于得到光流法基本方程；

（2）小运动，这个也必须满足，就是时间的变化不会引起位置的剧烈变化，这样灰度才能对位置求偏导（换句话说，小运动情况下我们才能用前后帧之间单位位置变化引起的灰度变化去近似灰度对位置的偏导数），这也是光流法不可或缺的假定；

（3）空间一致，一个场景上邻近的点投影到图像上也是邻近点，且邻近点速度一致。这是Lucas-Kanade光流法特有的假定，因为光流法基本方程约束只有一个，而要求x，y方向的速度，有两个未知变量。我们假定特征点邻域内做相似运动，就可以连立n多个方程求取x，y方向的速度（n为特征点邻域总点数，包括该特征点）。

2.方程求解

多个方程求两个未知变量，又是线性方程，很容易就想到用最小二乘法，事实上opencv也是这么做的。其中，最小误差平方和为最优化指标。

3.前面说到了小运动这个假定，目标速度很快的话用多尺度能解决这个问题。首先，对每一帧建立一个高斯金字塔，最大尺度图片在最顶层，原始图片在底层。然后，从顶层开始估计下一帧所在位置，作为下一层的初始位置，沿着金字塔向下搜索，重复估计动作，直到到达金字塔的底层。聪明的你肯定发现了：这样搜索不仅可以解决大运动目标跟踪，也可以一定程度上解决孔径问题（相同大小的窗口能覆盖大尺度图片上尽量多的角点，而这些角点无法在原始图片上被覆盖）。

在计算机视觉中，Lucas–Kanade光流算法是一种两帧差分的光流估计算法。它由Bruce D. Lucas 和 Takeo Kanade提出。

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

二维图像的移动相对于观察者而言是三维物体移动的在图像平面的投影。有序的图像可以估计出二维图像的瞬时图像速率或离散图像转移。

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

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

目标跟踪:

光流法是用于目标跟踪：对于一个连续的视频帧序列进行处理；针对每一个视频序列，利用一定目标检测方法，检测可能出现的前景目标；如果某一帧出现了前景目标，找到其具有代表性的关键特征点；对之后的任意两个相邻视频而言，寻找上一帧中出现的关键特征点在当前帧中的最佳位置，从而得到前景目标在当前帧中的位置坐标；如此迭代进行，便可实现目标的跟踪。

Search Results: Optical Flow type:"Example" - File Exchange - MATLAB Central https://www.mathworks.com/matlabcentral/fileexchange/?term=Optical+Flow+type%3A%22Example%22

一些代码：用帮助文档中vision.OpticalFlow中的例子

Optical Flow using Matlab - Free Open Source Codes - CodeForge.com http://www.codeforge.com/article/237408

【Matlab】运动目标检测之“光流法”

调用系统对象vision.OpticalFlow后产生的混合矩阵数据如何处理 – MATLAB中文论坛 http://www.ilovematlab.cn/thread-292544-1-1.html

% 利用系统对象 L-K光流法检测运动
%
clc;clear all; close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 创建各类系统对象
%用于读入视频的系统对象；
hvfr= vision.VideoFileReader(\'SampleVideo.avi\',...
    \'ImageColorSpace\',\'Intensity\',\'VideoOutputDataType\',\'uint8\');

% 用于图像数据类型转换的系统对象；
hidtc= vision.ImageDataTypeConverter;

%用于光流法检测的系统对象；
hof=vision.OpticalFlow(\'ReferenceFrameDelay\',1);
hof.OutputValue=\'Horizontal and vertical components in complex form\';

%用于在图像中绘制标记;
hsi=vision.ShapeInserter(\'Shape\',\'Lines\',\'BorderColor\',\'Custom\',...
    \'CustomBorderColor\',255);

%用于播放视频图像的系统对象；
hvp = vision.VideoPlayer (\'Name\',\'Motion Vector\');

while ~ isDone(hvfr)
    frame=step (hvfr);  %读入视频
    im = step (hidtc,frame); %将图像的每帧视频图像转换成单精度型
    of = step (hof,im); %用光流法对视频中的每一帧图像进行处理
    lines=videooptflowlines(of,20); %产生坐标点
    if ~ isempty(lines)
        out = step(hsi,im,lines);  %标记出光流
        step(hvp,out);  %观看检测效果
    end
end

%释放系统对象
release(hvp);
release(hvfr);

图示识别光流Optical Flow - CSDN博客 https://blog.csdn.net/jinxinsummer/article/details/77345035

clc
clear
close all
% 读取文件对象
mp4_name=\'Mp4\sub01\sub01_01.mp4\';
videoReader = vision.VideoFileReader(mp4_name,\'ImageColorSpace\',\'Intensity\',\'VideoOutputDataType\',\'uint8\');
% 类型转化对象
converter = vision.ImageDataTypeConverter;
% 光流对象
opticalFlow = vision.OpticalFlow(\'ReferenceFrameDelay\', 1);
opticalFlow.OutputValue = \'Horizontal and vertical components in complex form\';
if 0 % 使用的算法
opticalFlow.Method = \'Lucas-Kanade\';
opticalFlow.NoiseReductionThreshold = 0.001; % 默认是0.0039
else
opticalFlow.Method = \'Horn-Schunck\';
opticalFlow.Smoothness = 0.5; % 默认是1
end
% 显示对象
frame = step(videoReader);
figure
subplot(121)
himg = imshow(frame);
subplot(122)
hof = imshow(frame);
 
%用于在图像中绘制标记;
hsi=vision.ShapeInserter(\'Shape\',\'Lines\',\'BorderColor\',\'Custom\',...
 \'CustomBorderColor\',255);
 
%用于播放视频图像的系统对象；
hvp = vision.VideoPlayer (\'Name\',\'Motion Vector\');
 
% 开始播放
ii=1;
while ~isDone(videoReader)
% 得到一帧
frame = step(videoReader);
% 格式转化
im = step(converter, frame);
% 计算光流
of = step(opticalFlow, im);
%产生坐标点
lines=videooptflowlines(of,20); 
   if ~ isempty(lines)
   out = step(hsi,im,lines); %标记出光流
   step(hvp,out);%观看检测效果
   end
 
% 光流图转化
ofI = computeColor(real(of), imag(of));
path=strcat(\'opticalFlow_img\sub01\01\\',num2str(ii),\'.jpg\');
imwrite(ofI,path);
ii=ii+1;
% figure, imshow(ofI);
% 显示
set(himg, \'cdata\', frame)
set(hof, \'cdata\', ofI)
drawnow
end
release(videoReader);
release(hvp);

其中，computerColor()函数如下：

function img = computeColor(u,v)
 
%   computeColor color codes flow field U, V
 
%   According to the c++ source code of Daniel Scharstein 
%   Contact: schar@middlebury.edu
 
%   Author: Deqing Sun, Department of Computer Science, Brown University
%   Contact: dqsun@cs.brown.edu
%   $Date: 2007-10-31 21:20:30 (Wed, 31 Oct 2006) $
 
% Copyright 2007, Deqing Sun.
%
%                         All Rights Reserved
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for any purpose other than its incorporation into a
% commercial product is hereby granted without fee, provided that the
% above copyright notice appear in all copies and that both that
% copyright notice and this permission notice appear in supporting
% documentation, and that the name of the author and Brown University not be used in
% advertising or publicity pertaining to distribution of the software
% without specific, written prior permission.
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
% INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR ANY
% PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR BROWN UNIVERSITY BE LIABLE FOR
% ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
% WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
% ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
% OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. 
 
nanIdx = isnan(u) | isnan(v);
u(nanIdx) = 0;
v(nanIdx) = 0;
 
colorwheel = makeColorwheel();
ncols = size(colorwheel, 1);
 
rad = sqrt(u.^2+v.^2);          
 
a = atan2(-v, -u)/pi;
 
fk = (a+1) /2 * (ncols-1) + 1;  % -1~1 maped to 1~ncols
   
k0 = floor(fk);                 % 1, 2, ..., ncols
 
k1 = k0+1;
k1(k1==ncols+1) = 1;
 
f = fk - k0;
 
for i = 1:size(colorwheel,2)
    tmp = colorwheel(:,i);
    col0 = tmp(k0)/255;
    col1 = tmp(k1)/255;
    col = (1-f).*col0 + f.*col1;   
   
    idx = rad &lt;= 1;   
    col(idx) = 1-rad(idx).*(1-col(idx));    % increase saturation with radius
    
    col(~idx) = col(~idx)*0.75;             % out of range
    
    img(:,:, i) = uint8(floor(255*col.*(1-nanIdx)));         
end;    
 
%%
function colorwheel = makeColorwheel()
 
%   color encoding scheme
 
%   adapted from the color circle idea described at
%   http://members.shaw.ca/quadibloc/other/colint.htm
 
 
RY = 15;
YG = 6;
GC = 4;
CB = 11;
BM = 13;
MR = 6;
 
ncols = RY + YG + GC + CB + BM + MR;
 
colorwheel = zeros(ncols, 3); % r g b
 
col = 0;
%RY
colorwheel(1:RY, 1) = 255;
colorwheel(1:RY, 2) = floor(255*(0:RY-1)/RY)\';
col = col+RY;
%YG
colorwheel(col+(1:YG), 1) = 255 - floor(255*(0:YG-1)/YG)\';
colorwheel(col+(1:YG), 2) = 255;
col = col+YG;
 
%GC
colorwheel(col+(1:GC), 2) = 255;
colorwheel(col+(1:GC), 3) = floor(255*(0:GC-1)/GC)\';
col = col+GC;
%CB
colorwheel(col+(1:CB), 2) = 255 - floor(255*(0:CB-1)/CB)\';
colorwheel(col+(1:CB), 3) = 255;
col = col+CB;
 
%BM
colorwheel(col+(1:BM), 3) = 255;
colorwheel(col+(1:BM), 1) = floor(255*(0:BM-1)/BM)\';
col = col+BM;
%MR
colorwheel(col+(1:MR), 3) = 255 - floor(255*(0:MR-1)/MR)\';
colorwheel(col+(1:MR), 1) = 255;