(

而兩幀圖像之間的變化，就是t方向的梯度值，可以理解為當(dāng)前像素點(diǎn)沿著光流方向運(yùn)動(dòng)而得到的，所以我們可以得到上邊的這個(gè)式子。令：

)

這一部分《learing opencv》一書的第10章Lucas-Kanade光流部分寫得非常詳細(xì)，推薦大家看書。

另外我對(duì)這一部分附上一些個(gè)人的看法(謬誤之處還望不吝指正)：

1.首先是假設(shè)條件：

(1)亮度恒定，就是同一點(diǎn)隨著時(shí)間的變化，其亮度不會(huì)發(fā)生改變。這是基本光流法的假定(所有光流法變種都必須滿足)，用于得到光流法基本方程；

(2)小運(yùn)動(dòng)，這個(gè)也必須滿足，就是時(shí)間的變化不會(huì)引起位置的劇烈變化，這樣灰度才能對(duì)位置求偏導(dǎo)(換句話說，小運(yùn)動(dòng)情況下我們才能用前后幀之間單位位置變化引起的灰度變化去近似灰度對(duì)位置的偏導(dǎo)數(shù))，這也是光流法不可或缺的假定；

(3)空間一致，一個(gè)場(chǎng)景上鄰近的點(diǎn)投影到圖像上也是鄰近點(diǎn)，且鄰近點(diǎn)速度一致。這是Lucas-Kanade光流法特有的假定，因?yàn)楣饬鞣ɑ痉匠碳s束只有一個(gè)，而要求x，y方向的速度，有兩個(gè)未知變量。我們假定特征點(diǎn)鄰域內(nèi)做相似運(yùn)動(dòng)，就可以連立n多個(gè)方程求取x，y方向的速度(n為特征點(diǎn)鄰域總點(diǎn)數(shù)，包括該特征點(diǎn))。

2.方程求解

多個(gè)方程求兩個(gè)未知變量，又是線性方程，很容易就想到用最小二乘法，事實(shí)上opencv也是這么做的。其中，最小誤差平方和為最優(yōu)化指標(biāo)。

3.好吧，前面說到了小運(yùn)動(dòng)這個(gè)假定，聰明的你肯定很不爽了，目標(biāo)速度很快那這貨不是二掉了。幸運(yùn)的是多尺度能解決這個(gè)問題。首先，對(duì)每一幀建立一個(gè)高斯金字塔，最大尺度圖片在最頂層，原始圖片在底層。然后，從頂層開始估計(jì)下一幀所在位置，作為下一層的初始位置，沿著金字塔向下搜索，重復(fù)估計(jì)動(dòng)作，直到到達(dá)金字塔的底層。聰明的你肯定發(fā)現(xiàn)了：這樣搜索不僅可以解決大運(yùn)動(dòng)目標(biāo)跟蹤，也可以一定程度上解決孔徑問題(相同大小的窗口能覆蓋大尺度圖片上盡量多的角點(diǎn)，而這些角點(diǎn)無(wú)法在原始圖片上被覆蓋)。

三.opencv中的光流法函數(shù)

opencv2.3.1中已經(jīng)實(shí)現(xiàn)了基于光流法的特征點(diǎn)位置估計(jì)函數(shù)(當(dāng)前幀位置已知，前后幀灰度已知)，介紹如下(摘自opencv2.3.1參考手冊(cè))：calcOpticalFlowPyrLK

Calculates an optical flow for a sparse feature set using the iterative Lucas-Kanade method with pyramids.

void calcOpticalFlowPyrLK(InputArray prevImg, InputArray nextImg, InputArray prevPts,

InputOutputArray nextPts, OutputArray status, OutputArray err,

Size winSize=Size(15,15), int maxLevel=3, TermCriteria crite-

ria=TermCriteria(TermCriteria::COUNT+TermCriteria::EPS, 30, 0.01),

double derivLambda=0.5, int flags=0 )

Parameters

prevImg – First 8-bit single-channel or 3-channel input image.

nextImg – Second input image of the same size and the same type as prevImg .

prevPts – Vector of 2D points for which the flow needs to be found. The point coordinates

must be single-precision floating-point numbers.

nextPts – Output vector of 2D points (with single-precision floating-point coordinates)

containing the calculated new positions of input features in the second image. When

OPTFLOW_USE_INITIAL_FLOW flag is passed, the vector must have the same size as in the

input.

status – Output status vector. Each element of the vector is set to 1 if the flow for the

corresponding features has been found. Otherwise, it is set to 0.

err – Output vector that contains the difference between patches around the original and

moved points.

winSize – Size of the search window at each pyramid level.

maxLevel – 0-based maximal pyramid level number. If set to 0, pyramids are not used

(single level). If set to 1, two levels are used, and so on.

criteria – Parameter specifying the termination criteria of the iterative search algorithm

(after the specified maximum number of iterations criteria.maxCount or when the search

window moves by less than criteria.epsilon .

derivLambda – Not used.

flags – Operation flags:

– OPTFLOW_USE_INITIAL_FLOW Use initial estimations stored in nextPts . If the

flag is not set, then prevPts is copied to nextPts and is considered as the initial estimate.

四.用類封裝基于光流法的目標(biāo)跟蹤方法

廢話少說，附上代碼，包括特征點(diǎn)提取，跟蹤特征點(diǎn)，標(biāo)記特征點(diǎn)等。//幀處理基類

class FrameProcessor{

public:

virtual void process(Mat &input,Mat &ouput)=0;

};

//特征跟蹤類，繼承自幀處理基類

class FeatureTracker : public FrameProcessor{

Mat gray; //當(dāng)前灰度圖

Mat gray_prev; //之前的灰度圖

vector points[2];//前后兩幀的特征點(diǎn)

vector initial;//初始特征點(diǎn)

vector features;//檢測(cè)到的特征

int max_count; //要跟蹤特征的最大數(shù)目

double qlevel; //特征檢測(cè)的指標(biāo)

double minDist;//特征點(diǎn)之間最小容忍距離

vector status; //特征點(diǎn)被成功跟蹤的標(biāo)志

vector err; //跟蹤時(shí)的特征點(diǎn)小區(qū)域誤差和

public:

FeatureTracker():max_count(500),qlevel(0.01),minDist(10.){}

void process(Mat &frame,Mat &output){

//得到灰度圖

cvtColor (frame,gray,CV_BGR2GRAY);

frame.copyTo (output);

//特征點(diǎn)太少了，重新檢測(cè)特征點(diǎn)

if(addNewPoint()){

detectFeaturePoint ();

//插入檢測(cè)到的特征點(diǎn)

points[0].insert (points[0].end (),features.begin (),features.end ());

initial.insert (initial.end (),features.begin (),features.end ());

}

//第一幀

if(gray_prev.empty ()){

gray.copyTo (gray_prev);

}

//根據(jù)前后兩幀灰度圖估計(jì)前一幀特征點(diǎn)在當(dāng)前幀的位置

//默認(rèn)窗口是15*15

calcOpticalFlowPyrLK (

gray_prev,//前一幀灰度圖

gray,//當(dāng)前幀灰度圖

points[0],//前一幀特征點(diǎn)位置

points[1],//當(dāng)前幀特征點(diǎn)位置

status,//特征點(diǎn)被成功跟蹤的標(biāo)志

err);//前一幀特征點(diǎn)點(diǎn)小區(qū)域和當(dāng)前特征點(diǎn)小區(qū)域間的差，根據(jù)差的大小可刪除那些運(yùn)動(dòng)變化劇烈的點(diǎn)

int k = 0;

//去除那些未移動(dòng)的特征點(diǎn)

for(int i=0;i

if(acceptTrackedPoint (i)){

initial[k]=initial[i];

points[1][k++] = points[1][i];

}

}

points[1].resize (k);

initial.resize (k);

//標(biāo)記被跟蹤的特征點(diǎn)

handleTrackedPoint (frame,output);

//為下一幀跟蹤初始化特征點(diǎn)集和灰度圖像

matlab 參考程序 如下：%%% Usage: Lucas_Kanade('1.bmp','2.bmp',10)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%file1：輸入圖像1

%%%file2：輸入圖像2

%%%density：顯示密度

function Lucas_Kanade(file1, file2, density)

%%% Read Images %% 讀取圖像

img1 = im2double (imread (file1));

%%% Take alternating rows and columns %% 按行分成奇偶

[odd1, ~] = split (img1);

img2 = im2double (imread (file2));

[odd2, ~] = split (img2);

%%% Run Lucas Kanade %% 運(yùn)行光流估計(jì)

[Dx, Dy] = Estimate (odd1, odd2);

%%% Plot %% 繪圖

figure;

[maxI,maxJ] = size(Dx);

Dx = Dx(1:density:maxI,1:density:maxJ);

Dy = Dy(1:density:maxI,1:density:maxJ);

quiver(1:density:maxJ, (maxI):(-density):1, Dx, -Dy, 1);

axis square;

%%% Run Lucas Kanade on all levels and interpolate %% 光流

function [Dx, Dy] = Estimate (img1, img2)

level = 4;%%%金字塔層數(shù)

half_window_size = 2;

% [m, n] = size (img1);

G00 = img1;

G10 = img2;

if (level > 0)%%%從零到level

G01 = reduce (G00); G11 = reduce (G10);

end

if (level>1)

G02 = reduce (G01); G12 = reduce (G11);

end

if (level>2)

G03 = reduce (G02); G13 = reduce (G12);

end

if (level>3)

G04 = reduce (G03); G14 = reduce (G13);

end

l = level;

for i = level: -1 :0,

if (l == level)

switch (l)

case 4, Dx = zeros (size (G04)); Dy = zeros (size (G04));

case 3, Dx = zeros (size (G03)); Dy = zeros (size (G03));

case 2, Dx = zeros (size (G02)); Dy = zeros (size (G02));

case 1, Dx = zeros (size (G01)); Dy = zeros (size (G01));

case 0, Dx = zeros (size (G00)); Dy = zeros (size (G00));

end

else

Dx = expand (Dx);

Dy = expand (Dy);

Dx = Dx .* 2;

Dy = Dy .* 2;

end

switch (l)

case 4,

W = warp (G04, Dx, Dy);

[Vx, Vy] = EstimateMotion (W, G14, half_window_size);

case 3,

W = warp (G03, Dx, Dy);

[Vx, Vy] = EstimateMotion (W, G13, half_window_size);

case 2,

W = warp (G02, Dx, Dy);

[Vx, Vy] = EstimateMotion (W, G12, half_window_size);

case 1,

W = warp (G01, Dx, Dy);

[Vx, Vy] = EstimateMotion (W, G11, half_window_size);

case 0,

W = warp (G00, Dx, Dy);

[Vx, Vy] = EstimateMotion (W, G10, half_window_size);

end

[m, n] = size (W);

Dx(1:m, 1:n) = Dx(1:m,1:n) + Vx; Dy(1:m, 1:n) = Dy(1:m, 1:n) + Vy;

smooth (Dx);

smooth (Dy);

l = l - 1;

end

%%% Lucas Kanade on the image sequence at pyramid step %%

function [Vx, Vy] = EstimateMotion (W, G1, half_window_size)

[m, n] = size (W);

Vx = zeros (size (W)); Vy = zeros (size (W));

N = zeros (2*half_window_size+1, 5);

for i = 1:m,

l = 0;

for j = 1-half_window_size:1+half_window_size,

l = l + 1;

N (l,:) = getSlice (W, G1, i, j, half_window_size);

end

replace = 1;

for j = 1:n,

t = sum (N);

[v, d] = eig ([t(1) t(2);t(2) t(3)]);

namda1 = d(1,1); namda2 = d(2,2);

if (namda1 > namda2)

tmp = namda1; namda1 = namda2; namda2 = tmp;

tmp1 = v (:,1); v(:,1) = v(:,2); v(:,2) = tmp1;

end

if (namda2 < 0.001)

Vx (i, j) = 0; Vy (i, j) = 0;

elseif (namda2 > 100 * namda1)

n2 = v(1,2) * t(4) + v(2,2) * t(5);

Vx (i,j) = n2 * v(1,2) / namda2;

Vy (i,j) = n2 * v(2,2) / namda2;

else

n1 = v(1,1) * t(4) + v(2,1) * t(5);

n2 = v(1,2) * t(4) + v(2,2) * t(5);

Vx (i,j) = n1 * v(1,1) / namda1 + n2 * v(1,2) / namda2;

Vy (i,j) = n1 * v(2,1) / namda1 + n2 * v(2,2) / namda2;

end

N (replace, :) = getSlice (W, G1, i, j+half_window_size+1, half_window_size);

replace = replace + 1;

if (replace == 2 * half_window_size + 2)

replace = 1;

end

end

end

%%% The Reduce Function for pyramid %%金字塔壓縮

function result = reduce (ori)

[m,n] = size (ori);

mid = zeros (m, n);

%%%行列值的一半

m1 = round (m/2);

n1 = round (n/2);

%%%輸出即為輸入圖像的1/4

result = zeros (m1, n1);

%%%濾波

%%%0.05 0.25 0.40 0.25 0.05

w = generateFilter (0.4);

for j = 1 : m,

tmp = conv([ori(j,n-1:n) ori(j,1:n) ori(j,1:2)], w);

mid(j,1:n1) = tmp(5:2:n+4);

end

for i=1:n1,

tmp = conv([mid(m-1:m,i); mid(1:m,i); mid(1:2,i)]', w);

result(1:m1,i) = tmp(5:2:m+4)';

end

%%% The Expansion Function for pyramid %%金字塔擴(kuò)展

function result = expand (ori)

[m,n] = size (ori);

mid = zeros (m, n);

%%%行列值的兩倍

m1 = m * 2;

n1 = n * 2;

%%%輸出即為輸入圖像的4倍

result = zeros (m1, n1);

%%%濾波

%%%0.05 0.25 0.40 0.25 0.05

w = generateFilter (0.4);

for j=1:m,

t = zeros (1, n1);

t(1:2:n1-1) = ori (j,1:n);

tmp = conv ([ori(j,n) 0 t ori(j,1) 0], w);

mid(j,1:n1) = 2 .* tmp (5:n1+4);

end

for i=1:n1,

t = zeros (1, m1);

t(1:2:m1-1) = mid (1:m,i)';

tmp = conv([mid(m,i) 0 t mid(1,i) 0], w);

result(1:m1,i) = 2 .* tmp (5:m1+4)';

end

function filter = generateFilter (alpha)%%%濾波系數(shù)

filter = [0.25-alpha/2, 0.25, alpha, 0.25, 0.25-alpha/2];

function [N] = getSlice (W, G1, i, j, half_window_size)

N = zeros (1, 5);

[m, n] = size (W);

for y = -half_window_size:half_window_size,

Y1 = y +i;

if (Y1 < 1)

Y1 = Y1 + m;

elseif (Y1 > m)

Y1 = Y1 - m;

end

X1 = j;

if (X1 < 1)

X1 = X1 + n;

elseif (X1 > n)

X1 = X1 - n;

end

DeriX = Derivative (G1, X1, Y1, 'x'); %%%計(jì)算x、y方向梯度

DeriY = Derivative (G1, X1, Y1, 'y');

N = N + [ DeriX * DeriX, ...

DeriX * DeriY, ...

DeriY * DeriY, ...

DeriX * (G1 (Y1, X1) - W (Y1, X1)), ...

DeriY * (G1 (Y1, X1) - W (Y1, X1))];

end

function result = smooth (img)

result = expand (reduce (img));%%%太碉

function [odd, even] = split (img)

[m, ~] = size (img);%%%按行分成奇偶

odd = img (1:2:m, :);

even = img (2:2:m, :);

%%% Interpolation %% 插值

function result = warp (img, Dx, Dy)

[m, n] = size (img);

[x,y] = meshgrid (1:n, 1:m);

x = x + Dx (1:m, 1:n);

y = y + Dy (1:m,1:n);

for i=1:m,

for j=1:n,

if x(i,j)>n

x(i,j) = n;

end

if x(i,j)<1

x(i,j) = 1;

end

if y(i,j)>m

y(i,j) = m;

end

if y(i,j)<1

y(i,j) = 1;

end

end

end

result = interp2 (img, x, y, 'linear');%%%二維數(shù)據(jù)內(nèi)插值

%%% Calculates the Fx Fy %% 計(jì)算x、y方向梯度

function result = Derivative (img, x, y, direction)

[m, n] = size (img);

switch (direction)

case 'x',

if (x == 1)

result = img (y, x+1) - img (y, x);

elseif (x == n)

result = img (y, x) - img (y, x-1);

else

result = 0.5 * (img (y, x+1) - img (y, x-1));

end

case 'y',

if (y == 1)

result = img (y+1, x) - img (y, x);

elseif (y == m)

result = img (y, x) - img (y-1, x);

else

result = 0.5 * (img (y+1, x) - img (y-1, x));

end

end

總結(jié)

以上是生活随笔為你收集整理的lucas–kanade_Lucas–Kanade光流算法学习的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。