关于粒子平流与FTLE场计算的代码说明注释

首先，粒子平流需要事先计算好每两帧图像的光流信息，然后利用滑动窗技术计算某段时间内的平均光流信息，利用平均光流信息进行粒子平流，最后，可以利用粒子平流的结果计算FTLE场，后者据说能更好地描述光流的运动信息。关于上述技术，具体可参考“S . Ali and M . Shah. A lagrangian particle dynamics approach for crowd flow segmentation and  stability analysis .“与"Pa r a l l e l Pa r t i c l e Adve c t i on a nd FTLE Com put a t i on f or Ti m e - Va r yi ng Fl ow Fi e l ds".

下面主要标注下粒子平流与FTLE场计算的代码：（PS：该段代码仅仅计算了某段时间内的粒子平流与FTLE场，粒子平流可以得到多个结果，譬如时间段为30，就能得到大致29个结果，FTLE场仅能得到一个结果，譬如FTLE1-30）

 

PS:FTLE场的作用：把群体运动在某段视频内看作为静止的，这样就能够从视频中分割出群体，适用于那些群体密度等级非常大，而且人群流经的路径也比较稳定的场合。

%% 下面这段代码实现的功能是从预存mat文件中读取每两帧图像的光流信息；

optical_flow_file_names = dir([optical_flow_folder, '\*.mat']);

matMotionFileName = fullfile ( optical_flow_folder, optical_flow_file_names(start_frame).name );

load(matMotionFileName);

%% um和vm中保存的就是u和v，分别是x方向和y方向的光流信息；

um = nan2zeros(u);
vm = nan2zeros(v);

counter = 1;

%% 下面这段代码实现的功能是计算某段时间内的平均光流

for i = start_frame + 1  : end_frame

    matMotionFileName = fullfile ( optical_flow_folder, optical_flow_file_names(i).name );

    load(matMotionFileName);

    um = um + nan2zeros(u);
    vm = vm + nan2zeros(v);

    counter = counter + 1;

end

%% 对叠加的粒子流进行了平均
um = um ./ counter;
vm = vm ./ counter;

%% 下面这段代码实现的功能是设定粒子平流所需要的参数

%% x、y是带标号的网格，尺寸与输入图像相同，x是列相同的1、2、3...，y是行相同的1、2、3...。
x1 = min(x(:));
x2 = max(x(:));
y1 = min(y(:));
y2 = max(y(:));

xmesh = x;

ymesh = y;

dt       = ftle_options.step_size;  %%每秒8帧，因此每帧1/8秒；

frame_rate = ftle_options.frame_rate;

t_length = end_frame - start_frame + 1;

T_span  = (end_frame/frame_rate) - (start_frame/frame_rate) + (1/frame_rate);

fprintf('Forward particle advection --- Frame range: %d....%d \n', start_frame, end_frame);

 

%%

%%% Assuming that the mean field is representative of the motion in this
%%% block of frames. Since algorithm is working on a sliding window any
%%% changes in the dynamic behavior of the underlying flow will be captured
%%% by later blocks. This step also helps in reducing the effect of noise
%%% from the optical flow algorithm.
%% 输入的um和vm是平均光流场；
u = um;
v = vm;

for t_integration = 1 : t_length - 1;

    index = start_frame + t_integration - 1;

    if ftle_options.directional_segmentation == true

        [u, v] = normalize_magnitude(u,v);  %%进行标准化处理

        u = nan2zeros(u);

        v = nan2zeros(v);
    end

%% 以下代表进行粒子平流

%%%%%%Advection（进行粒子平流）%%%%%%%%%%
    if t_integration == 1   %% 初始情况的处理；

        xflowmap{t_integration} = xmesh;

        yflowmap{t_integration} = ymesh;

    else     %%利用线性插值完成粒子平流的计算
%   Vq = INTERP2(X,Y,V,Xq,Yq) interpolates to find Vq, the values of the
%   underlying 2-D function V at the query points in matrices Xq and Yq.
%   Matrices X and Y specify the points at which the data V is given.
%   对于上面这段注释，有V=F（X,Y），即V与XY有某种函数关系，通过XY与XqYq的近似关系获取V的插值Vq；
%   V和XY有一个函数关系，根据这个函数关系，估计Xq,Yq对应的值Vq；
%   dt表示每帧1/8秒；
%   根据平均光流信息u来计算插值信息；
       
        xflowmap{t_integration} = xflowmap{t_integration-1} + dt*interp2(xmesh, ymesh, u, xflowmap{t_integration-1}, yflowmap{t_integration-1}, 'linear', 0);

        yflowmap{t_integration} = yflowmap{t_integration-1} + dt*interp2(xmesh, ymesh, v, xflowmap{t_integration-1}, yflowmap{t_integration-1}, 'linear', 0);

    end

end

%% 根据粒子平流的结果计算FTLE场，貌似只用到了粒子平流的最后一组结果

fprintf('Computing forward FTLE \n');
%% 计算FTLE场，也就是计算李雅普诺夫熵；
[xFX,xFY] = gradient(xflowmap{end}, xmesh(1,2) - xmesh(1,1));
[yFX,yFY] = gradient(yflowmap{end},  ymesh(2,1)- ymesh(1,1));

sigma = zeros(size(xFX));

for i =1 : size(xFX,1)

    for j = 1 : size(xFX,2)

        A11= xFX(i, j);
        A12= xFY(i, j);
        A21= yFX(i, j);
        A22= yFY(i, j);

        A=[A11 A12;A21 A22];

        B=A'*A;

        delta=max(eig(B));

        sigma(i,j) = log(delta)/(2*T_span);

    end

end

%% 保存FTLE场的结果

%%%%%%%Save FTLE data
[pfx_crowd_folder,garbage] = fileparts(optical_flow_folder);

ftle_folder     = fullfile(pfx_crowd_folder, 'FTLE');

if ~exist([ftle_folder])
    mkdir(ftle_folder);
end

matFTLEFileName = fullfile ( ftle_folder, sprintf('ForwardFTLE%04d-%04d.mat', start_frame, end_frame ));
save(matFTLEFileName, 'ftle_options', 'sigma');