Lucas–Kanade光流演算法學習

阿新 • • 發佈：2018-11-10

Lucas–Kanade光流演算法是一種兩幀差分的光流估計演算法。它由Bruce D. Lucas 和 Takeo Kanade提出。

光流(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）位置的體素。
我們假設移動足夠的小，那麼對影象約束方程使用泰勒公式，我們可以得到：

H.O.T. 指更高階，在移動足夠小的情況下可以忽略。從這個方程中我們可以得到：

或者

我們得到：

V x ,V y ,V z 分別是I(x,y,z,t)的光流向量中x，y，z的組成。 , , 和則是影象在(x

,y ,z ,t )這一點向相應方向的差分。
所以

I x V x + I y V y + I z V z = − I t。

寫做：

這個方程有三個未知量，尚不能被解決，這也就是所謂光流演算法的光圈問題。那麼要找到光流向量則需要另一套解決的方案。而Lucas-Kanade演算法是一個非迭代的演算法：

假設流(Vx,Vy,Vz)在一個大小為m*m*m(m>1)的小窗中是一個常數，那麼從畫素1...n , n = m 3 中可以得到下列一組方程：（也就是說，對於這多個點，它們三個方向的速度是一樣的）即，假設了在一個畫素的周圍的一個小的視窗內的所有畫素點的光流大小是相同的。

（

而兩幀影象之間的變化，就是t方向的梯度值，可以理解為當前畫素點沿著光流方向運動而得到的，所以我們可以得到上邊的這個式子。令：

）

三個未知數但是有多於三個的方程，這個方程組自然是個超定方程，也就是說方程組內有冗餘，方程組可以表示為：

記作：

為了解決這個超定問題，我們採用最小二乘法：

得到：

其中的求和是從1到n。

這也就是說尋找光流可以通過在四維上影象導數的分別累加得出。我們還需要一個權重函式W(i, j,k) ，來突出視窗中心點的座標。高斯函式做這項工作是非常合適的，

這個演算法的不足在於它不能產生一個密度很高的流向量，例如在運動的邊緣和黑大的同質區域中的微小移動方面流資訊會很快的褪去。它的優點在於有噪聲存在的魯棒性還是可以的。

補充：opencv裡實現的看上去蠻複雜，現在還不是太明白。其中LK經典演算法也是迭代法，是由高斯迭代法解線性方程組進行迭代的。

參考資料：http://www.cnblogs.com/gnuhpc/archive/2012/12/04/2802124.html

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

另外我對這一部分附上一些個人的看法（謬誤之處還望不吝指正）：

1.首先是假設條件：

（1）亮度恆定，就是同一點隨著時間的變化，其亮度不會發生改變。這是基本光流法的假定（所有光流法變種都必須滿足），用於得到光流法基本方程；

（2）小運動，這個也必須滿足，就是時間的變化不會引起位置的劇烈變化，這樣灰度才能對位置求偏導（換句話說，小運動情況下我們才能用前後幀之間單位位置變化引起的灰度變化去近似灰度對位置的偏導數），這也是光流法不可或缺的假定；

（3）空間一致，一個場景上鄰近的點投影到影象上也是鄰近點，且鄰近點速度一致。這是Lucas-Kanade光流法特有的假定，因為光流法基本方程約束只有一個，而要求x，y方向的速度，有兩個未知變數。我們假定特徵點鄰域內做相似運動，就可以連立n多個方程求取x，y方向的速度（n為特徵點鄰域總點數，包括該特徵點）。

2.方程求解

多個方程求兩個未知變數，又是線性方程，很容易就想到用最小二乘法，事實上opencv也是這麼做的。其中，最小誤差平方和為最優化指標。

3.好吧，前面說到了小運動這個假定，聰明的你肯定很不爽了，目標速度很快那這貨不是二掉了。幸運的是多尺度能解決這個問題。首先，對每一幀建立一個高斯金字塔，最大尺度圖片在最頂層，原始圖片在底層。然後，從頂層開始估計下一幀所在位置，作為下一層的初始位置，沿著金字塔向下搜尋，重複估計動作，直到到達金字塔的底層。聰明的你肯定發現了：這樣搜尋不僅可以解決大運動目標跟蹤，也可以一定程度上解決孔徑問題（相同大小的視窗能覆蓋大尺度圖片上儘量多的角點，而這些角點無法在原始圖片上被覆蓋）。

三.opencv中的光流法函式

opencv2.3.1中已經實現了基於光流法的特徵點位置估計函式（當前幀位置已知，前後幀灰度已知），介紹如下（摘自opencv2.3.1參考手冊）：

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.

四.用類封裝基於光流法的目標跟蹤方法

廢話少說，附上程式碼，包括特徵點提取，跟蹤特徵點，標記特徵點等。

<span style="font-size:18px;">//幀處理基類

class FrameProcessor{

public:

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

};

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

class FeatureTracker : public FrameProcessor{

Mat gray; //當前灰度圖

Mat gray_prev; //之前的灰度圖

vector<Point2f> points[2];//前後兩幀的特徵點

vector<Point2f> initial;//初始特徵點

vector<Point2f> features;//檢測到的特徵

int max_count; //要跟蹤特徵的最大數目

double qlevel; //特徵檢測的指標

double minDist;//特徵點之間最小容忍距離

vector<uchar> status; //特徵點被成功跟蹤的標誌

vector<float> err; //跟蹤時的特徵點小區域誤差和

public:

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

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

//得到灰度圖

cvtColor (frame,gray,CV_BGR2GRAY);

frame.copyTo (output);

//特徵點太少了，重新檢測特徵點

if(addNewPoint()){

detectFeaturePoint ();

//插入檢測到的特徵點

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);

}

//根據前後兩幀灰度圖估計前一幀特徵點在當前幀的位置

//預設視窗是15*15

calcOpticalFlowPyrLK (

gray_prev,//前一幀灰度圖

gray,//當前幀灰度圖

points[0],//前一幀特徵點位置

points[1],//當前幀特徵點位置

status,//特徵點被成功跟蹤的標誌

err);//前一幀特徵點點小區域和當前特徵點小區域間的差，根據差的大小可刪除那些運動變化劇烈的點

int k = 0;

//去除那些未移動的特徵點

for(int i=0;i<points[1].size ();i++){

if(acceptTrackedPoint (i)){

initial[k]=initial[i];

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

}

points[1].resize (k);

initial.resize (k);

//標記被跟蹤的特徵點

handleTrackedPoint (frame,output);

//為下一幀跟蹤初始化特徵點集和灰度影象

轉載自：http://www.xuebuyuan.com/2023445.html

matlab 參考程式如下：

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

%%% 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 %% 執行光流估計

[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;%%%金字塔層數

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

%%% 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 %%金字塔擴充套件

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)%%%濾波係數

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'); %%%計算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

result = interp2 (img, x, y, 'linear');%%%二維資料內插值

%%% Calculates the Fx Fy %% 計算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

好文要頂關注我收藏該文

Lucas–Kanade光流演算法學習

Lucas–Kanade光流演算法學習

【演算法分析】Lucas–Kanade光流演算法

光流法學習

INDEMIND帶你玩轉OpenCV4.0（一）：DIS光流演算法解析

玩轉OpenCV 4.0（一）：呼叫Dis光流演算法示例

光流法學習--optical flow

基於OpenCV的三種光流演算法實現原始碼及測試結果

OpenCV 光流演算法加速---使用GPU來計算光流

Farneback 光流演算法詳解與 calcOpticalFlowFarneback 原始碼分析

OpenCV21（金字塔LK光流演算法）

運動目標跟蹤-光流演算法3

Lucas-Kanade 演算法原理以及應用，正向、反向、additive、Compositional光流法

Python+OpenCV學習（15）---Lucas Kanade 角點光流軌跡跟蹤

【OpenCV筆記】光流法之金字塔Lucas-Kanade

OpenCV學習筆記（二十六）——小試SVM演算法ml OpenCV學習筆記（二十七）——基於級聯分類器的目標檢測objdect OpenCV學習筆記（二十八）——光流法對運動目標跟蹤Video Ope

光流（四）--Kanade-Lucas-Tomasi(KLT)目標跟蹤

OpenCV中光流法跟蹤器的使用方法學習

opencv學習筆記四十二：稠密光流跟蹤

光流 | 基於對偶方法的變分光流改進演算法

光流法-運動目標的檢測（opencv學習）

Lucas–Kanade光流演算法學習

相關推薦