光流法optical flow

标签：光流法

光流是图像亮度的运动信息描述。光流法计算最初是由Horn和Schunck于1981年提出的，创造性地将二维速度场与灰度相联系，引入光流约束方程，得到光流计算的基本算法．光流计算基于物体移动的光学特性提出了2个假设：

①运动物体的灰度在很短的间隔时间内保持不变；
②给定邻域内的速度向量场变化是缓慢的。

算法原理

假设图像上一个像素点(x,y)，在t时刻的亮度为E(x+Δx,y+Δy,t+Δt)，同时用u(x,y0和v(x,y)来表示该点光流在水平和垂直方向上的移动分量：

u=dx/dt

v=dy/dt

在经过一段时间间隔Δt后该点对应点亮度为E(x+Δx,y+Δy,t+Δt),当Δt很小趋近于0时，我们可以认为该点亮度不变，所以可以有：
E(x,y,t)=E(x+Δx,y+Δy,t+Δt)
当该点的亮度有变化时，将移动后点的亮度由Taylor公式展幵，可得：

忽略其二阶无穷小，由于Δt趋近于0时，有：

式中w=(u,v)，所以上式就是基本的光流约束方程。
其中令表示图像中像素点灰度沿x，y，t方向的梯度，可将上式改写成：

Lucas-Kanade改进算法

金字塔的建立

令I0 = I 是第 0 层的图像，它是金字塔图像中分辨率最高的图像，图像的宽度和高度分别定义为nx0 = nx 和 ny0 = ny 。以一种递归的方式建立金字塔：从I0中计算I1，从I1中计算I2 ，···。令L =1, 2,...代表金字塔的层数，L通常取2,3,4。IL?1 是第L?1层的图像，nxL?1 和 nyL?1分别是图像IL?1 的宽度和高度。图像IL可按如下方式由IL?1 求得：

即用一个[0.25 0.5 0.25]的低通滤波器对IL-1进行卷积。

金字塔跟踪

总体来讲，金字塔特征跟踪算法描述如下：首先，光流和仿射变换矩阵在最高一层的图像上计算出；将上一层的计算结果作为初始值传递给下一层图像。在初始值的基础上，计算这一层的光流和仿射变化矩阵；再将这一层的光流和仿射矩阵作为初始值传递给下一层图像，直到传递给最后一层，即原始图像层，这一层计算出来的光流和仿射变换矩阵作为最后的光流和仿射变换矩阵的结果。
对于L=0,1,2,…L，定义是图像中像素点u在第L层对应点的坐标。根据上一步中图像金字塔的定义，可以计算出

我们用数学的思想重新描述在L层和L+1层迭代运算，假设对于最上层的光流计算初值为，并且对于最上层的变换矩阵猜测

为了在L层上计算光流和仿射变换矩阵，需要重新定义
在L层上的匹配误差ξL：

其中图像和是原始图像在L层上采样出来的图像，基于这层中的光流和仿射矩阵初值gL和GL可以计算出两个对应图像和：

通过观察匹配误差可知，对于匹配误差比较小的时候，改变窗口大小对于特征点追踪是可行的。接下来就是计算该层上的光流dL和变换矩阵AL，我们将在下一步中谈论。现在，假设在这一层上的光流和变换矩阵己经计算出来。接着将结果传递给下一层，计算出下一层的假设初值：

将gL-1和GL-1作为初值，重新循环上面的步骤，直到最上一层，计算出光流d和仿射变换矩阵A。设顶层时的初始为：

这种算法最明显的优势在于对于每一层的光流都会保持很小，但是最终计算来的光流可以进行累积，便于有效地跟踪特征点。

迭代过程

这一步是算法的核心步骤。在金字塔的每一层，目标是计算出光流dL和仿射变换矩阵AL从而使误差ξL最小。由于每一层的迭代过程是相同的，所以我们就描述从一层到下一层的迭代过程。首先将上一层的光流u和A传给这一层，计算这一帧图像中像素点的光照，同时计算出图像在该点x方向和y方向上的偏导

Ix=[I(x+1,y)-I(x-1,y)]/2

Iy=[I(x,y+1)-I(x,y-1)]/2

在此基础上，计算出空间梯度矩阵：

更新光流v=2*v

迭代过程：计算后一帧图像中对应像素点的灰度，计算两
帧图像间相同位置点的灰度值之差，在计算图像之间的误差
向量：

最后计算针对仿射光流

，
更新跟踪结果

直到某个阈值，结束在这一层的迭代过程。

特征点选择

因此，可按照以下的步骤选择特征点：
1、计算图像 I 中每一个像素的矩阵G和最小特征值λm。
2、寻找整副图像中最小特征值 λm 中的最大特征值λmax。
3、保留最小特征值 λm 大于给定阈值的像素点。阈值通常取5% λmax ~10% λmax 。
4、保留 λm 局部最大值的像素：像素特征值 λm 大于其3*3 邻域中其他像素的特征值 λm 。
5、剔除像素密集区域中的一些像素，确保图像中相邻像素的距离都大于给定的阈值（常取5~10 pixels）。
上述操作完成后，图像 I 中剩下的像素即为选择的特征点，并作为跟踪特征点。特征点选择算法的步骤5 确保了特征点间的最小距离。

没有必要取一个大的综合窗口选择特征点（或计算矩阵G）。大量实验证明，wx = wy =1的 3*3 大小的综合窗口能够取得满意的效果。

算法流程

%%%  Project Title: Lucas Kanade Motion Estimation Using Pyramids (Level 4) 
%%%  Note: The project specification says use density of 10 in plotting 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
%%% Usage: Lucas_Kanade('1.bmp','2.bmp',10) 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function Lucas_Kanade(file1,file2,density); 
%% Read Images %% 
img1 = im2double (imread (file1)); 
%% Take alternating rows and columns %% 
[odd1, even1] = split (img1); 
img2 = im2double (imread (file2)); 
[odd2, even2] = 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) 
    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); 
result = zeros (m1, n1); 
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; 
result = zeros (m1, n1); 
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'); 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, n] = 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 %% 
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 

光流法optical flow

标签：光流法

原文地址：http://blog.csdn.net/u014568921/article/details/46638557