python光流法测速 光流法matlab

光流法指的是一种简单实用的图像运动的表达方式，通常定义为一个图像序列中的图像亮度模式的表观运动，即空间物体表面上的点的运动。光流的研究是利用图像序列中的像素强度数据的时域变化和相关性来确定各自像素位置的"运动"，即研究图像灰度在时间上的变化与景象中物体结构及其运动的关系。速度在视觉传感器的成像平面上的表达。

真正提出有效光流计算方法还归功于Horn和Schunck在1981年创造性地将二维速度场与灰度相联系，引入光流约束方程的算法，是光流算法发展的基石。

本文使用的就是这种算法。

原视频

识别结果

图中红色小箭头即是通过HS算法得出的像素点矢量。

附代码：

main.m

注：GPU在图像处理（单精度运算）的性能远超CPU，因此使用gpuarray（）函数将矩阵置于GPU运算大幅提高程序运行速度，此功能仅N卡支持，若设备无N卡将该函数注释即可。

% 基于Horn-Schunck算法的计算在视频中显示光流
% alpha是权重因子
% N是迭代次数
close all;

clc;%清除命令窗口

videopath='C:\Users\Administrator\Desktop\work2\video\';%视频帧路径
alpha = 0.5;
N = 50;%值越大精度越高，计算速度越慢
% 将视频读取到对象
vidObj = VideoReader('sample.avi');

% 获取第一帧
if hasFrame(vidObj) % matlab自带函数，确定是否有视频帧读取
    Fr1 = rgb2gray(readFrame(vidObj));%读取一个视频帧并转换为灰度图
    Fr1 = im2single(Fr1);%将图像转换为单精度值
    Fr1 = gpuArray(Fr1);%将矩阵数据置于显卡中运算
else
    return
end
k=0;
while hasFrame(vidObj)%当有视频帧时
    Fr2 = rgb2gray(readFrame(vidObj)); % 读取下一帧并转换为灰度图
    Fr2 = im2single(Fr2);%将图像转换为单精度值
    Fr2 = gpuArray(Fr2);%将矩阵数据置于显卡中运算
    k=k+1;
    % 通过HS算法计算U、V
    [U, V] = HS(Fr1, Fr2, alpha, N);
        
    % 显示框架和光流
   figure(1);
   
   imshow(Fr1,[]); %显示一帧
   hold on;
   axis image; % 将子框架设置为相同大小
   showOF(U, V); %函数showOF显示流速矩阵U、V的光流

   k1=num2str(k);
   filename=[k1,'.jpg'];
   saveas(1,[videopath,filename]);

   hold off;
   Fr1 = Fr2; % Fr2成为下一序列的第一帧

end

cd(videopath); %读取所有的jpg图片
allnames = struct2cell(dir('*.jpg'));
[difgrayFrame,len]=size(allnames);
aviobj = VideoWriter('C:\Users\Administrator\Desktop\work2\save_video2.avi');%视频存储位置
aviobj.FrameRate = 25; %设置帧率

open(aviobj) %制作视频
for i = 1:k
     name = [num2str(i),'.jpg'];
     frame = imread(name);
     writeVideo(aviobj,frame);
end
close(aviobj)

 derivative.m

function [Ex, Ey, Et] = derivative(Im1, Im2)
% 函数DERIVATIVE计算偏导数Ex、Ey、Et
% 双类别的两幅图像Im1、Im2的序列。

% 卷积核
Kx = 0.25 * [-1 1; -1 1];
Ky = 0.25 * [-1 -1; 1 1];
Kt = 0.25 * [-1 -1; -1 -1]; % kt1 = Kt, kt2 = -Kt

% 计算偏导数
Ex = conv2(Im1, Kx, 'same') + conv2(Im2, Kx, 'same');
Ey = conv2(Im1, Ky, 'same') + conv2(Im2, Ky, 'same');
Et = conv2(Im1, Kt, 'same') + conv2(Im2, -Kt, 'same');

HS.m(算法核心)

function [U, V] = HS(Im1, Im2, alpha, N)

%函数HS算法计算两幅图像序列的流速U、V
%基于Horn-Schunck算法的双类Im1、Im2。
%alpha是加权因子和N是迭代次数。

% 计算偏导数
[Ex, Ey, Et] = derivative(Im1, Im2);

% 初始 U, V
[l, c] = size(Im1);%获取数组大小l行c列
U = zeros(l, c);%零矩阵
V = zeros(l, c);

K = [1/12 1/6 1/12; 1/6 -1 1/6; 1/12 1/6 1/12]; % 拉普拉斯卷积核
A = alpha^2 + Ex.^2 + Ey.^2;

for i = 1:N
    % 计算U、V平均值
    U_avg = conv2(U, K, 'same');%返回U与K卷积中大小与U相同的中心部分
    V_avg = conv2(V, K, 'same');
    B = (Ex.*U_avg + Ey.*V_avg + Et);
    
    % 计算当前迭代的U、V
    U = U_avg - Ex.*B./A;
    V = V_avg - Ey.*B./A;
end

 showOF.m

function showOF(U, V)
% 函数showOF显示流速矩阵U、V的光流

[l,c] = size(U);
d = 5; % 步幅
IndX = 1:d:l;
IndY = 1:d:c;
[X,Y] = meshgrid(IndY, IndX);%二维网格
U1 = U(IndX, IndY);
V1 = V(IndX, IndY);
%在由 X 和 Y 指定的笛卡尔坐标上绘制具有定向分量 U 和 V 的箭头。
quiver(X, 720-Y, U1(end:-1:1,:), V1(end:-1:1,:), 3,Color='r');

识别结果：

matlab光流法实例