光流场计算与HS算法-CSDN博客

本文链接：https://blog.csdn.net/jpc20144055069/article/details/96475084

前言

光流是指空间运动物体在观测成像面上的像素运动的瞬时速度，它利用图像序列像素强度的时域变化和相关性来确定各自像素位置的“运动”，即反映图像灰度在时间上的变化与景物中物体结构与其运动的关系。将二维图像平面特定坐标上的灰度瞬时变化率定义为光流矢量。

光流的物理意义

当人的眼睛观察运动物体时，物体的景象在人眼的视网膜上形成一系列连续变化的图像，这一系列连续变化的信息不断“流过”视网膜（即图像平面），好像一种光的“流”，故称之为光流。光流表达了图像的变化，由于它包含了目标运动的信息，因此可被观察者用来确定目标的运动情况。

光流场

在空间中，运动可以用运动场描述，而在一个图像平面上，物体的运动往往是通过图像序列中不同图像灰度分布的不同体现的，从而，空间中的运动场转移到图像上就表示为光流场（optical flow field）。
　　光流场是一个二维矢量场，它反映了图像上每一点灰度的变化趋势，可看成是带有灰度的像素点在图像平面上运动而产生的瞬时速度场。它包含的信息即是各像点的瞬时运动速度矢量信息。
　　研究光流场的目的就是为了从序列图像中近似计算不能直接得到的运动场。光流场在理想情况下，光流场对应于运动场。

三维空间的矢量场及其在二维平面内的投影

光流场计算的基本原理

设在t 时刻，像素点 (x,y) 处的灰度值为I(x,y,t)；在(t+ $\Delta$ t)，该点运动到新的位置，它在图像上的位置变为(x+ $\Delta$ x, y+ $\Delta$ y)，灰度值记为I(x+ $\Delta$ x, y+ $\Delta$ y)。根据图像一致性假设，即图像沿着运动轨迹的亮度保持不变，满足 $\frac{dI(x,y,t)}{dt}$ =0，则： $I(x,y,t)=I(x+\Delta x, y+\Delta y + t+\Delta t)\tag {1}$ 　　设u和 $\upsilon$ 分别为该点的光流矢量沿x和y方向的两个分量，且 u= $\frac{dx}{dt}$ ， $\upsilon$ = $\frac{dy}{dt}$ ，将式(1)等号右边用泰勒公式展开看，得到： $I(x+\Delta x, y+\Delta y, t+ \Delta t) = I(x,y,t)+\frac{\partial I}{\partial x}\Delta x+\frac{\partial I}{\partial y}\Delta y+\frac{\partial I}{\partial t}\Delta t+\varepsilon \tag {2}$ 　　忽略二阶以上的高次项，则有： $\frac{\partial I}{\partial x}\Delta x+\frac{\partial I}{\partial y}\Delta y+\frac{\partial I}{\partial t}\Delta t=0 \tag{3}$ 　　由于 $\Delta t\to 0$ ，于是有： $\frac{\partial I}{\partial x}\frac{dx}{dt}+\frac{\partial I}{\partial y}\frac{dy}{dt}+\frac{\partial I}{\partial t} = 0\tag{4}$ 　　也即： $I_xu+I_y\upsilon+I_t= 0\tag{5}$ 　　式（5）是光流基本等式。设 $I_x、I_y和I_t$ 分别为参考点像素的灰度值沿x,y,t这三个方向的偏导数。
　　Horn-Schunck算法引用的附加约束条件的基本思想是：在求解光流时，要求光流本身尽可能地平滑，即引入对光流的整体平滑性约束求解光流方程病态问题。所谓平滑，就是在给定的领域内 $(\nabla^2u+\nabla^2\upsilon)$ 应尽量地小，这就是求条件极值时的约束条件。对 $u、\upsilon$ 的附加条件如下： $min\{[\frac{\partial u}{\partial x}]^2+[\frac{\partial u}{\partial y}]^2]+[\frac{\partial \upsilon}{\partial x}]^2+[\frac{\partial \upsilon}{\partial y}]^2\}\tag{6}$ 　　式中， $\nabla ^2u=[\frac{\partial u}{\partial x}]^2+[\frac{\partial u}{\partial y}]^2$ 是 $u$ 的拉普拉斯算子， $\nabla ^2\upsilon=[\frac{\partial u}{\partial x}]^2+[\frac{\partial u}{\partial y}]^2$ 是 $\upsilon$ 的拉普拉斯算子。综合式(5)和式(6)，Horn-Schunck算法将光流 $u、\upsilon$ 的计算归结为如下问题： $min\{\iint (I_xu+I_y\upsilon + I_t)^2 + \alpha ^ 2 [[\frac{\partial u}{\partial x}]^2+[\frac{\partial u}{\partial y}]^2+[\frac{\partial \upsilon}{\partial x}]^2+[\frac{\partial \upsilon}{\partial y}]^2]\}\tag{7}$ 　　因而，可以得到其相应的欧拉-拉格朗日方程，并利用高斯-塞德尔方法进行求解，得到图像每个次置第一次至第(n+1)次迭代估计 $(u^{n+1}, \upsilon^{n+1})$ 为：

$u^{n+1} = \overline{u^n}-\overline{I_x^2}\frac{\overline{I_x}\overline{u^n}+\overline{I_y}\overline{\upsilon^n}+I_t}{\alpha^2 +\overline{I_x}^2+ \overline{I_y}^2}\tag{8}$ $\upsilon^{n+1} = \overline{u^n}-\overline{I_y^2}\frac{\overline{I_x}\overline{u^n}+\overline{I_y}\overline{\upsilon^n}+I_t}{\alpha^2 +\overline{I_x}^2+ \overline{I_y}^2}\tag{9}$

求解过程要得到稳定的解通常需要上百次的迭代。整个迭代过程既与图像尺寸有关，又与每次的传递量（速度的改变量）有关。由迭代公式可以发现，在一些缺乏特征且较为平坦的区域（梯度为0或较小），其速度由迭代公式的第一项决定，该点的速度信息需要从特征较为丰富的区域传递过来。为了加快算法的收敛速度，一方面可以用金字塔的层次结构来减小图像的尺寸扩散，另一方面可以采用增加扩散量的方法来加速算法。

Horn-Schunck算法MATLAB实现程序如下：

1.求导函数 computeDerivatices

function [fx, fy, ft] = computeDerivatives(im1, im2)
% 功能：求输入图像参考像素点的像素值沿三轴方向的偏导数
% 输入：
%     im1-输入图像1
%     im2-输入图像2
% 输出：
%     fx-参考像素点的灰度值沿x方向的偏导数
%     fy-参考像素点的灰度值沿y方向的偏导数
%     fz-参考像素点的灰度值沿z方向的偏导数
if size(im2,1)==0
    im2=zeros(size(im1));
end
% 利用标准模板求得式（5）中的偏导数Ix, Iy, It
fx = conv2(im1,0.25* [-1 1; -1 1],'same') + conv2(im2, 0.25*[-1 1; -1 1],'same');
fy = conv2(im1, 0.25*[-1 -1; 1 1], 'same') + conv2(im2, 0.25*[-1 -1; 1 1], 'same');
ft = conv2(im1, 0.25*ones(2),'same') + conv2(im2, -0.25*ones(2),'same');

2.高斯滤波函数 gaussFilter

function G=gaussFilter(segma,kSize)
% 功能：实现高斯滤波
% 输入：
%     sigma-高斯分布的概率密度函数的方差
%     kSize-高斯向量的模板尺寸大小
% 输出：
%     G-方差为segma，大小为kSize的一维高斯向量模板
if nargin<1
    segma=1;
end
if nargin<2
    kSize=2*(segma*3);
end
% 利用均值为0，方差为segma高斯分布概率密度函数求解一维高斯向量模板
G=(1/(sqrt(2*pi)*segma)) * exp (-(x.^2)/(2*segma^2));

3.平滑性约束条件函数 smoothImg

function smoothedImg=smoothImg(img,segma)
% 功能：实现平滑性约束条件
% 输入：
%     img-数字图像
%     sigma-高斯分布的方差
% 输出：
%     smoothedImg-经高斯滤波的图像矩阵

if nargin<2
    segma=1;
end
%调用高斯滤波函数
G=gaussFilter(segma);  
%根据（6）对图像进行平滑约束
smoothedImg=conv2(img,G,'same');  
% 二次平滑
smoothedImg=conv2(smoothedImg,G','same');

4.求解光流场函数 HS

function [u, v] = HS(im1, im2, alpha, ite, uInitial, vInitial, displayFlow, displayImg)
% 功能：求解光流场
% 输入：
%     im1-输入图像1
%     im2-输入图像2
%     alpha-反映HS光流算法的平滑性约束条件的参数
%     ita-（8）和（9）式中的迭代次数
%     uInitial-光流横向分量初始值
%     vInitial-光流纵向分量初始值
%     displayFlow-光流场显示参数，其值为1时显示，为0时不显示
%     displayImg-显示光流场的指定图像，如果为空矩阵，则无指定图像输出
% 输出：
%     u-横向光流矢量
%     v-纵向光流矢量

%  初始化参数
if nargin<1 || nargin<2
    im1=imread('HS1.tif');
    im2=imread('HS2.tif');
end
if nargin<3
    alpha=1;
end
if nargin<4
    ite=100;
end
if nargin<5 || nargin<6
    uInitial = zeros(size(im1(:,:,1)));
    vInitial = zeros(size(im2(:,:,1)));
elseif size(uInitial,1) ==0 || size(vInitial,1)==0
    uInitial = zeros(size(im1(:,:,1)));
    vInitial = zeros(size(im2(:,:,1)));
end
if nargin<7
    displayFlow=1;
end
if nargin<8
    displayImg=im1;
end
%  将RGB图像转化为灰度图像
if size(size(im1),2)==3
    im1=rgb2gray(im1);
end
if size(size(im2),2)==3
    im2=rgb2gray(im2);
end
im1=double(im1);
im2=double(im2);
% 调用平滑性约束函数对图像进行平滑
im1=smoothImg(im1,1);
im2=smoothImg(im2,1);
tic;
% 为光流矢量设置初始值
u = uInitial;
v = vInitial;
[fx, fy, ft] = computeDerivatives(im1, im2); % 调用求导函数对时间分量和空间分量进行求导
kernel_1=[1/12 1/6 1/12;1/6 0 1/6;1/12 1/6 1/12]; % 均值模板
% 根据式（3.15.9）迭代求解，迭代次数为100
for i=1:ite
    % 计算光流矢量的局部均值
    uAvg=conv2(u,kernel_1,'same');
    vAvg=conv2(v,kernel_1,'same');
    % 根据式（3.15.9）用迭代法求解光流矢量
    u= uAvg - ( fx .* ( ( fx .* uAvg ) + ( fy .* vAvg ) + ft ) ) ./ ( alpha^2 + fx.^2 + fy.^2); 
    v= vAvg - ( fy .* ( ( fx .* uAvg ) + ( fy .* vAvg ) + ft ) ) ./ ( alpha^2 + fx.^2 + fy.^2);
end
u(isnan(u))=0;
v(isnan(v))=0;
% 画图
if displayFlow==1
   plotFlow(u, v, displayImg, 5, 5);  % 调用画图函数
end

5.画图函数 plotFlow

function plotFlow(u, v, imgOriginal, rSize, scale)
% 功能：绘制光流场图
% 输入：u-横向光流矢量     v-纵向光流矢量     imgOriginal-光流场显示的图像
%     rSize-可见光流矢量区域的尺寸             scale-光流场规模
figure();
if nargin>2
    if sum(sum(imgOriginal))~=0
        imshow(imgOriginal,[0 255]);
        hold on;
    end
end
if nargin<4
    rSize=5;
end
if nargin<5
    scale=3;
end
for i=1:size(u,1)
    for j=1:size(u,2)
        if floor(i/rSize)~=i/rSize || floor(j/rSize)~=j/rSize
            u(i,j)=0;
            v(i,j)=0;
        end
    end
end
quiver(u, v, scale, 'color', 'b', 'linewidth', 2);
set(gca,'YDir','reverse');