MATLAB实现的金字塔光流算法

MATLAB实现的金字塔Lucas-Kanade（LK）光流算法的代码

function [u, v] = pyramidLK(img1, img2, levels, windowSize)
% 金字塔Lucas-Kanade光流算法
% 输入:
%   img1, img2 - 连续两帧灰度图像
%   levels - 金字塔层数
%   windowSize - 光流计算窗口大小（奇数）
% 输出:
%   u, v - 光流的水平和垂直分量

% 转换为双精度浮点以进行计算
img1 = im2double(img1);
img2 = im2double(img2);

% 创建图像金字塔
pyramid1 = cell(1, levels);
pyramid2 = cell(1, levels);
pyramid1{1} = img1;
pyramid2{1} = img2;

% 构建高斯金字塔
for i = 2:levels
    pyramid1{i} = impyramid(pyramid1{i-1}, 'reduce');
    pyramid2{i} = impyramid(pyramid2{i-1}, 'reduce');
end

% 初始化顶层光流
u = zeros(size(pyramid1{levels}));
v = zeros(size(pyramid1{levels}));

% 从金字塔顶层开始向下迭代
for level = levels:-1:1
    % 获取当前层图像
    I1 = pyramid1{level};
    I2 = pyramid2{level};
    
    % 如果当前层不是顶层，则上采样光流
    if level < levels
        % 上采样并缩放光流
        u = 2 * imresize(u, size(I1), 'bilinear');
        v = 2 * imresize(v, size(I1), 'bilinear');
    end
    
    % 根据当前光流对I2进行扭曲
    warpedI2 = warpImage(I2, u, v);
    
    % 计算当前层的光流增量
    [du, dv] = lkFlow(I1, warpedI2, windowSize);
    
    % 更新光流估计
    u = u + du;
    v = v + dv;
end

% 绘制光流场
plotFlow(u, v, img1);
end

function warped = warpImage(img, u, v)
% 使用光流(u,v)对图像进行扭曲
[xx, yy] = meshgrid(1:size(img,2), 1:size(img,1));
warped = interp2(xx, yy, img, xx+u, yy+v, 'linear', 0);
end

function [u, v] = lkFlow(img1, img2, windowSize)
% 单层Lucas-Kanade光流计算
% 计算图像梯度
[Ix, Iy] = gradient(img1);
It = img2 - img1;  % 时间导数

% 初始化光流
u = zeros(size(img1));
v = zeros(size(img1));

% 窗口半径
halfWin = floor(windowSize/2);

% 遍历图像中的每个像素
for i = halfWin+1:size(img1,1)-halfWin
    for j = halfWin+1:size(img1,2)-halfWin
        % 提取当前窗口
        winIx = Ix(i-halfWin:i+halfWin, j-halfWin:j+halfWin);
        winIy = Iy(i-halfWin:i+halfWin, j-halfWin:j+halfWin);
        winIt = It(i-halfWin:i+halfWin, j-halfWin:j+halfWin);
        
        % 构造矩阵A和向量b
        A = [winIx(:), winIy(:)];
        b = -winIt(:);
        
        % 求解最小二乘问题 A^T A [u;v] = A^T b
        if rank(A'*A) >= 2
            solution = (A'*A) \ (A'*b);
            u(i,j) = solution(1);
            v(i,j) = solution(2);
        end
    end
end

% 中值滤波去除异常值
u = medfilt2(u, [5,5]);
v = medfilt2(v, [5,5]);
end

function plotFlow(u, v, img)
% 绘制光流场
figure;
imshow(img);
hold on;

% 每隔10个像素绘制一个箭头
step = 10;
[x, y] = meshgrid(1:step:size(img,2), 1:step:size(img,1));
u_show = u(1:step:end, 1:step:end);
v_show = v(1:step:end, 1:step:end);

% 归一化箭头长度用于可视化
max_flow = max(sqrt(u_show(:).^2 + v_show(:).^2));
if max_flow > 0
    scale = 3 / max_flow;
else
    scale = 1;
end

quiver(x(:), y(:), u_show(:)*scale, v_show(:)*scale, 0, 'y', 'LineWidth', 1.5);
title('光流场可视化');
hold off;
end

使用

% 读取连续两帧图像
img1 = imread('frame1.jpg');
img2 = imread('frame2.jpg');

% 转换为灰度图
if size(img1, 3) == 3
    img1 = rgb2gray(img1);
    img2 = rgb2gray(img2);
end

% 设置参数并计算光流
levels = 3;        % 金字塔层数
windowSize = 15;   % 窗口大小（奇数）
[u, v] = pyramidLK(img1, img2, levels, windowSize);

说明：

1. 金字塔结构：

   • 通过impyramid函数构建高斯金字塔
   • 从顶层（低分辨率）开始计算，逐层向下细化
   • 每层光流上采样后作为下一层的初始估计

2. 核心LK光流计算：

   • 空间梯度计算：使用gradient函数计算x和y方向导数
   • 时间导数：相邻帧之间的像素差
   • 最小二乘求解：在局部窗口内求解光流方程

   \(\begin{bmatrix} \sum I_x^2 & \sum I_x I_y \\ \sum I_x I_y & \sum I_y^2 \end{bmatrix} \begin{bmatrix} u \\ v \end{bmatrix} = \begin{bmatrix} -\sum I_x I_t \\ -\sum I_y I_t \end{bmatrix}\)

3. 图像扭曲：

   • 使用双线性插值(interp2)根据当前光流估计扭曲第二帧图像
   • 迭代优化光流估计

4. 后处理：

   • 中值滤波(medfilt2)去除异常值
   • 光流场可视化：通过quiver函数绘制箭头图

参考代码 金字塔LK光流计算的代码 www.youwenfan.com/contentcnd/98164.html

参数选择建议：

1. 金字塔层数(levels)：

   • 通常3-5层
   • 对于大位移场景需要更多层

2. 窗口大小(windowSize)：

   • 典型值：15×15或21×21
   • 太小：对噪声敏感
   • 太大：模糊运动边界

3. 迭代次数：

   • 本代码每层单次迭代（实际应用中可增加迭代次数）
   • 可在lkFlow函数中添加迭代循环

改进方向：

1. 增加迭代：在每层金字塔内添加多次迭代
2. 亚像素精度：使用更精细的插值方法
3. 鲁棒估计：引入M估计或双向一致性检查
4. 平滑约束：加入全局平滑项（类似Horn-Schunck方法）

此实现适用于中等运动场景，对于快速运动或低纹理区域，可能需要结合其他技术（如特征点匹配）来提高鲁棒性。