光流法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)，所以上式就是基本的光流约束方程。
其中令
Lucas-Kanade改进算法
在实际场景中图像I和图像J可以代表前后两帧图像。对于图像特征点金字塔跟踪来说的目的是：对于前一帧的图像I上一点u(ux,uy)，要在后一帧图像J上找到一点v(ux+dx,uy+dy)与之相匹配，即灰度值最接近。那么向量d=[dx,dy]就是图像在点u处的运动速度，也就是所说像素点u的光流。为了进一步说明向量d的含义。我们假设前一帧图像经历了仿射变换到后一帧图像，定义变换矩阵为
定义误差
典型的wx和wy取值为7，8，10，20个像素。
金字塔的建立
令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，定义
我们用数学的思想重新描述在L层和L+1层迭代运算，假设对于最上层的光流计算初值为
为了在L层上计算光流和仿射变换矩阵，需要重新定义
在L层上的匹配误差ξL：
其中图像
通过观察匹配误差可知，对于匹配误差比较小的时候，改变窗口大小对于特征点追踪是可行的。接下来就是计算该层上的光流dL和变换矩阵AL，我们将在下一步中谈论。现在，假设在这一层上的光流和变换矩阵己经计算出来。接着将结果传递给下一层，计算出下一层的假设初值：
将gL-1和GL-1作为初值，重新循环上面的步骤，直到最上一层，计算出光流d和仿射变换矩阵A。设顶层时的初始为：
这种算法最明显的优势在于对于每一层的光流都会保持很小，但是最终计算来的光流可以进行累积，便于有效地跟踪特征点。
迭代过程
这一步是算法的核心步骤。在金字塔的每一层，目标是计算出光流dL和仿射变换矩阵AL从而使误差ξL最小。由于每一层的迭代过程是相同的，所以我们就描述从一层到下一层的迭代过程。首先将上一层的光流u和A传给这一层，计算这一帧图像中像素点的光照
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