Lucas–Kanade光流算法學習

（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 參考程序以下：

%%% 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
    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 %%金字塔擴展
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
    end
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
end

相關標籤/搜索

每日一句

每一个你不满意的现在，都有一个你没有努力的曾经。