Skip to content
GitLab
Commit f13ff8ab authored by jml1g18's avatar jml1g18
Browse files

first commit of 2D PIV code

parents
Hide whitespace changes
Inline Side-by-side
PIV_2D.m 0 → 100644
View file @ f13ff8ab
function piv_result = PIV_2D(A, B, piv_setup)
% PIV_2D(A, B, piv_setup)
%
% basic 2D PIV cross-correlation of images A and B
% which must be the same size. the result is returned in piv_result
% only a single pass of PIV cross-correlation is performed
% images must be dewarped prior to
%
% piv_setup is a struct containing the basics of the PIV setup
% its members are:
%
% im_x 1 x Nx vector vector describing x axis of image
% im_y 1 x Ny vector vector describing y axis of image
% wsize 1 x 2 vector window size, in pixels
% overlap 1 x 2 vector window overlap, as a percentage
% peak_finder string character string describing peak
% finding routine
% weight_fun string character string describing weighting
% function
% dt scalar PIV dt, default is 1
% median_filter boolean boolean indicating whether or not to
% apply normalised median filtering (see
% Raffel 1998)
% epsilon_0 scalar typical PIV error, parameter for median
% filtering, usually 0.1-0.2 px
% epsilon_thresh scalar normalised threshold for vector
% rejection using median filter
%
% piv_result is a struct containing the result of 2D PIV
% cross-correlation. its members are:
%
% n_windows 1 x 2 vector number of PIV cross-correlation windows
% in either direction (x and y)
% win_x 1 x n_windows(1) x-axis of grid for PIV window centres
% win_y 1 x n_windows(2) y-axis of grid for PIV window centres
% ux matrix size n_windows(1) x n_windows(2)
% x-component of displacement
% uy matrix y-component of displacement
% nan_mask bool matrix boolean matrix identifying bad vectors
% peak_mag matrix chosen peak magnitude from peak fitting
% Q matrix Q factor from peak fitting
%
%% convert images to single precision
A = single(A);
B = single(B);
%% calculate grid for window centres
im_x = piv_setup.im_x;
im_y = piv_setup.im_y;
wsize = piv_setup.wsize;
overlap = piv_setup.overlap;
im_dx = im_x(2) - im_x(1);
im_dy = im_y(2) - im_y(1);
piv_dt = piv_setup.dt;
Nx = length(im_x);
Ny = length(im_y);
win_axis_x = -wsize(1)/2+0.5 : +wsize(1)/2-0.5;
win_axis_y = -wsize(2)/2+0.5 : +wsize(2)/2-0.5;
win_spacing_x = round((1 - overlap/100) * wsize(1));
win_spacing_y = round((1 - overlap/100) * wsize(2));
win_ctrs_x = 0.5 + wsize(1)/2 : win_spacing_x : Nx - wsize(1)/2 + 0.5;
win_ctrs_y = 0.5 + wsize(2)/2 : win_spacing_y : Ny - wsize(2)/2 + 0.5;
n_windows = [length(win_ctrs_x) length(win_ctrs_y)];
%% calculate window weighting function
% and correlation plane weighting function
% According to Raffel et al, the correlation plane weight factors
% can be obtained by convolving the image weighting function with itself
window_weight = window_weight_fun(wsize);
samplefun = zeros(wsize*2, 'single');
samplefun(wsize(1)/2 + (1:wsize(1)), wsize(2)/2 + (1:wsize(2))) ...
= window_weight;
weight = fastxcorr2(samplefun, samplefun) / prod(wsize);
weight = weight(wsize(1)/2 + (1:wsize(1)), wsize(2)/2 + (1:wsize(2)));
weight = weight / max(weight(:));
correl_weight = 1 ./ weight;
%% allocate memory for cross-correlation result
ux_mat = zeros(n_windows);
uy_mat = zeros(n_windows);
Q_mat = zeros(n_windows);
peak_mag_mat = zeros(n_windows);
nan_mask = false(n_windows);
xx_cut = zeros(wsize);
xx_xrange = wsize(1)/4 + (1 : wsize(1)/2);
xx_yrange = wsize(2)/4 + (1 : wsize(2)/2);
%% PIV cross-correlation
% translate integer codes for peak finder
peak_finder = 3;
if isfield(piv_setup, 'peak_finder')
if strcmp('3-point', piv_setup.peak_finder)
peak_finder = 3;
elseif strcmp('gauss4', piv_setup.peak_finder)
peak_finder = 4;
elseif strcmp('gauss5', piv_setup.peak_finder)
peak_finder = 5;
elseif strcmp('gauss6', piv_setup.peak_finder)
peak_finder = 6;
end
end
% perform cross correlation
[peak_loc_x, peak_loc_y, peak_height] = ...
PIV_2d_cross_correlate( A, B, ...
single(win_ctrs_x-1), single(win_ctrs_y-1), ...
single(wsize), single(window_weight), ...
1, peak_finder);
% format results
ux_mat = squeeze(peak_loc_x(1, :, :)); % first peak
uy_mat = squeeze(peak_loc_y(1, :, :));
peak_mag_mat = squeeze(peak_height(1,:, :));
Q = inf(n_windows);
%% inpaint NaNs
nan_mask = isnan(ux_mat) | isnan(uy_mat);
if any(nan_mask(:))
ux_mat = inpaint_nans(double(ux_mat));
uy_mat = inpaint_nans(double(uy_mat));
end
%% normalised median test filter
if piv_setup.median_filter
% apply median filter to velocity field
% the median filter replaces each vector with the median of its
% eight neighbours
ux_med = medfilt2(ux_mat, [3 3], 'symmetric');
uy_med = medfilt2(uy_mat, [3 3], 'symmetric');
% calculate median residual between any given median filtered vector and its
% neighbours
ux_pad = padarray(ux_mat, [1 1], 'symmetric', 'both');
uy_pad = padarray(uy_mat, [1 1], 'symmetric', 'both');
r_mn = zeros([n_windows 3 3]);
for m = 1 : 3
for n = 1 : 3
r_mn(:,:,m,n) = (ux_med - ux_pad(m-1 + (1:n_windows(1)),n-1 + (1:n_windows(2)))).^2 ...
+ (uy_med - uy_pad(m-1 + (1:n_windows(1)),n-1 + (1:n_windows(2)))).^2;
end
end
r_mn = reshape(sqrt(r_mn), [n_windows 9]);
r_mn = r_mn(:,:,[1 2 3 4 6 7 8 9]);
r_med = median(r_mn, 3);
% calculate the norm of the residual between median
% and unfiltered vector field
r_2d = sqrt((ux_med - ux_mat).^2 + (uy_med - uy_mat).^2);
median_filter = r_2d ./ (r_med + piv_setup.epsilon_0) > piv_setup.epsilon_thresh;
% replace anomalous vectors with inpainted values
if any(median_filter(:))
ux_mat(median_filter) = nan;
uy_mat(median_filter) = nan;
ux_mat = inpaint_nans(double(ux_mat));
uy_mat = inpaint_nans(double(uy_mat));
end
end
%% prepare results structure for output
piv_result = struct('n_windows', n_windows, ...
'win_x', interp1(1:Nx, im_x, win_ctrs_x), ...
'win_y', interp1(1:Ny, im_y, win_ctrs_y), ...
'ux', ux_mat * im_dx / piv_dt, ...
'uy', uy_mat * im_dy / piv_dt, ...
'nan_mask', nan_mask, ...
'Q', Q_mat, ...
'peak_mag', peak_mag_mat);
end
%% cross-correlation weighting window
function [weight] = window_weight_fun(window_size)
% Calculates the cross-correlation weighting window
m = window_size(1);
n = window_size(2);
dist_to_cen_max = ( ((m+1)/2)^2 + ((n+1)/2)^2 )^0.5;
weight = zeros(m, n);
% coefficients for blackman window
a0 = 0.42659;
a1 = 0.49656;
a2 = 0.076849;
for ii = 1 : m
for jj = 1 : n
%% blackman window
dist_to_cen = ( ((ii-((m+1)/2))^2) + ((jj-((n+1)/2))^2) )^0.5;
theta = dist_to_cen / dist_to_cen_max * pi + pi;
weight(ii,jj) = a0 - a1*cos(theta) + a2*cos(2*theta);
end
end
% delete very small coefficients
weight(weight <= 0) = 0.01;
end
PIV_2D_matlab.m 0 → 100644
View file @ f13ff8ab
function piv_result = PIV_2D(A, B, piv_setup)
% PIV_2D(A, B, piv_setup)
%
% basic 2D PIV cross-correlation of images A and B
% which must be the same size. the result is returned in piv_result
% only a single pass of PIV cross-correlation is performed
% images must be dewarped prior to
%
% piv_setup is a struct containing the basics of the PIV setup
% its members are:
%
% im_x 1 x Nx vector vector describing x axis of image
% im_y 1 x Ny vector vector describing y axis of image
% wsize 1 x 2 vector window size, in pixels
% overlap 1 x 2 vector window overlap, as a percentage
% peak_finder string character string describing peak
% finding routine
% weight_fun string character string describing weighting
% function
% dt scalar PIV dt, default is 1
% median_filter boolean boolean indicating whether or not to
% apply normalised median filtering (see
% Raffel 1998)
% epsilon_0 scalar typical PIV error, parameter for median
% filtering, usually 0.1-0.2 px
% epsilon_thresh scalar normalised threshold for vector
% rejection using median filter
%
% piv_result is a struct containing the result of 2D PIV
% cross-correlation. its members are:
%
% n_windows 1 x 2 vector number of PIV cross-correlation windows
% in either direction (x and y)
% win_x 1 x n_windows(1) x-axis of grid for PIV window centres
% win_y 1 x n_windows(2) y-axis of grid for PIV window centres
% ux matrix size n_windows(1) x n_windows(2)
% x-component of displacement
% uy matrix y-component of displacement
% nan_mask bool matrix boolean matrix identifying bad vectors
% peak_mag matrix chosen peak magnitude from peak fitting
% Q matrix Q factor from peak fitting
%
%% convert images to single precision
A = single(A);
B = single(B);
%% calculate grid for window centres
im_x = piv_setup.im_x;
im_y = piv_setup.im_y;
wsize = piv_setup.wsize;
overlap = piv_setup.overlap;
im_dx = im_x(2) - im_x(1);
im_dy = im_y(2) - im_y(1);
piv_dt = piv_setup.dt;
Nx = length(im_x);
Ny = length(im_y);
win_axis_x = -wsize(1)/2+0.5 : +wsize(1)/2-0.5;
win_axis_y = -wsize(2)/2+0.5 : +wsize(2)/2-0.5;
win_spacing_x = round((1 - overlap/100) * wsize(1));
win_spacing_y = round((1 - overlap/100) * wsize(2));
win_ctrs_x = 0.5 + wsize(1)/2 : win_spacing_x : Nx - wsize(1)/2 + 0.5;
win_ctrs_y = 0.5 + wsize(2)/2 : win_spacing_y : Ny - wsize(2)/2 + 0.5;
n_windows = [length(win_ctrs_x) length(win_ctrs_y)];
%% calculate window weighting function
% and correlation plane weighting function
% According to Raffel et al, the correlation plane weight factors
% can be obtained by convolving the image weighting function with itself
window_weight = window_weight_fun(wsize);
samplefun = zeros(wsize*2, 'single');
samplefun(wsize(1)/2 + (1:wsize(1)), wsize(2)/2 + (1:wsize(2))) ...
= window_weight;
weight = fastxcorr2(samplefun, samplefun) / prod(wsize);
weight = weight(wsize(1)/2 + (1:wsize(1)), wsize(2)/2 + (1:wsize(2)));
weight = weight / max(weight(:));
correl_weight = 1 ./ weight;
%% allocate memory for cross-correlation result
ux_mat = zeros(n_windows);
uy_mat = zeros(n_windows);
Q_mat = zeros(n_windows);
peak_mag_mat = zeros(n_windows);
nan_mask = false(n_windows);
xx_cut = zeros(wsize);
xx_xrange = wsize(1)/4 + (1 : wsize(1)/2);
xx_yrange = wsize(2)/4 + (1 : wsize(2)/2);
%% PIV cross-correlation
% translate integer codes for peak finder
peak_finder = 3;
if isfield(piv_setup, 'peak_finder')
if strcmp('3-point', piv_setup.peak_finder)
peak_finder = 3;
elseif strcmp('gauss4', piv_setup.peak_finder)
peak_finder = 4;
elseif strcmp('gauss5', piv_setup.peak_finder)
peak_finder = 5;
elseif strcmp('gauss6', piv_setup.peak_finder)
peak_finder = 6;
end
end
% perform cross correlation
[peak_loc_x, peak_loc_y, peak_height] = ...
PIV_2d_cross_correlate( A, B, ...
single(win_ctrs_x-1), single(win_ctrs_y-1), ...
single(wsize), single(window_weight), ...
1, peak_finder);
% format results
ux_mat = squeeze(peak_loc_x(1, :, :)); % first peak
uy_mat = squeeze(peak_loc_y(1, :, :));
peak_mag_mat = squeeze(peak_height(1,:, :));
Q = inf(n_windows);
for ix = 1 : n_windows(1)
for iy = 1 : n_windows(2)
%% get weighted correlation window
winA = A(win_axis_x+win_ctrs_x(ix), win_axis_y+win_ctrs_y(iy));
winB = B(win_axis_x+win_ctrs_x(ix), win_axis_y+win_ctrs_y(iy));
% apply weighting
winA = winA .* window_weight;
winB = winB .* window_weight;
%% cross correlate
% normalised cross correlation
winA = winA - mean(winA(:));
winB = winB - mean(winB(:));
winA_en = sum(winA(:).^2);
winB_en = sum(winB(:).^2);
xx = xcorr2(winB, winA);
xx = xx(wsize(1)/2 + (1:wsize(1)), wsize(2)/2 + (1:wsize(2)));
xx = xx / sqrt(winA_en * winB_en);
% weighted correlation matrix helps remove bias
xx = xx .* correl_weight;
%% peak find
xx_cut(xx_xrange, xx_yrange) = xx(xx_xrange, xx_yrange); % remove edges
peak_mag = max(xx_cut(:));
peak_loc = PIV_2d_peaklocate(xx_cut);
peak_loc = peak_loc(:) - wsize(:) / 2;
%% save peak
ux_mat(ix, iy) = peak_loc(1);
uy_mat(ix, iy) = peak_loc(2);
peak_mag_mat(ix,iy) = peak_mag;
Q_mat(ix, iy) = inf;
end
end
%% inpaint NaNs
nan_mask = isnan(ux_mat) | isnan(uy_mat);
if any(nan_mask(:))
ux_mat = inpaint_nans(double(ux_mat));
uy_mat = inpaint_nans(double(uy_mat));
end
%% normalised median test filter
if piv_setup.median_filter
% apply median filter to velocity field
% the median filter replaces each vector with the median of its
% eight neighbours
ux_med = medfilt2(ux_mat, [3 3], 'symmetric');
uy_med = medfilt2(uy_mat, [3 3], 'symmetric');
% calculate median residual between any given median filtered vector and its
% neighbours
ux_pad = padarray(ux_mat, [1 1], 'symmetric', 'both');
uy_pad = padarray(uy_mat, [1 1], 'symmetric', 'both');
r_mn = zeros([n_windows 3 3]);
for m = 1 : 3
for n = 1 : 3
r_mn(:,:,m,n) = (ux_med - ux_pad(m-1 + (1:n_windows(1)),n-1 + (1:n_windows(2)))).^2 ...
+ (uy_med - uy_pad(m-1 + (1:n_windows(1)),n-1 + (1:n_windows(2)))).^2;
end
end
r_mn = reshape(sqrt(r_mn), [n_windows 9]);
r_mn = r_mn(:,:,[1 2 3 4 6 7 8 9]);
r_med = median(r_mn, 3);
% calculate the norm of the residual between median
% and unfiltered vector field
r_2d = sqrt((ux_med - ux_mat).^2 + (uy_med - uy_mat).^2);
median_filter = r_2d ./ (r_med + piv_setup.epsilon_0) > piv_setup.epsilon_thresh;
% replace anomalous vectors with inpainted values
if any(median_filter(:))
ux_mat(median_filter) = nan;
uy_mat(median_filter) = nan;
ux_mat = inpaint_nans(double(ux_mat));
uy_mat = inpaint_nans(double(uy_mat));
end
end
%% prepare results structure for output
piv_result = struct('n_windows', n_windows, ...
'win_x', interp1(1:Nx, im_x, win_ctrs_x), ...
'win_y', interp1(1:Ny, im_y, win_ctrs_y), ...
'ux', ux_mat * im_dx / piv_dt, ...
'uy', uy_mat * im_dy / piv_dt, ...
'nan_mask', nan_mask, ...
'Q', Q_mat, ...
'peak_mag', peak_mag_mat);
end
%% cross-correlation weighting window
function [weight] = window_weight_fun(window_size)
% Calculates the cross-correlation weighting window
m = window_size(1);
n = window_size(2);
dist_to_cen_max = ( ((m+1)/2)^2 + ((n+1)/2)^2 )^0.5;
weight = zeros(m, n);
% coefficients for blackman window
a0 = 0.42659;
a1 = 0.49656;
a2 = 0.076849;
for ii = 1 : m
for jj = 1 : n
%% blackman window
dist_to_cen = ( ((ii-((m+1)/2))^2) + ((jj-((n+1)/2))^2) )^0.5;
theta = dist_to_cen / dist_to_cen_max * pi + pi;
weight(ii,jj) = a0 - a1*cos(theta) + a2*cos(2*theta);
end
end
% delete very small coefficients
weight(weight <= 0) = 0.01;
end
PIV_2D_outlier.m 0 → 100644
View file @ f13ff8ab
function b_filter = PIV_2D_outlier(ux, uy, epsilon, threshold)
% Universal outlier detection for 2D PIV data
% based on normalised median outlier test in [1]
%
% [1] Jerry Westerweel & Fulvio Scarano, Universal outlier detection
% for PIV data (2005), Experiments In Fluids, 39 1096-1100
n_windows = size(ux);
ui = cat(3, ux, uy);
%% iterate over components
r_0p = zeros([n_windows, 2]);
n_neighbours= zeros([n_windows, 2]);
for c = 1 : 2
% get this component
U = ui(:,:,c);
% fill out neighbours
U_pad = padarray(U, [1 1], nan, 'both');
U_nn = zeros([n_windows 8]);
U_nn(:,:,1) = U_pad(1:end-2, 1:end-2);
U_nn(:,:,2) = U_pad(1:end-2, 2:end-1);
U_nn(:,:,3) = U_pad(1:end-2, 3:end );
U_nn(:,:,4) = U_pad(2:end-1, 1:end-2);
U_nn(:,:,5) = U_pad(2:end-1, 3:end );
U_nn(:,:,6) = U_pad(3:end , 3:end );
U_nn(:,:,7) = U_pad(3:end , 2:end-1);
U_nn(:,:,8) = U_pad(3:end , 1:end-2);
% median of neighbours
U_med = nanmedian(U_nn, 3);
% median residual of neighbours
r_0 = abs(U_med - U);
r_i = abs(U_med - U_nn);
r_m = nanmedian(r_i, 3);
% measure number of valid neighbours
n_neighbours(:,:,c) = conv2(~isnan(U), ones(3,3), 'same');
% measure normalised residual, equation (2)
r_0p(:,:,c) = r_0 ./ (r_m + epsilon);
end
%% apply detection criterion
r_0_combined= max(r_0p, [], 3);
b_filter = r_0_combined > threshold | any(isnan(r_0p), 3) | any(n_neighbours<6,3);
end
\ No newline at end of file
PIV_2D_wdef.m 0 → 100644
View file @ f13ff8ab
function piv_result = PIV_2D_wdef(A, B, mask, piv_setup)
% PIV_2D(A, B, piv_setup)
%
% basic 2D PIV cross-correlation of images A and B
% which must be the same size. the result is returned in piv_result
% only a single pass of PIV cross-correlation is performed
% images must be dewarped prior to
%
% piv_setup is a struct containing the basics of the PIV setup
% its members are:
%
% im_x 1 x Nx vector vector describing x axis of image
% im_y 1 x Ny vector vector describing y axis of image
% wsize 1 x 2 vector window size, in pixels
% overlap 1 x 2 vector window overlap, as a percentage
% peak_finder string character string describing peak
% finding routine
% weight_fun string character string describing weighting
% function
% dt scalar PIV dt, default is 1
% median_filter boolean boolean indicating whether or not to
% apply normalised median filtering (see
% Raffel 1998)
% epsilon_0 scalar typical PIV error, parameter for median
% filtering, usually 0.1-0.2 px
% epsilon_thresh scalar normalised threshold for vector
% rejection using median filter
% mask_thresh scalar maximum acceptable fraction of masked pixels
% within an interrogation window
%
% piv_result is a struct containing the result of 2D PIV
% cross-correlation. its members are:
%
% n_windows 1 x 2 vector number of PIV cross-correlation windows
% in either direction (x and y)
% win_x 1 x n_windows(1) x-axis of grid for PIV window centres
% win_y 1 x n_windows(2) y-axis of grid for PIV window centres
% ux matrix size n_windows(1) x n_windows(2)
% x-component of displacement
% uy matrix y-component of displacement
% nan_mask bool matrix boolean matrix identifying bad vectors
% peak_mag matrix chosen peak magnitude from peak fitting
% Q matrix Q factor from peak fitting
%
%% convert images to single precision
A = single(A);
B = single(B);
%% mask out invalid regions
A(mask) = 0;
B(mask) = 0;
%% extract main parameters
n_passes = piv_setup.n_passes;
n_peaks = piv_setup.n_peaks;
im_x = piv_setup.im_x;
im_y = piv_setup.im_y;
im_dx = im_x(2) - im_x(1);
im_dy = im_y(2) - im_y(1);
piv_dt = piv_setup.dt;
Nx = length(im_x);
Ny = length(im_y);
im_size = size(A);
im_px = prod(im_size);
im_i = 1 : im_size(1);
im_j = 1 : im_size(2);
[im_imat,im_jmat]=ndgrid(im_i, im_j);
im_mesh = single(cat(3, im_imat, im_jmat));
% image interpolation for window deformation
wdef_ksize = 4;
wdef_kernel = 'lanczos';
%% pre-allocate results structure
empty = cell(n_passes,1);
piv_result = struct( 'n_windows', empty, ...
'win_x', empty, 'win_y', empty, ...
'ux', empty, 'uy', empty, ...
'nan_mask', empty, 'Q', empty, ...
'peak_mag', empty, 'peak_choice', empty);
%% ITERATE OVER PASSES
for pass = 1 : n_passes
%% report
fprintf('PASS %d of %d\n', pass, n_passes);
%% calculate grid for window centres in images A and B
%
% /---/ +---+ \---\
% / / -> | | -> \ \
% /---/ +---+ \---\
% A (0) B
%
% the deformation between A and B is symmetric
% one supposes a pair of identical deformations that maps
% A forwards and B backwards onto the hypothetical (0), which is cuboidal
%
% a window is defined by:
% - the regular mesh of points in the window
% - the deformation matrix applied to the window
wsize = piv_setup.wsize(pass, :);
wtype = piv_setup.wtype{pass};
overlap = piv_setup.overlap(pass);
%win_axis_x = -wsize(1)/2+0.5 : +wsize(1)/2-0.5;
%win_axis_y = -wsize(2)/2+0.5 : +wsize(2)/2-0.5;
win_spacing_x = round((1 - overlap/100) * wsize(1));
win_spacing_y = round((1 - overlap/100) * wsize(2));
win_ctrs_x = 0.5 + wsize(1)/2 : win_spacing_x : Nx - wsize(1)/2 + 0.5;
win_ctrs_y = 0.5 + wsize(2)/2 : win_spacing_y : Ny - wsize(2)/2 + 0.5;
[win_x,win_y] = ndgrid(win_ctrs_x, win_ctrs_y);
n_windows = [length(win_ctrs_x) length(win_ctrs_y)];
% calculate window weighting function
window_weight = window_weight_fun(wsize, wtype);
window_weight_norm = window_weight / sum(sum(window_weight));
%% build mask
% vectors are masked out if there are too few unmasked pixels in
% an interrogation window
% TODO: speed up this routine, it takes just as long as
% cross-correlation!
f_mask = conv2( mask , ones(wsize(1),1)/wsize(1), 'same');
f_mask = conv2(f_mask.', ones(wsize(2),1)/wsize(2), 'same').';
b_mask = interpn(im_x, im_y, f_mask, win_x, win_y, 'nearest') > piv_setup.mask_thresh;
%% build predictor displacement field and apply window deformation