Uncertainty Quantification using Image Matching - Cylinder
Description
Based on paper and source code: - Sciacchitano, A., Wieneke, B., & Scarano, F. (2013). PIV uncertainty quantification by image matching. Measurement Science and Technology, 24 (4). https://doi.org/10.1088/0957-0233/24/4/045302. - http://piv.de/uncertainty/?page_id=221
Step 1: Particle peak detection
As described by (Sciacchitano et al., 2013, Eq. 1), the image intensity product \(\Pi\) from image matching intensities is defined as:
The peaks image \(\varphi\) is defined as:
Step 2: Disparity vector computation
The sub-pixel peak position estimator adopted here is the standard 3-point Gaussian fit. The particle positions of times \(t_1\) and \(t_2\) are defined as \(\boldsymbol{X}^1 = \{x^1_1,x^1_2,...,x^1_N\}\) and \(\boldsymbol{X}^2 = \{x^2_1,x^2_2,...,x^2_N\}\). Discrete disparity vectors are defined as:
Step 3: Disparity ensemble statistics inside window
The mean of disparity set inside a window is defined as:
where \(c_i = \sqrt{\Pi(x_i)}\) for \(i=1,2,...,N\).
The standard deviation of disparity set inside a window is defined as:
Finally, the instantaneous error (estimate) vector is defined as:
Setup
Packages
[1]:
%reload_ext autoreload
%autoreload 2
[2]:
import matplotlib.pyplot as plt
import numpy as np
import os
import pivuq
Load images
[3]:
parent_path = "./data/cylinder_wake/"
image_pair = np.array(
[plt.imread(os.path.join(parent_path + ipath)).astype("float") for ipath in ["frameA.tif", "frameB.tif"]]
)
Load reference velocity
[4]:
data = np.load(os.path.join(parent_path + "vectors_OF.npz"))
X_ref = data["X"]
Y_ref = data["Y"]
U_ref = data["U"]
Uncertainity quantificiation using image matching
Method: ILK
[5]:
%%time
X_ilk, Y_ilk, D = pivuq.disparity.ilk(
image_pair,
U_ref,
window_size=16,
window="gaussian",
velocity_upsample_kind="linear",
warp_direction="center",
warp_order=-1,
warp_nsteps=1,
)
CPU times: user 8.57 s, sys: 126 ms, total: 8.7 s
Wall time: 7.71 s
Method: SWS
[6]:
%%time
X_sws, Y_sws, delta, N, mu, sigma = pivuq.disparity.sws(
image_pair,
U_ref,
window_size=16,
window="gaussian",
radius=1,
sliding_window_subtraction=True,
ROI=None,
velocity_upsample_kind="linear",
warp_direction="center",
warp_order=-1,
warp_nsteps=1,
)
CPU times: user 15.1 s, sys: 18.6 ms, total: 15.1 s
Wall time: 4.38 s
Plot: Instantaneous error map
[7]:
fig, axes = plt.subplots(ncols=3, sharex=True, sharey=True, figsize=(15, 8))
# References velocity
ax = axes[0]
im = ax.contourf(X_ref, Y_ref, np.linalg.norm(U_ref, axis=0), np.linspace(0, 20, 11), extend="max")
fig.colorbar(im, ax=ax)
ax.set(title="$|U_{ref}|$")
# Uncertainty: ILK
ax = axes[1]
im = ax.contourf(X_ilk, Y_ilk, np.linalg.norm(D, axis=0), np.linspace(0, 1, 11), extend="max")
fig.colorbar(im, ax=ax)
ax.set(title="ILK method: $D$")
# Uncertainty: SWS
ax = axes[2]
im = ax.contourf(X_sws, Y_sws, np.linalg.norm(delta, axis=0), np.linspace(0, 1, 11), extend="max")
fig.colorbar(im, ax=ax)
ax.set(title="SWS method: $|\delta|$")
for ax in axes.ravel():
ax.set(xlim=(50, 400), ylim=(0, 600), xlabel="$j$", ylabel="$i$")