当前位置：首页 > news >正文

MATLAB实现的金字塔光流算法

news 2025/11/7 7:56:27

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:levelspyramid1{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);
endfunction 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);
endfunction [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)-halfWinfor 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和向量bA = [winIx(:), winIy(:)];b = -winIt(:);% 求解最小二乘问题 A^T A [u;v] = A^T bif rank(A'*A) >= 2solution = (A'*A) \ (A'*b);u(i,j) = solution(1);v(i,j) = solution(2);endend
end% 中值滤波去除异常值
u = medfilt2(u, [5,5]);
v = medfilt2(v, [5,5]);
endfunction 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 > 0scale = 3 / max_flow;
elsescale = 1;
endquiver(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) == 3img1 = rgb2gray(img1);img2 = rgb2gray(img2);
end% 设置参数并计算光流
levels = 3;        % 金字塔层数
windowSize = 15;   % 窗口大小（奇数）
[u, v] = pyramidLK(img1, img2, levels, windowSize);

说明：

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

核心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}