由於工程需要用到 Lucas-Kanade 光流，在此進行一下簡單整理（后續還會陸續整理關於KCF，PCA，SVM，最小二乘、嶺回歸、核函數、dpm等等）：
光流，簡單說也就是畫面移動過程中，圖像上每個像素的x，y位移量，比如第t幀的時候A點的位置是(x1, y1)，那么我們在第t+1幀的時候再找到A點，假如它的位置是(x2,y2)，那么我們就可以確定A點的運動了：(u, v) = (x2, y2) - (x1,y1)
1、假設原圖是I(x,y,z) (這里是擴展到三維空間的，所以還有個z值)，移動后的圖像是I(x+δx,y+δy,z+δz,t+δt),兩者滿足:
2、其中圖像移動可以認為I (x ,y ,z ,t ) = I (x + δx ,y + δy ,z + δz ,t + δt )
也就是說：
（ H.O.T. 指更高階，在移動足夠小的情況下可以忽略)
3、從這個方程中我們可以得到：
其中Vx = u， Vy=v，也就是光流的值（二維圖像沒有z），
則是圖像在(x ,y,z ,t )這一點的梯度 （
就是兩幀圖像塊之間差值） 。
4、假設流(Vx,Vy,Vz)在一個大小為m*m*m(m>1)的小窗中是一個常數，那么從像素1...n , n = m*m*m 中可以得到下列一組方程：
三個未知數但是有多於三個的方程，這個方程組自然是個超定方程，也就是說方程組內有冗余，方程組可以表示為：
也就是：
采用最小二乘法：
5、另外，由於LK算法假設是小位移，為了解決大位移問題，需要在多層圖像縮放金字塔上求解，每一層的求解結果乘以2后加到下一層：
6、具體就見matlab代碼：
其中求解最小二乘的行列式求解只有2維所以計算量尚可容忍
%Data acquisition
im1= ((imread('1.png')));
im2= ((imread('2.png')));
im1=single(im1);
im2=single(im2);
[result,corner_count,ptx,pty] = harris(im1); //harris角點是求光流的關鍵點
imagesc(result);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%parameters : levels number, window size, iterations number, regularization
numLevels= 4;
window= 10;
iterations=3;
alpha = 0.001;
hw = floor(window/2);
t0 = clock;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%pyramids creation
pyramid1 = im1;
pyramid2 = im2;
%init
for i=2:numLevels
im1 = impyramid(im1, 'reduce');
im2 = impyramid(im2, 'reduce');
pyramid1(1:size(im1,1), 1:size(im1,2), i) = im1;
pyramid2(1:size(im2,1), 1:size(im2,2), i) = im2;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Processing all levels
for p = 1:numLevels
%current pyramid
im1 = pyramid1(1:(size(pyramid1,1)/(2^(numLevels - p))), 1:(size(pyramid1,2)/(2^(numLevels - p))), (numLevels - p)+1);
im2 = pyramid2(1:(size(pyramid2,1)/(2^(numLevels - p))), 1:(size(pyramid2,2)/(2^(numLevels - p))), (numLevels - p)+1);
%init
if p==1
u = zeros(size(im1));
v = zeros(size(im1));
else
%resizing
u = 2 * imresize(u,size(u)*2,'bilinear');
v = 2 * imresize(v,size(v)*2,'bilinear');
end
%refinment loop
for r = 1:iterations
u=round(u);
v=round(v);
%every pixel loop
for i = 1+hw:size(im1,1)-hw
for j = 1+hw:size(im2,2)-hw
patch1 = im1(i-hw:i+hw, j-hw:j+hw);
%moved patch
lr = i-hw+v(i,j);
hr = i+hw+v(i,j);
lc = j-hw+u(i,j);
hc = j+hw+u(i,j);
if (lr < 1)||(hr > size(im1,1))||(lc < 1)||(hc > size(im1,2))
%Regularized least square processing
else
patch2 = im2(lr:hr, lc:hc);
fx = conv2(patch1, 0.25* [-1 1; -1 1]) + conv2(patch2, 0.25*[-1 1; -1 1]);
fy = conv2(patch1, 0.25* [-1 -1; 1 1]) + conv2(patch2, 0.25*[-1 -1; 1 1]);
ft = conv2(patch1, 0.25*ones(2)) + conv2(patch2, -0.25*ones(2));
Fx = fx(2:window-1,2:window-1)';
Fy = fy(2:window-1,2:window-1)';
Ft = ft(2:window-1,2:window-1)';
A = [Fx(:) Fy(:)];
G=A'*A;
G(1,1)=G(1,1)+alpha; G(2,2)=G(2,2)+alpha;
U=1/(G(1,1)*G(2,2)-G(1,2)*G(2,1))*[G(2,2) -G(1,2);-G(2,1) G(1,1)]*A'*-Ft(:);
u(i,j)=u(i,j)+U(1); v(i,j)=v(i,j)+U(2);
end
end
end
end
etime(clock,t0)
end