【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

調用系統對象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 <= 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;

 


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM