# DPIV_ALGORITHM
import numpy as np
from scipy import ndimage
from scipy import interpolate
from scipy.io import savemat
import copy
import cv2
import time
from dpivsoft.Classes import Parameters
from dpivsoft.Classes import GPU
from dpivsoft.Classes import grid
[docs]
def processing(Img1, Img2):
"""
Run the complete two-pass PIV algorithm on an image pair
(CPU implementation). Processing options are taken from the
global Parameters class.
Parameters
----------
Img1, Img2 : 2d np.ndarray
First and second frames of the PIV pair.
Returns
-------
x_2, y_2 : 2d np.ndarray
Coordinates of the final (second-pass) grid in pixels.
u_2, v_2 : 2d np.ndarray
Displacement field on the final grid in pixels.
"""
# Generate x-y mesh for the PIV
grid.generate_mesh(Img1.shape[1], Img1.shape[0])
# Apply mask to images
if Parameters.mask:
Img1, Img2 = masking(Img1, Img2)
# First cross-correlation
if Parameters.direct_calc:
# direct cross correlation
x1, y1, u1, v1 = corrDirect1(Img1, Img2)
else:
# FFT cross correlation
x1, y1, u1, v1 = corrFFT1(Img1, Img2)
# Iterate on the first grid if specified
if Parameters.no_iter_1 > 1:
x1, y1, u1, v1 = corrFFT1bis(Img1, Img2, x1, y1, u1, v1)
# Second cross-correlation
x2, y2, u2, v2 = corrFFT2(Img1, Img2, x1, y1, u1, v1)
return x2, y2, u2, v2
[docs]
def corrFFT1(Img1, Img2):
"""
First PIV pass: FFT-based cross-correlation on the coarse grid,
without window deformation.
Parameters
----------
Img1, Img2 : 2d np.ndarray
Image pair.
Returns
-------
x_1, y_1 : 2d np.ndarray
First-pass grid coordinates in pixels.
u_1, v_1 : 2d np.ndarray
Displacement field on the first-pass grid in pixels.
"""
# Definition of Parameters to reduce length
box_size_x = Parameters.box_size_1_x
box_size_y = Parameters.box_size_1_y
no_boxes_1_x = Parameters.no_boxes_1_x
no_boxes_1_y = Parameters.no_boxes_1_y
window_x = Parameters.window_1_x
window_y = Parameters.window_1_y
Height, Width = Img1.shape
Parameters.Data.height = Height
Parameters.Data.width = Width
# Initialize all matrices
box_origin_x_1 = grid.box_origin_x_1
box_origin_y_1 = grid.box_origin_y_1
x_1 = grid.x_1
y_1 = grid.y_1
u_1 = np.zeros([no_boxes_1_y, no_boxes_1_x])
v_1 = np.zeros([no_boxes_1_y, no_boxes_1_x])
# Prevent the search window from being larger than the box
window_x = np.amin([2 * np.round(window_x / 2), box_size_x - 4])
window_y = np.amin([2 * np.round(window_y / 2), box_size_y - 4])
if Parameters.weighting:
i_matrix, j_matrix = np.meshgrid(np.arange(box_size_x),
np.arange(box_size_y))
Weighting_Function = weight_function(i_matrix, j_matrix,
box_size_x, box_size_y)
# Apply Gaussian filter to images only for first iteration
if Parameters.gaussian_size:
Img1, Img2 = gaussian_filter(Img1, Img2, Parameters.gaussian_size)
for i in range(0, no_boxes_1_x):
for j in range(0, no_boxes_1_y):
# Define sub images to work
box_o_y = int(box_origin_y_1[j, i])
box_o_x = int(box_origin_x_1[j, i])
SubImg1 = (Img1[box_o_y: box_o_y + box_size_y,
box_o_x: box_o_x + box_size_x])
SubImg2 = (Img2[box_o_y: box_o_y + box_size_y,
box_o_x: box_o_x + box_size_x])
# Replace masked pixels by the unmasked mean intensity
if Parameters.mask:
SubImg1, SubImg2 = change_mask(SubImg1, SubImg2)
if Parameters.weighting:
SubImg1 = np.multiply(SubImg1, Weighting_Function)
SubImg2 = np.multiply(SubImg2, Weighting_Function)
# Center the image intensity on its mean
SubImg1 = SubImg1 - np.sum(SubImg1) / (box_size_y * box_size_x)
SubImg2 = SubImg2 - np.sum(SubImg2) / (box_size_y * box_size_x)
# Intensity norms used to normalize the correlation
Sigma1 = max(0.1, np.sqrt(np.sum(SubImg1**2)))
Sigma2 = max(0.1, np.sqrt(np.sum(SubImg2**2)))
# Cross-correlation of the image pair, normalized by
# Sigma1 and Sigma2
correlation = (np.fft.fftshift(np.real(np.fft.ifft2(
np.multiply(np.conj(np.fft.fft2(SubImg1)),
np.fft.fft2(SubImg2))))) / (Sigma1 * Sigma2))
# Find peaks of the correlation function
epsilon_x, epsilon_y, col_idx, row_idx = find_peaks(
correlation, window_x, window_y)
u_1[j, i] = epsilon_x + col_idx - box_size_x / 2
v_1[j, i] = epsilon_y + row_idx - box_size_y / 2
return x_1, y_1, u_1, v_1
[docs]
def corrFFT1bis(Img1, Img2, x_1, y_1, u_1, v_1):
"""
Optional extra iterations of the first pass (Parameters.no_iter_1 > 1):
median filtering plus window deformation, refining the corrFFT1
result on the same grid.
Parameters
----------
Img1, Img2 : 2d np.ndarray
Image pair.
x_1, y_1, u_1, v_1 : 2d np.ndarray
Grid and displacement field from corrFFT1.
Returns
-------
x_1, y_1, u_1, v_1 : 2d np.ndarray
Refined displacement field on the same first-pass grid.
"""
# Definition of Parameters to reduce length
box_size_x = Parameters.box_size_1_x
box_size_y = Parameters.box_size_1_y
no_boxes_1_x = Parameters.no_boxes_1_x
no_boxes_1_y = Parameters.no_boxes_1_y
window_x = Parameters.window_1_x
window_y = Parameters.window_1_y
median_limit = Parameters.median_limit
Height, Width = Img1.shape
# Define origin of boxes
box_origin_x = x_1 - box_size_x / 2
box_origin_y = y_1 - box_size_y / 2
i_matrix, j_matrix = np.meshgrid(np.arange(box_size_x),
np.arange(box_size_y))
if Parameters.weighting:
Weighting_Function = weight_function(i_matrix, j_matrix,
box_size_x, box_size_y)
# Prevent the search window from being larger than the box
window_x = np.amin([2 * np.round(window_x / 2), box_size_x - 4])
window_y = np.amin([2 * np.round(window_y / 2), box_size_y - 4])
for calc in range(1, Parameters.no_iter_1):
# Median filter
u_1, v_1, err_vect = median_filter(u_1, v_1, median_limit)
# If masked, zero vectors inside the mask to prevent bleeding
if Parameters.mask:
u_1, v_1 = check_mask(u_1, v_1, grid.mask_1)
# Compute the velocity gradients (Jacobian) on the grid
du_dx, du_dy, dv_dx, dv_dy = jacobian_matrix(u_1, v_1, x_1, y_1,
no_boxes_1_x, no_boxes_1_y)
for j in range(0, no_boxes_1_y):
for i in range(0, no_boxes_1_x):
# Obtain deformed image.
SubImg1, SubImg2, u_index, v_index = deform_image(
Img1, Img2, Width, Height, box_origin_x, box_origin_y,
i_matrix, j_matrix, box_size_x, box_size_y, u_1, v_1,
du_dx, du_dy, dv_dx, dv_dy, i, j)
# Replace masked pixels by the unmasked mean intensity
if Parameters.mask:
SubImg1, SubImg2 = change_mask(SubImg1, SubImg2)
# Weighting if required
if Parameters.weighting:
SubImg1 = np.multiply(SubImg1, Weighting_Function)
SubImg2 = np.multiply(SubImg2, Weighting_Function)
# Intensity norms used to normalize the correlation
Sigma1 = max(0.1, np.sqrt(np.sum(SubImg1**2)))
Sigma2 = max(0.1, np.sqrt(np.sum(SubImg2**2)))
# Cross-correlation of the image pair, normalized by
# Sigma1 and Sigma2
correlation = (np.fft.fftshift(np.abs(np.fft.ifft2(
np.multiply(np.conj(np.fft.fft2(SubImg1)),
np.fft.fft2(SubImg2))))) / (Sigma1 * Sigma2))
# Find peaks of the correlation function
epsilon_x, epsilon_y, col_idx, row_idx = find_peaks(
correlation, window_x, window_y)
u_1[j, i] = (u_index[row_idx, col_idx] + epsilon_x +
col_idx - box_size_x / 2)
v_1[j, i] = (v_index[row_idx, col_idx] + epsilon_y +
row_idx - box_size_y / 2)
return x_1, y_1, u_1, v_1
[docs]
def corrFFT2(Img1, Img2, x_1, y_1, u_1, v_1):
"""
Second PIV pass: interpolate the first-pass result onto the finer
second grid, then iterate FFT cross-correlation with window
deformation until sub-pixel convergence (at most
Parameters.no_iter_2 iterations per point).
Parameters
----------
Img1, Img2 : 2d np.ndarray
Image pair.
x_1, y_1, u_1, v_1 : 2d np.ndarray
Grid and displacement field from the first pass.
Returns
-------
x_2, y_2 : 2d np.ndarray
Second-pass grid coordinates in pixels.
u_2, v_2 : 2d np.ndarray
Displacement field on the second-pass grid in pixels.
"""
# Definition of Parameters to reduce length
no_boxes_1_x = Parameters.no_boxes_1_x
no_boxes_1_y = Parameters.no_boxes_1_y
box_size_2_x = Parameters.box_size_2_x
box_size_2_y = Parameters.box_size_2_y
no_boxes_2_x = Parameters.no_boxes_2_x
no_boxes_2_y = Parameters.no_boxes_2_y
window_x = Parameters.window_2_x
window_y = Parameters.window_2_y
median_limit = Parameters.median_limit
Height, Width = Img1.shape
# Define index matrices
i_matrix, j_matrix = np.meshgrid(np.arange(0, box_size_2_x),
np.arange(0, box_size_2_y))
if Parameters.weighting:
Weighting_Function = weight_function(i_matrix, j_matrix,
box_size_2_x, box_size_2_y)
# Prevent the search window from being larger than the box
window_x = np.amin([2 * np.round(window_x / 2), box_size_2_x - 4])
window_y = np.amin([2 * np.round(window_y / 2), box_size_2_y - 4])
# Second grid is placed completely inside first one
x_2 = grid.x_2[0, :]
y_2 = grid.y_2[:, 0]
# Calculate Jacobian Matrix
du_dx, du_dy, dv_dx, dv_dy = jacobian_matrix(u_1, v_1, x_1, y_1,
no_boxes_1_x, no_boxes_1_y)
# Interpolate First Run Results on second grid
du_dx, du_dy, dv_dx, dv_dy, u_2, v_2, x_2, y_2 = interpolations(
du_dx, du_dy, dv_dx, dv_dy, u_1, v_1, x_1, y_1, x_2, y_2,
no_boxes_1_x * no_boxes_1_y)
# Define origin of boxes without translation
box_origin_x_2 = x_2 - box_size_2_x / 2
box_origin_y_2 = y_2 - box_size_2_y / 2
for j in range(0, no_boxes_2_y):
for i in range(0, no_boxes_2_x):
k = 0
epsilon_x = 1
epsilon_y = 1
while ((np.abs(epsilon_x > 0.5) or np.abs(epsilon_y > 0.5)) and
k < Parameters.no_iter_2):
k = k + 1
# Obtain deformed image.
SubImg1, SubImg2, u_index, v_index = deform_image(
Img1, Img2, Width, Height, box_origin_x_2, box_origin_y_2,
i_matrix, j_matrix, box_size_2_x, box_size_2_y,
u_2, v_2, du_dx, du_dy, dv_dx, dv_dy, i, j)
# Replace masked pixels by the unmasked mean intensity
if Parameters.mask:
SubImg1, SubImg2 = change_mask(SubImg1, SubImg2)
# Weighting if required
if Parameters.weighting:
SubImg1 = np.multiply(SubImg1, Weighting_Function)
SubImg2 = np.multiply(SubImg2, Weighting_Function)
# Intensity norms used to normalize the correlation
Sigma1 = max(0.1, np.sqrt(np.sum(SubImg1**2)))
Sigma2 = max(0.1, np.sqrt(np.sum(SubImg2**2)))
# Cross-correlation of the image pair, normalized by
# Sigma1 and Sigma2
correlation = (np.fft.fftshift(np.abs(np.fft.ifft2(
np.multiply(np.conj(np.fft.fft2(SubImg1)),
np.fft.fft2(SubImg2))))) / (Sigma1 * Sigma2))
# Find peaks of the correlation function
epsilon_x, epsilon_y, col_idx, row_idx = find_peaks(
correlation, window_x, window_y)
u_2[j, i] = (u_index[row_idx, col_idx] + epsilon_x +
col_idx - box_size_2_x / 2)
v_2[j, i] = (v_index[row_idx, col_idx] + epsilon_y +
row_idx - box_size_2_y / 2)
u_2, v_2, err_vect = median_filter(u_2, v_2, median_limit)
# If masked, zero vectors inside the mask to prevent bleeding
if Parameters.mask:
u_2, v_2 = check_mask(u_2, v_2, grid.mask_2)
return x_2, y_2, u_2, v_2
[docs]
def corrDirect1(Img1, Img2):
"""
First PIV pass using direct (spatial) cross-correlation instead of
FFT: the correlation is evaluated only inside the search window.
Parameters
----------
Img1, Img2 : 2d np.ndarray
Image pair.
Returns
-------
x_1, y_1 : 2d np.ndarray
First-pass grid coordinates in pixels.
u_1, v_1 : 2d np.ndarray
Displacement field on the first-pass grid in pixels.
box_origin_x_1 : 2d np.ndarray
x origin of each correlation box in pixels.
"""
# Definition of Parameters to reduce length
box_size_x = Parameters.box_size_1_x
box_size_y = Parameters.box_size_1_y
no_boxes_1_x = Parameters.no_boxes_1_x
no_boxes_1_y = Parameters.no_boxes_1_y
window_x = Parameters.window_x_1
window_y = Parameters.window_y_1
Height, Width = Img1.shape
# Initialize all matrices
box_origin_x_1 = np.zeros([no_boxes_1_y, no_boxes_1_x])
box_origin_y_1 = np.zeros([no_boxes_1_y, no_boxes_1_x])
x_1 = np.zeros([no_boxes_1_y, no_boxes_1_x])
y_1 = np.zeros([no_boxes_1_y, no_boxes_1_x])
u_1 = np.zeros([no_boxes_1_y, no_boxes_1_x])
v_1 = np.zeros([no_boxes_1_y, no_boxes_1_x])
correlation = np.zeros([window_y + 1, window_x + 1])
if Parameters.weighting:
i_matrix, j_matrix = np.meshgrid(np.arange(box_size_x),
np.arange(box_size_y))
Weighting_Function = weight_function(i_matrix, j_matrix,
box_size_x, box_size_y)
# Apply Gaussian filter to images only for first iteration
if Parameters.gaussian_size:
Img1, Img2 = gaussian_filter(Img1, Img2, Parameters.gaussian_size)
for j in range(0, no_boxes_1_y):
for i in range(0, no_boxes_1_x):
x_1[j, i] = (1 + round((window_x / 2 + box_size_x) / 2) +
round((i) * (Width - box_size_x - window_x / 2 - 4) / (no_boxes_1_x - 1)))
y_1[j, i] = (1 + round((window_y / 2 + box_size_y) / 2) +
round((j) * (Height - box_size_y - window_y / 2 - 4) / (no_boxes_1_y - 1)))
box_origin_x_1[j, i] = x_1[j, i] + 1 - round(box_size_x / 2)
box_origin_y_1[j, i] = y_1[j, i] + 1 - round(box_size_y / 2)
for jj in range(-round(window_y / 2), round(window_y / 2) + 1):
for ii in range(-round(window_x / 2), round(window_x / 2) + 1):
box_o_y = int(box_origin_y_1[j, i] + round(jj / 2))
box_o_x = int(box_origin_x_1[j, i] + round(ii / 2))
SubImg1 = (Img1[box_o_y - jj:box_o_y - jj + box_size_y,
box_o_x - ii:box_o_x - ii + box_size_x])
SubImg2 = (Img2[box_o_y:box_o_y + box_size_y,
box_o_x:box_o_x + box_size_x])
if Parameters.weighting:
SubImg1 = np.multiply(SubImg1, Weighting_Function)
SubImg2 = np.multiply(SubImg2, Weighting_Function)
# Center the image intensity on its mean
SubImg1 = SubImg1 - np.sum(SubImg1) / (box_size_y * box_size_x)
SubImg2 = SubImg2 - np.sum(SubImg2) / (box_size_y * box_size_x)
SubImg1 = np.divide(SubImg1, max(0.1,
np.sqrt(np.sum(SubImg1**2))))
SubImg2 = np.divide(SubImg2, max(0.1,
np.sqrt(np.sum(SubImg2**2))))
temp1 = jj + round(window_y / 2)
temp2 = ii + round(window_x / 2)
correlation[temp1, temp2] = (
np.sum(np.multiply(SubImg1, SubImg2))
)
# Find peaks of the correlation function
epsilon_x, epsilon_y, col_idx, row_idx = find_peaks(
correlation, window_x, window_y, 0)
u_1[j, i] = epsilon_x + col_idx - (window_x / 2)
v_1[j, i] = epsilon_y + row_idx - (window_y / 2)
return x_1, y_1, u_1, v_1, box_origin_x_1
[docs]
def gauss_subpixel(a, b, c):
"""Three-point Gaussian sub-pixel peak estimator.
Given the peak sample ``b`` and its two neighbours ``a`` (left/below) and
``c`` (right/above), return the sub-pixel offset of the peak relative to
``b``. Returns 0.0 when any sample is non-positive or the curvature is
degenerate (estimator undefined).
"""
if a <= 0 or b <= 0 or c <= 0:
return 0.0
denom = np.log(a) + np.log(c) - 2 * np.log(b)
if denom == 0:
return 0.0
return 0.5 * (np.log(a) - np.log(c)) / denom
[docs]
def find_peaks(correlation, window_x, window_y, westerweel=1,
peak_ratio=None, return_valid=False):
"""
Locate the correlation peak with sub-pixel accuracy.
The two highest peaks are found; the one with the larger surrounding
correlation sum is kept and refined with a three-point Gaussian
estimator (see gauss_subpixel).
Parameters
----------
correlation : 2d np.ndarray
Cross-correlation map of one interrogation window.
window_x, window_y : int
Size in pixels of the search window around the map center.
westerweel : int, optional
If nonzero (default), apply the Westerweel bias correction to
the peak neighbourhood before the sub-pixel fit.
peak_ratio : float, optional
Maximum allowed second/first peak-height ratio; defaults to
Parameters.peak_ratio.
return_valid : bool, optional
If True, also return whether the peak passed the peak-ratio
test. Used by the stereo disparity correction, which NaNs
ambiguous windows; the main pipeline leaves it False and relies
on median_filter instead.
Returns
-------
epsilon_x, epsilon_y : float
Sub-pixel offsets of the peak.
max_col, max_row : int
Integer peak position (column, row).
valid : bool
Only returned if return_valid is True.
"""
box_size_y, box_size_x = correlation.shape
# Default to the global setting; callers (e.g. disparity) may override.
if peak_ratio is None:
peak_ratio = Parameters.peak_ratio
# Indices of searching windows
ini_cor_y = int(np.round(box_size_y / 2 - window_y / 2))
end_cor_y = int(np.round(box_size_y / 2 + window_y / 2 + 1))
ini_cor_x = int(np.round(box_size_x / 2 - window_x / 2))
end_cor_x = int(np.round(box_size_x / 2 + window_x / 2 + 1))
# Find first peak
maxcor1 = (np.amax(correlation[ini_cor_y:end_cor_y, ini_cor_x:end_cor_x]))
max_row1, max_col1 = np.where(correlation == maxcor1)
max_row1 = max_row1[0].astype(int)
max_col1 = max_col1[0].astype(int)
if max_row1 == 0:
max_row1 = int(box_size_y / 2)
max_col1 = int(box_size_x / 2)
# Extract the 3x3 neighbourhood of the first peak
lm = max(1, int(np.round(box_size_x / 16)))
matmax1 = correlation[np.round(max_row1 - 1):np.round(max_row1 + 2),
np.round(max_col1 - 1):np.round(max_col1 + 2)]
correlation_z = copy.copy(correlation)
correlation_z[max_row1 - lm:max_row1 + lm + 1, max_col1 - lm:max_col1 + lm + 1] = 0
# Find second peak
maxcor2 = (np.amax(correlation_z[ini_cor_y:end_cor_y, ini_cor_x:end_cor_x]))
max_row2, max_col2 = np.where(correlation_z == maxcor2)
max_row2 = max_row2[0].astype(int)
max_col2 = max_col2[0].astype(int)
# Extract the 3x3 neighbourhood of the second peak
matmax2 = (correlation[np.round(max_row2 - 1):np.round(max_row2 + 2),
np.round(max_col2 - 1):np.round(max_col2 + 2)])
correlation_z[max_row2 - lm:max_row2 + lm + 1, max_col2 - lm:max_col2 + lm + 1] = 0
sum_cor_1 = np.sum(np.sum(matmax1)) - maxcor1
sum_cor_2 = np.sum(np.sum(matmax2)) - maxcor2
if sum_cor_1 > sum_cor_2:
fit_peak = matmax1
max_row = max_row1
max_col = max_col1
else:
fit_peak = matmax2
max_row = max_row2
max_col = max_col2
rejected = (maxcor2 / (maxcor1 + 1e-10) > peak_ratio
or fit_peak.shape != (3, 3))
if rejected:
fit_peak = np.zeros([3, 3])
# Define weighting functions
weight_i, weight_j = np.meshgrid(np.arange(-1, 2), np.arange(-1, 2))
# Westerweel bias correction
if westerweel:
fit_peak = (np.divide(np.divide(
fit_peak, 1 - np.abs(max_row - box_size_y / 2 + weight_j) / box_size_y),
1 - np.abs(max_col - box_size_x / 2 + weight_i) / box_size_x))
fit_peak[np.where(fit_peak < 0.001)] = 0.001
# Get sub-pixel accuracy with the Gaussian estimator
if fit_peak[1, 1] == np.amin(fit_peak[1, :]):
epsilon_x = 0
max_col = box_size_x // 2
else:
epsilon_x = gauss_subpixel(fit_peak[1, 0], fit_peak[1, 1], fit_peak[1, 2])
if fit_peak[1, 1] == np.amin(fit_peak[:, 1]):
epsilon_y = 0
max_row = box_size_y // 2
else:
epsilon_y = gauss_subpixel(fit_peak[0, 1], fit_peak[1, 1], fit_peak[2, 1])
if return_valid:
return epsilon_x, epsilon_y, max_col, max_row, not rejected
return epsilon_x, epsilon_y, max_col, max_row
[docs]
def gaussian_filter(Image_1, Image_2, gaussian_size):
"""
Smooth both images with the Gaussian kernel built by
gaussian_kernel. Used only on the first pass.
"""
B = gaussian_kernel(gaussian_size)
# Apply Gaussian filter to particles in the image
Image_1 = cv2.filter2D(Image_1, -1, B, borderType=cv2.BORDER_CONSTANT)
Image_2 = cv2.filter2D(Image_2, -1, B, borderType=cv2.BORDER_CONSTANT)
return Image_1, Image_2
[docs]
def gaussian_kernel(gaussian_size):
"""
Build the normalized 2d Gaussian convolution kernel of width
gaussian_size (pixels) used to pre-filter the images.
"""
x = np.arange(1, 2 * np.ceil(2 * gaussian_size).astype(int) + 2)
x = np.exp(-(x - np.ceil(2 * gaussian_size) - 1)**2 / gaussian_size**2)
B = np.outer(np.transpose(x), x)
B = B / np.sum(np.sum(B))
return B
[docs]
def jacobian_matrix(u, v, x, y, no_box_x, no_box_y):
"""
Compute the velocity gradients (du/dx, du/dy, dv/dx, dv/dy) on the
PIV grid: centered differences in the interior, one-sided at the
boundaries, smoothed with a 3x3 box filter.
Returns
-------
du_dx, du_dy, dv_dx, dv_dy : 2d np.ndarray
Velocity gradients on the same grid.
"""
# Initialize matrices
du_dy = np.zeros([no_box_y, no_box_x])
dv_dy = np.zeros([no_box_y, no_box_x])
du_dx = np.zeros([no_box_y, no_box_x])
dv_dx = np.zeros([no_box_y, no_box_x])
# Grid spacing
delta_x = x[0, 1] - x[0, 0]
delta_y = y[1, 0] - y[0, 0]
du_dy[1:no_box_y - 1, :] = (u[2:no_box_y, :] - u[0:no_box_y - 2, :]) / 2 / delta_y
dv_dy[1:no_box_y - 1, :] = (v[2:no_box_y, :] - v[0:no_box_y - 2, :]) / 2 / delta_y
du_dy[0, :] = (u[1, :] - u[0, :]) / delta_y
dv_dy[0, :] = (v[1, :] - v[0, :]) / delta_y
du_dy[no_box_y - 1, :] = (u[no_box_y - 1, :] - u[no_box_y - 2, :]) / delta_y
dv_dy[no_box_y - 1, :] = (v[no_box_y - 1, :] - v[no_box_y - 2, :]) / delta_y
du_dx[:, 1:no_box_x - 1] = (u[:, 2:no_box_x] - u[:, 0:no_box_x - 2]) / (2 * delta_x)
dv_dx[:, 1:no_box_x - 1] = (v[:, 2:no_box_x] - v[:, 0:no_box_x - 2]) / (2 * delta_x)
du_dx[:, 0] = (u[:, 1] - u[:, 0]) / delta_x
dv_dx[:, 0] = (v[:, 1] - v[:, 0]) / delta_x
du_dx[:, no_box_x - 1] = (u[:, no_box_x - 1] - u[:, no_box_x - 2]) / delta_x
dv_dx[:, no_box_x - 1] = (v[:, no_box_x - 1] - v[:, no_box_x - 2]) / delta_x
# Smooth the gradients with a 3x3 box filter
k = 1 / 9 * np.ones([3, 3]) # Convolution Kernel
du_dx = cv2.filter2D(du_dx, -1, k, borderType=cv2.BORDER_REPLICATE)
du_dy = cv2.filter2D(du_dy, -1, k, borderType=cv2.BORDER_REPLICATE)
dv_dx = cv2.filter2D(dv_dx, -1, k, borderType=cv2.BORDER_REPLICATE)
dv_dy = cv2.filter2D(dv_dy, -1, k, borderType=cv2.BORDER_REPLICATE)
return du_dx, du_dy, dv_dx, dv_dy
[docs]
def interpolations(du_dx, du_dy, dv_dx, dv_dy, u_1, v_1, x_1, y_1, x_2, y_2, no_box):
"""
Interpolate the first-pass velocity field and its gradients onto
the second-pass grid (linear scattered-data interpolation).
Returns
-------
du_dx, du_dy, dv_dx, dv_dy, u_2, v_2 : 2d np.ndarray
Fields interpolated onto the second grid.
x_2, y_2 : 2d np.ndarray
Second grid as meshgrid arrays.
"""
x_1 = x_1.reshape(no_box, order='F')
y_1 = y_1.reshape(no_box, order='F')
u_1 = u_1.reshape(no_box, order='F')
v_1 = v_1.reshape(no_box, order='F')
du_dx = du_dx.reshape(no_box, order='F')
du_dy = du_dy.reshape(no_box, order='F')
dv_dx = dv_dx.reshape(no_box, order='F')
dv_dy = dv_dy.reshape(no_box, order='F')
x_2, y_2 = np.meshgrid(x_2, y_2)
u_2 = interpolate.griddata((x_1, y_1), u_1, (x_2, y_2), method='linear')
v_2 = interpolate.griddata((x_1, y_1), v_1, (x_2, y_2), method='linear')
du_dx = interpolate.griddata((x_1, y_1), du_dx, (x_2, y_2), method='linear')
du_dy = interpolate.griddata((x_1, y_1), du_dy, (x_2, y_2), method='linear')
dv_dx = interpolate.griddata((x_1, y_1), dv_dx, (x_2, y_2), method='linear')
dv_dy = interpolate.griddata((x_1, y_1), dv_dy, (x_2, y_2), method='linear')
return du_dx, du_dy, dv_dx, dv_dy, u_2, v_2, x_2, y_2
[docs]
def translated_pixels(i_index, j_index, u_index, v_index, Width, Height,
box_size_x, box_size_y):
"""
Compute the symmetrically shifted (+-displacement/2) pixel positions
used to deform both sub-images. If a shift would fall outside the
image, the displacement is discarded (no shift) in that direction.
Returns
-------
i_frac_1, j_frac_1, i_frac_2, j_frac_2 : 2d np.ndarray
Fractional pixel offsets for the bilinear interpolation.
j_index_1, i_index_1, j_index_2, i_index_2 : 2d np.ndarray
Shifted pixel positions for sub-images 1 and 2.
"""
# Define position of translated pixel
if ((np.max(np.max(i_index + np.abs(u_index / 2))) < Width - 2) and
(np.min(np.min(i_index - np.abs(u_index / 2))) > 3)):
i_index_1 = i_index - u_index / 2
i_index_2 = i_index + u_index / 2
else:
i_index_1 = i_index
i_index_2 = i_index
u_index = np.zeros([box_size_y, box_size_x])
if ((np.max(np.max(j_index + np.abs(v_index / 2))) < Height - 2) and
(np.min(np.min(j_index - np.abs(v_index / 2))) > 3)):
j_index_1 = j_index - v_index / 2
j_index_2 = j_index + v_index / 2
else:
j_index_1 = j_index
j_index_2 = j_index
v_index = np.zeros([box_size_y, box_size_x])
# Define pixels position on translated Sub_image1
i_frac_1 = (i_index_1 - np.floor(i_index_1))
j_frac_1 = (j_index_1 - np.floor(j_index_1))
# Define pixels position on translated Sub_image2
i_frac_2 = (i_index_2 - np.floor(i_index_2))
j_frac_2 = (j_index_2 - np.floor(j_index_2))
return (i_frac_1, j_frac_1, i_frac_2, j_frac_2, j_index_1, i_index_1,
j_index_2, i_index_2)
[docs]
def weight_function(i_matrix, j_matrix, box_size_x, box_size_y):
"""
Separable quadratic weighting function that de-emphasizes pixels
near the edges of the interrogation window.
"""
# Define weighting function
zeta = i_matrix / box_size_x - 0.5 - 0.5 / box_size_x
eta = j_matrix / box_size_y - 0.5 - 0.5 / box_size_y
Weighting_Func = (9 * np.multiply(4 * np.power(zeta, 2) - 4 * np.abs(zeta) + 1,
4 * np.power(eta, 2) - 4 * np.abs(eta) + 1))
return Weighting_Func
[docs]
def masking(Img1, Img2):
"""Apply the binary mask (Parameters.Data.mask) to both images."""
Img1 = Parameters.Data.mask * Img1
Img2 = Parameters.Data.mask * Img2
return Img1, Img2
[docs]
def change_mask(SubImg1, SubImg2):
"""
Replace masked (zero) pixels of each sub-image by the mean
intensity of the unmasked pixels, so the mask does not bias the
correlation.
"""
unmasked_pixels = np.count_nonzero(SubImg1)
SubImg1[SubImg1 == 0] = sum(sum(SubImg1)) / max(unmasked_pixels, 1)
unmasked_pixels = np.count_nonzero(SubImg1)
SubImg2[SubImg2 == 0] = sum(sum(SubImg2)) / max(unmasked_pixels, 1)
return SubImg1, SubImg2
[docs]
def check_mask(u, v, mask):
"""
Set the velocity to 0 where the center of the correlation box lies
inside the mask. Enforces the no-slip condition and prevents the
median filter from bleeding values across the mask edge.
"""
u = u * mask
v = v * mask
return u, v
[docs]
def load_images(name_img_1, name_img_2):
"""
Load an image pair for PIV processing as grayscale float32 arrays.
A value of 1 is added so that no pixel outside the mask is
exactly 0.
"""
Img1 = 1 + np.asarray(cv2.cvtColor(cv2.imread(name_img_1),
cv2.COLOR_BGR2GRAY)).astype(np.float32)
Img2 = 1 + np.asarray(cv2.cvtColor(cv2.imread(name_img_2),
cv2.COLOR_BGR2GRAY)).astype(np.float32)
return Img1, Img2
[docs]
def save(x, y, u, v, filename, option='dpivsoft', Matlab=False, param=False):
"""
Save the flow field to disk, scaled with Parameters.calibration
and Parameters.delta_t.
Parameters
----------
x, y, u, v : 2d np.ndarray
Grid and velocity field to save.
filename : str
Output file name (extension added by NumPy/SciPy where needed).
option : {'dpivsoft', 'openpiv'}
'dpivsoft': save a python .npz file using the original
DPIVSoft (MATLAB) variable naming.
'openpiv': save the field in an ASCII file compatible with
OpenPIV.
Matlab : bool
If True, additionally save a MATLAB .mat file with all
processing parameters.
param : bool
If True, include the processing parameters in the .npz file.
"""
# Scale results
x = x * Parameters.calibration
y = y * Parameters.calibration
u = u * Parameters.calibration / Parameters.delta_t
v = v * Parameters.calibration / Parameters.delta_t
if Matlab:
mdic = {"x": x * 1.0, "y": y * 1.0, "u": u * 1.0, "v": v * 1.0,
"calibration": float(Parameters.calibration),
"delta_t": float(Parameters.delta_t),
"median_limit": float(Parameters.median_limit),
"no_calculation_1": float(Parameters.no_iter_1),
"no_calculation_2": float(Parameters.no_iter_2),
"box_size_1_x": float(Parameters.box_size_1_x),
"box_size_1_y": float(Parameters.box_size_1_y),
"box_size_2_x": float(Parameters.box_size_2_x),
"box_size_2_y": float(Parameters.box_size_2_y),
"no_boxes_1_x": float(Parameters.no_boxes_1_x),
"no_boxes_1_y": float(Parameters.no_boxes_1_y),
"no_boxes_2_x": float(Parameters.no_boxes_2_x),
"no_boxes_2_y": float(Parameters.no_boxes_2_y),
"no_calculation": float(Parameters.no_iter_1),
"direct_calculation": float(Parameters.direct_calc),
"gaussian_size": float(Parameters.gaussian_size),
"window_1_x": float(Parameters.window_1_x),
"window_1_y": float(Parameters.window_1_y),
"window_2_x": float(Parameters.window_2_x),
"window_2_y": float(Parameters.window_2_y),
"weighting": float(Parameters.weighting),
"peak_ratio": float(Parameters.peak_ratio),
"image_width": float(Parameters.Data.width),
"image_height": float(Parameters.Data.height),
"mask": float(Parameters.mask)}
savemat(filename + '.mat', mdic)
if option == 'dpivsoft':
if param:
np.savez(filename, x=x, y=y, u=u, v=v,
calibration=Parameters.calibration,
delta_t=Parameters.delta_t,
median_limit=Parameters.median_limit,
gaussian_size=Parameters.gaussian_size,
no_calculation_1=Parameters.no_iter_1,
no_calculation_2=Parameters.no_iter_2,
box_size_1_x=Parameters.box_size_1_x,
box_size_1_y=Parameters.box_size_1_y,
box_size_2_x=Parameters.box_size_2_x,
box_size_2_y=Parameters.box_size_2_y,
window_1_x=Parameters.window_1_x,
window_1_y=Parameters.window_1_y,
window_2_x=Parameters.window_2_x,
window_2_y=Parameters.window_2_y,
weighting=Parameters.weighting,
peak_ratio=Parameters.peak_ratio,
mask=Parameters.mask,
direct_calc=Parameters.direct_calc)
else:
np.savez(filename, x=x, y=y, u=u, v=v,
calibration=Parameters.calibration)
elif option == 'openpiv':
fmt = "%8.4f"
delimiter = "\t"
# Build output array
out = np.vstack([m.flatten() for m in [x, y, u, v, grid.mask_2]])
np.savetxt(
filename,
out.T,
fmt=fmt,
delimiter=delimiter,
header="x"
+ delimiter
+ "y"
+ delimiter
+ "u"
+ delimiter
+ "v"
+ delimiter
+ "mask",
)
else:
sys.exit("Saving option not found")