filterDisplacements_2D.m

function u = filterDisplacements_2D(u0,filterSize,z)
% I = filterDisplacements(I0,filterSize,z) applies a low-pass convolution
% filter to the displacement field to mitigate divergence based on
%
% F. F. J. Schrijer and F. Scarano. Effect of predictor corrector filtering
% on the stability and spatial resolution of iterative PIV interrogation.
% Exp. Fluids, 45(5):927{941, May 2008. doi: 10.1007/s00348-008-0511-7
%
% INPUTS
% -------------------------------------------------------------------------
%   u0: displacement field vector defined at every meshgrid point with
%      spacing dm. Format: cell array, each containing a 3D matrix
%         (components in x,y,z)
%         u0{1} = displacement in x-direction
%         u0{2} = displacement in y-direction
%         u0{3} = magnitude
%   filterSize: size of the filter
%   z: filter strength
%
% OUTPUTS
% -------------------------------------------------------------------------
%   u: cell containing the filtered displacement field
%
% NOTES
% -------------------------------------------------------------------------
% none
%
% For more information please see section 2.2.
% If used please cite:
% Landauer, A.K., Patel, M., Henann, D.L. et al. Exp Mech (2018).
% https://doi.org/10.1007/s11340-018-0377-4

% Parse inputs and set defaults
if nargin < 3,  z = 0.0075; end
if ~iscell(u0), u0 = {u0}; end
u = cell(size(u0));

if z == 0
    u = cellfun(@double, u0, 'UniformOutput',0); % no filter
else
    rf = generateFilter(filterSize,z);

    % apply filter using convolution
    for i = 1:length(u0), u{i} = double(convn(u0{i}, rf,'same')); end
end

end

%% ========================================================================
function rf = generateFilter(filterSize,z)
% generates the filter convolution filter
l = filterSize;

[m{1}, m{2}] = ndgrid(-l(1)/2:l(1)/2,-l(2)/2:l(2)/2);
m = cellfun(@abs, m, 'UniformOutput', 0);


f1 = (l(1)/2)^z - (m{1}).^z;
f2 = (l(2)/2)^z - (m{2}).^z;

i{1} = (m{1} >= m{2});% ?? & m{1} >= m{3});
i{2} = (m{1} < m{2});% ?? & m{2} >= m{3});

rf0 = f1.*i{1}+f2.*i{2};

rf = rf0/sum(rf0(:));

end