L-K金字塔光流法matlab代码实现

L-K金字塔光流法matlab代码实现，参考1981年An iterative image registration technique with an application to stereo vision文章

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Name: Musawir Ali Shah
%%%  Assignment #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（file1file2density）;
%% Read Images %%
im11=imread （file1）;
if size（size（im11）2）==3
    im1=rgb2gray（im11）;
end
img1 = im2double （im1）;
%% Take alternating rows and columns %%
 %[odd1 even1] = split （img1）;
im22=imread （file2）;
if size（size（im22）2）==3
    im2=rgb2gray（im22）;
end
img2 = im2double（im2）;
 %[odd2 even2] = split （img2）;
%% Run Lucas Kanade %%
 %[Dx Dy] = Estimate （odd1 odd2）;
[Dx Dy] = Estimate （img1 img2）;
%% Plot %%
figure; 
imshow（imread（file1））;
hold on;
[maxImaxJ]=size（Dx）;
Dx=Dx（1:density:maxI1:density:maxJ）;
Dy=Dy（1:density:maxI1:density:maxJ）;
quiver（1:density:maxJ（maxI）:（-density）:1Dx-Dy1）;
hold on;
[maxImaxJ]=size（Dx1）;
Dx1=Dx1（1:density:maxI1:density:maxJ）;
Dy1=Dy1（1:density:maxI1:density:maxJ）;
quiver（1:density:maxJ（maxI）:（-density）:1Dx1-Dy11）;
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:m1:n） + Vx; Dy（1:m 1:n） = Dy（1:m 1:n） + Vy;
    smooth （Dx）; smooth （Dy）;
    l = l - 1;
end


%% Lucas Kanad

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     目录           0  2019-02-25 22:14  Lucas_Kanade\
     文件       64078  2003-04-08 18:03  Lucas_Kanade\1.bmp
     文件       81110  2009-04-09 13:04  Lucas_Kanade\11.bmp
     文件       64078  2003-04-08 18:03  Lucas_Kanade\2.bmp
     文件       82302  2009-04-09 13:04  Lucas_Kanade\22.bmp
     文件       21958  2003-04-07 14:03  Lucas_Kanade\3.bmp
     文件       21958  2003-04-07 14:03  Lucas_Kanade\4.bmp
     文件       82302  2009-04-09 14:07  Lucas_Kanade\a.bmp
     文件       82302  2009-04-09 14:08  Lucas_Kanade\b.bmp
     文件      132186  2019-02-25 20:46  Lucas_Kanade\frame10.jpg
     文件      403108  2007-07-28 03:16  Lucas_Kanade\frame10.png
     文件      128038  2019-02-25 20:46  Lucas_Kanade\frame11.jpg
     文件      402536  2007-07-28 03:16  Lucas_Kanade\frame11.png
     文件      169984  2004-04-18 03:57  Lucas_Kanade\Lucas.doc
     文件        7312  2019-02-25 22:35  Lucas_Kanade\Lucas_Kanade.m
     文件         310  2019-02-25 22:06  Lucas_Kanade\tests.m
     文件          31  2012-03-15 16:52  Lucas_Kanade\Untitled.m